Merge pull request #49 from JesseLivezey/docs

JesseLivezey · web-flow · commit f1ffc7afbbbc · 2020-06-15T14:51:48.000-07:00
docs update
diff --git a/codecov.yaml b/codecov.yaml
@@ -3,3 +3,5 @@ ignore:
   - "tests"
   - "dca/plotting.py"
   - "dca/style.py"
+  - "dca/knn.py"
+  - "dca/data_util.py"
diff --git a/dca/__init__.py b/dca/__init__.py
@@ -2,11 +2,9 @@
                                  SlowFeatureAnalysis,
                                  ForecastableComponentsAnalysis)
 from .dca import (DynamicalComponentsAnalysis,
-                  DynamicalComponentsAnalysisKNN,
                   DynamicalComponentsAnalysisFFT)
 
 __all__ = ['DynamicalComponentsAnalysis',
-           'DynamicalComponentsAnalysisKNN',
            'DynamicalComponentsAnalysisFFT',
            'GaussianProcessFactorAnalysis',
            'SlowFeatureAnalysis',
diff --git a/dca/cov_util.py b/dca/cov_util.py
@@ -3,11 +3,54 @@
 import collections
 import torch
 import functools
+from numpy.lib.stride_tricks import as_strided
 
-from .data_util import form_lag_matrix
 from sklearn.utils.extmath import randomized_svd
 
 
+def form_lag_matrix(X, T, stride=1, stride_tricks=True, writeable=False):
+    """Form the data matrix with `T` lags.
+
+    Parameters
+    ----------
+    X : ndarray (n_time, N)
+        Timeseries with no lags.
+    T : int
+        Number of lags.
+    stride : int
+        Number of original samples to move between lagged samples.
+    stride_tricks : bool
+        Whether to use numpy stride tricks to form the lagged matrix or create
+        a new array. Using numpy stride tricks can can lower memory usage, especially for
+        large `T`. If `False`, a new array is created.
+    writeable : bool
+        For testing. You should not need to set this to True. This function uses stride tricks
+        to form the lag matrix which means writing to the array will have confusing behavior.
+        If `stride_tricks` is `False`, this flag does nothing.
+
+    Returns
+    -------
+    X_with_lags : ndarray (n_lagged_time, N * T)
+        Timeseries with lags.
+    """
+    if not isinstance(stride, int) or stride < 1:
+        raise ValueError('stride should be an int and greater than or equal to 1.')
+    N = X.shape[1]
+    n_lagged_samples = (len(X) - T) // stride + 1
+    if n_lagged_samples < 1:
+        raise ValueError('T is too long for a timeseries of length {}.'.format(len(X)))
+    if stride_tricks:
+        X = np.asarray(X, dtype=float, order='C')
+        shape = (n_lagged_samples, N * T)
+        strides = (X.strides[0] * stride,) + (X.strides[-1],)
+        X_with_lags = as_strided(X, shape=shape, strides=strides, writeable=writeable)
+    else:
+        X_with_lags = np.zeros((n_lagged_samples, T * N))
+        for i in range(n_lagged_samples):
+            X_with_lags[i, :] = X[i * stride:i * stride + T, :].flatten()
+    return X_with_lags
+
+
 def rectify_spectrum(cov, epsilon=1e-6, verbose=False):
     """Rectify the spectrum of a covariance matrix.
 
@@ -299,15 +342,15 @@ def calc_cov_from_cross_cov_mats(cross_cov_mats):
     return cov
 
 
-def calc_pi_from_data(X, T):
-    """Calculates the mutual information ("predictive information"
-    or "PI") between variables  {1,...,T_pi} and {T_pi+1,...,2*T_pi}, which
-    are jointly Gaussian with covariance matrix cov_2_T_pi.
+def calc_pi_from_data(X, T, proj=None):
+    """Calculates the Gaussian Predictive Information between variables
+    {1,...,T_pi} and {T_pi+1,...,2*T_pi}..
 
     Parameters
     ----------
-    cov_2_T_pi : np.ndarray, shape (2*T_pi, 2*T_pi)
-        Covariance matrix.
+    T : int
+        This T should be 2 * T_pi. This T sets the joint window length not the
+        past or future window length.
 
     Returns
     -------
@@ -316,13 +359,12 @@ def calc_pi_from_data(X, T):
     """
     ccms = calc_cross_cov_mats_from_data(X, T)
 
-    return calc_pi_from_cross_cov_mats(ccms)
+    return calc_pi_from_cross_cov_mats(ccms, proj=proj)
 
 
 def calc_pi_from_cov(cov_2_T_pi):
-    """Calculates the mutual information ("predictive information"
-    or "PI") between variables  {1,...,T_pi} and {T_pi+1,...,2*T_pi}, which
-    are jointly Gaussian with covariance matrix cov_2_T_pi.
+    """Calculates the Gaussian Predictive Information between variables
+    {1,...,T_pi} and {T_pi+1,...,2*T_pi} with covariance matrix cov_2_T_pi.
 
     Parameters
     ----------
diff --git a/dca/data_util.py b/dca/data_util.py
@@ -4,50 +4,8 @@
 from scipy.interpolate import interp1d
 from scipy.signal import resample
 from scipy.ndimage import convolve1d
-from numpy.lib.stride_tricks import as_strided
-
-
-def form_lag_matrix(X, T, stride=1, stride_tricks=True, writeable=False):
-    """Form the data matrix with `T` lags.
-
-    Parameters
-    ----------
-    X : ndarray (n_time, N)
-        Timeseries with no lags.
-    T : int
-        Number of lags.
-    stride : int
-        Number of original samples to move between lagged samples.
-    stride_tricks : bool
-        Whether to use numpy stride tricks to form the lagged matrix or create
-        a new array. Using numpy stride tricks can can lower memory usage, especially for
-        large `T`. If `False`, a new array is created.
-    writeable : bool
-        For testing. You should not need to set this to True. This function uses stride tricks
-        to form the lag matrix which means writing to the array will have confusing behavior.
-        If `stride_tricks` is `False`, this flag does nothing.
-
-    Returns
-    -------
-    X_with_lags : ndarray (n_lagged_time, N * T)
-        Timeseries with lags.
-    """
-    if not isinstance(stride, int) or stride < 1:
-        raise ValueError('stride should be an int and greater than or equal to 1.')
-    N = X.shape[1]
-    n_lagged_samples = (len(X) - T) // stride + 1
-    if n_lagged_samples < 1:
-        raise ValueError('T is too long for a timeseries of length {}.'.format(len(X)))
-    if stride_tricks:
-        X = np.asarray(X, dtype=float, order='C')
-        shape = (n_lagged_samples, N * T)
-        strides = (X.strides[0] * stride,) + (X.strides[-1],)
-        X_with_lags = as_strided(X, shape=shape, strides=strides, writeable=writeable)
-    else:
-        X_with_lags = np.zeros((n_lagged_samples, T * N))
-        for i in range(n_lagged_samples):
-            X_with_lags[i, :] = X[i * stride:i * stride + T, :].flatten()
-    return X_with_lags
+
+from .cov_util import form_lag_matrix  # noqa:F401
 
 
 def sum_over_chunks(X, stride):
diff --git a/dca/dca.py b/dca/dca.py
@@ -8,11 +8,10 @@
 import torch.nn.functional as F
 
 from .cov_util import (calc_cross_cov_mats_from_data, calc_pi_from_cross_cov_mats,
-                       form_lag_matrix, calc_pi_from_cross_cov_mats_block_toeplitz)
+                       calc_pi_from_cross_cov_mats_block_toeplitz)
 
 __all__ = ['DynamicalComponentsAnalysis',
            'DynamicalComponentsAnalysisFFT',
-           'DynamicalComponentsAnalysisKNN',
            'ortho_reg_fn',
            'build_loss',
            'init_coef']
@@ -181,8 +180,9 @@ class DynamicalComponentsAnalysis(object):
     """Dynamical Components Analysis.
 
     Runs DCA on multidimensional timeseries data X to discover a projection
-    onto a d-dimensional subspace which maximizes the complexity, as defined by the Gaussian
-    Predictive Information (PI) of the d-dimensional dynamics over windows of length T.
+    onto a d-dimensional subspace of an N-dimensional space which maximizes the complexity, as
+    defined by the Gaussian Predictive Information (PI) of the d-dimensional dynamics over windows
+    of length T.
 
     Parameters
     ----------
@@ -226,6 +226,18 @@ class DynamicalComponentsAnalysis(object):
 
     Attributes
     ----------
+    T : int
+        Default T used for PI.
+    T_fit : int
+        T used for last cross covariance estimation.
+    d : int
+        Default d used for fitting the projection.
+    d_fit : int
+        d used for last projection fit.
+    cross covs : torch tensor
+        Cross covariance matrices from the last covariance estimation.
+    coef_ : ndarray (N, d)
+        Projection matrix from fit.
 
     """
     def __init__(self, d=None, T=None, init="random_ortho", n_init=1, stride=1, tol=1e-6,
@@ -278,7 +290,7 @@ def estimate_cross_covariance(self, X, T=None, regularization=None, reg_ops=None
         X : ndarray or list of ndarrays
             Data to estimate the cross covariance matrix.
         T : int
-            T for PI calculation (optional.)
+            T for PI calculation (optional).
         regularization : str
             Whether to regularize cross covariance estimation.
         reg_ops : dict
@@ -313,6 +325,10 @@ def fit_projection(self, d=None, T=None, n_init=None):
         ----------
         d : int
             Dimensionality of the projection (optional.)
+        T : int
+            T for PI calculation (optional). Default is `self.T`. If `T` is set here
+            it must be less than or equal to `self.T` or self.estimate_cross_covariance() must
+            be called with a larger `T`.
         n_init : int
             Number of random restarts (optional.)
         """
@@ -340,6 +356,10 @@ def _fit_projection(self, d=None, T=None, record_V=False):
         ----------
         d : int
             Dimensionality of the projection (optional.)
+        T : int
+            T for PI calculation (optional). Default is `self.T`. If `T` is set here
+            it must be less than or equal to `self.T` or self.estimate_cross_covariance() must
+            be called with a larger `T`.
         record_V : bool
             If True, saves a copy of V at each optimization step. Default is False.
         """
@@ -350,13 +370,15 @@ def _fit_projection(self, d=None, T=None, record_V=False):
         self.d_fit = d
         if T is None:
             T = self.T
+        if T < 1:
+            raise ValueError
         if (2 * T) > self.cross_covs.shape[0]:
             raise ValueError('T must less than or equal to the value when ' +
-                             '`estimate_cross_covariance` was called.')
+                             '`estimate_cross_covariance()` was called.')
         self.T_fit = T
 
         if self.cross_covs is None:
-            raise ValueError('Call estimate_cross_covariance() first.')
+            raise ValueError('Call `estimate_cross_covariance()` first.')
 
         c = self.cross_covs[:2 * T]
         N = c.shape[1]
@@ -706,117 +728,3 @@ def score(self, X):
         X : ndarray or list
         """
         return pi_fft(X, self.coef_, self.T)
-
-
-class DynamicalComponentsAnalysisKNN(object):
-    """Dynamical Components Analysis using MI estimation using the
-    k-nearest neighbors methos of Kraskov, et al. This estimator is not
-    differentiable and so numerical gradients are taken (very slow!).
-
-    WARNING: This code has not been used or tested and is still being developed.
-
-    Runs DCA on multidimensional timeseries data X to discover a projection
-    onto a d-dimensional subspace which maximizes the complexity of the d-dimensional
-    dynamics.
-
-    Parameters
-    ----------
-    d : int
-        Number of basis vectors onto which the data X are projected.
-    T : int
-        Size of time windows across which to compute mutual information.
-    init : string
-        Options: "random", "PCA"
-        Method for initializing the projection matrix.
-    """
-    def __init__(self, d=None, T=None, init="random", n_init=1, tol=1e-6,
-                 ortho_lambda=10., verbose=False, use_scipy=True):
-        self.d = d
-        self.T = T
-        self.init = init
-        self.n_init = n_init
-        self.tol = tol
-        self.ortho_lambda = ortho_lambda
-        self.verbose = verbose
-        self.coef_ = None
-
-    def fit(self, X, d=None, T=None, n_init=None):
-        self.mean_ = X.mean(axis=0, keepdims=True)
-        X -= self.mean_
-        if n_init is None:
-            n_init = self.n_init
-        pis = []
-        coefs = []
-        for ii in range(n_init):
-            coef, pi = self._fit_projection(X, d=d)
-            pis.append(pi)
-            coefs.append(coef)
-        idx = np.argmax(pis)
-        self.coef_ = coefs[idx]
-        return self
-
-    def _fit_projection(self, X, d=None):
-        from info_measures.continuous import kraskov_stoegbauer_grassberger as ksg
-        if d is None:
-            d = self.d
-
-        N = X.shape[1]
-        if type(self.init) == str:
-            if self.init == "random":
-                V_init = np.random.normal(0, 1, (N, d))
-            elif self.init == "random_ortho":
-                V_init = scipy.stats.ortho_group.rvs(N)[:, :d]
-            elif self.init == "uniform":
-                V_init = np.ones((N, d)) / np.sqrt(N)
-                V_init = V_init + np.random.normal(0, 1e-3, V_init.shape)
-            else:
-                raise ValueError
-        else:
-            raise ValueError
-        V_init /= np.linalg.norm(V_init, axis=0, keepdims=True)
-
-        callback = None
-        if self.verbose:
-
-            def callback(v_flat):
-                v = v_flat.reshape(N, d)
-                X_lag = form_lag_matrix(X.dot(v), 2 * self.T)
-                mi = ksg.MutualInformation(X_lag[:, :self.T * d],
-                                           X_lag[:, self.T * d:])
-                pi = mi.mutual_information()
-                reg_val = ortho_reg_fn(v, self.ortho_lambda)
-                print("PI: {} bits, reg: {}".format(str(np.round(pi, 4)),
-                                                    str(np.round(reg_val, 4))))
-            callback(V_init)
-
-        def f(v_flat):
-            v = v_flat.reshape(N, d)
-            X_lag = form_lag_matrix(X.dot(v), 2 * self.T)
-            mi = ksg.MutualInformation(X_lag[:, :self.T], X_lag[:, self.T:])
-            pi = mi.mutual_information()
-            reg_val = ortho_reg_fn(v, self.ortho_lambda)
-            loss = -pi + reg_val
-            return loss
-        opt = minimize(f, V_init.ravel(), method='L-BFGS-B', callback=callback)
-        v = opt.x.reshape(N, d)
-
-        # Orthonormalize the basis prior to returning it
-        V_opt = scipy.linalg.orth(v)
-        final_pi = self.score(X, V_opt)
-        return V_opt, final_pi
-
-    def transform(self, X):
-        return (X - self.mean_).dot(self.coef_)
-
-    def fit_transform(self, X, d=None, T=None):
-        self.fit(X, d=d, T=T)
-        return self.transform(X)
-
-    def score(self, X, coef=None):
-        if coef is None:
-            coef = self.coef_
-        from info_measures.continuous import kraskov_stoegbauer_grassberger as ksg
-        X_lag = form_lag_matrix(X.dot(coef), 2 * self.T)
-        mi = ksg.MutualInformation(X_lag[:, :self.T], X_lag[:, self.T:])
-        pi = mi.mutual_information()
-        return pi
diff --git a/dca/knn.py b/dca/knn.py
diff --git a/tests/test_dca.py b/tests/test_dca.py
