from collections import namedtuple

import numpy as np

from . import parallel
from . import util
from .py2to3 import range
from .vectorspace import VectorSpaceArrays, VectorSpaceHandles


def compute_DMD_arrays_snaps_method(
    vecs, adv_vecs=None, mode_indices=None, inner_product_weights=None,
    atol=1e-13, rtol=None, max_num_eigvals=None):
    """Computes DMD modes using data stored in arrays, using method of
    snapshots.

    Args:
        ``vecs``: Array whose columns are data vectors.

    Kwargs:
        ``adv_vecs``: Array whose columns are data vectors advanced in time.
        If not provided, then it is assumed that the vectors describe a
        sequential time-series. Thus ``vecs`` becomes ``vecs[:, :-1]`` and
        ``adv_vecs`` becomes ``vecs[:, 1:]``.

        ``mode_indices``: List of indices describing which modes to compute.
        Examples are ``range(10)`` or ``[3, 0, 6, 8]``.  If no mode indices are
        specified, then all modes will be computed.

        ``inner_product_weights``: 1D or 2D array of inner product weights.
        Corresponds to :math:`W` in inner product :math:`v_1^* W v_2`.

        ``atol``: Level below which eigenvalues of correlation array are
        truncated.

        ``rtol``: Maximum relative difference between largest and smallest
        eigenvalues of correlation array.  Smaller ones are truncated.

        ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
        computed.  This is enforced by truncating the basis onto which the
        approximating linear map is projected.  Computationally, this
        corresponds to truncating the eigendecomposition of the correlation
        array. If set to None, no truncation will be performed, and the
        maximum possible number of DMD eigenvalues will be computed.

    Returns:
        ``res``: Results of DMD computation, stored in a namedtuple with
        the following attributes:

        * ``eigvals``: 1D array of eigenvalues of approximating low-order linear
          map (DMD eigenvalues).

        * ``spectral_coeffs``: 1D array of DMD spectral coefficients, calculated
          as the magnitudes of the projection coefficients of first data vector.
          The projection is onto the span of the DMD modes using the
          (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``exact_modes``: Array whose columns are exact DMD modes.

        * ``proj_modes``: Array whose columns are projected DMD modes.

        * ``adjoint_modes``: Array whose columns are adjoint DMD modes.

        * ``proj_coeffs``: Array of projection coefficients for data vectors
          expressed as a linear combination of DMD modes.  Columns correspond to
          vector objects, rows correspond to DMD modes. The projection is onto
          the span of the DMD modes using the (biorthogonal) adjoint DMD modes.
          Note that this is the same as a least-squares projection onto the span
          of the DMD modes.

        * ``adv_proj_coeffs``: Array of projection coefficients for data vectors
          advanced in time, expressed as a linear combination of DMD modes.
          Columns correspond to vector objects, rows correspond to DMD
          modes. The projection is onto the span of the DMD modes using the
          (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``R_low_order_eigvecs``: Array of right eigenvectors of approximating
          low-order linear map.

        * ``L_low_order_eigvecs``: Array of left eigenvectors of approximating
          low-order linear map.

        * ``correlation_array_eigvals``: 1D array of eigenvalues of
          correlation array.

        * ``correlation_array_eigvecs``: Array of eigenvectors of
          correlation array.

        * ``correlation_array``: Correlation array; elements are inner products
          of data vectors with each other.

        * ``cross_correlation_array``: Cross-correlation array; elements are
          inner products of data vectors with data vectors advanced in time.
          Going down rows, the data vector changes; going across columns the
          advanced data vector changes.

        Attributes can be accessed using calls like ``res.exact_modes``.  To
        see all available attributes, use ``print(res)``.

    This uses the method of snapshots, which is faster than the direct method
    (see :py:func:`compute_DMD_arrays_direct_method`) when ``vecs`` has more
    rows than columns, i.e., when there are more elements in a vector than
    there are vectors. However, it "squares" this array and its singular
    values, making it slightly less accurate than the direct method.

    """
    if parallel.is_distributed():
        raise RuntimeError('Cannot run in parallel.')

    # Force data to be arrays (not matrices)
    vecs = np.array(vecs)
    if adv_vecs is not None:
        adv_vecs = np.array(adv_vecs)

    # Set up vector space (for inner products)
    vec_space = VectorSpaceArrays(weights=inner_product_weights)

    # Sequential dataset
    if adv_vecs is None:
        # Compute correlation array for all vectors.
        # This is more efficient because only one call is made to the inner
        # product routine, even though we don't need the last row and column
        # yet.
        expanded_correlation_array =\
            vec_space.compute_symm_inner_product_array(vecs)
        correlation_array = expanded_correlation_array[:-1, :-1]
        cross_correlation_array = expanded_correlation_array[:-1, 1:]
    # Non-sequential data
    else:
        # Check that dimensions match
        if vecs.shape != adv_vecs.shape:
            raise ValueError(('vecs and adv_vecs are not the same shape.'))

        # Compute the correlation array from the unadvanced snapshots only.
        correlation_array = vec_space.compute_symm_inner_product_array(
            vecs)
        cross_correlation_array = vec_space.compute_inner_product_array(
            vecs, adv_vecs)

    correlation_array_eigvals, correlation_array_eigvecs = util.eigh(
        correlation_array, is_positive_definite=True, atol=atol, rtol=rtol)

    # Truncate if necessary
    if max_num_eigvals is not None and (
        max_num_eigvals < correlation_array_eigvals.size):
        correlation_array_eigvals = correlation_array_eigvals[:max_num_eigvals]
        correlation_array_eigvecs = correlation_array_eigvecs[
            :, :max_num_eigvals]

    # Compute low-order linear map for sequential or non-sequential case.
    correlation_array_eigvals_sqrt_inv = np.diag(
        correlation_array_eigvals ** -0.5)
    low_order_linear_map = correlation_array_eigvals_sqrt_inv.dot(
        correlation_array_eigvecs.conj().T.dot(
            cross_correlation_array.dot(
                correlation_array_eigvecs.dot(
                    correlation_array_eigvals_sqrt_inv))))

    # Compute eigendecomposition of low-order linear map.
    eigvals, R_low_order_eigvecs, L_low_order_eigvecs =\
        util.eig_biorthog(low_order_linear_map, scale_choice='left')

    # Compute build coefficients
    build_coeffs_proj = correlation_array_eigvecs.dot(
        correlation_array_eigvals_sqrt_inv.dot(
            R_low_order_eigvecs))
    build_coeffs_exact = build_coeffs_proj.dot(np.diag(eigvals ** -1.))
    build_coeffs_adjoint = correlation_array_eigvecs.dot(
        correlation_array_eigvals_sqrt_inv.dot(
            L_low_order_eigvecs))

    # Compute spectral coefficients
    spectral_coeffs = np.abs(L_low_order_eigvecs.conj().T.dot(
        np.diag(correlation_array_eigvals ** 0.5).dot(
            correlation_array_eigvecs[0, :].conj().T)).squeeze())

    # Compute projection coefficients
    proj_coeffs = L_low_order_eigvecs.conj().T.dot(
        np.diag(correlation_array_eigvals ** 0.5).dot(
            correlation_array_eigvecs.conj().T))
    adv_proj_coeffs = L_low_order_eigvecs.conj().T.dot(
        np.diag(correlation_array_eigvals ** -0.5).dot(
            correlation_array_eigvecs.conj().T.dot(
                cross_correlation_array)))

    # For sequential data, user must provide one more vec than columns of
    # build_coeffs.
    if vecs.shape[1] - build_coeffs_exact.shape[0] == 1:
        exact_modes = vec_space.lin_combine(
            vecs[:, 1:], build_coeffs_exact,
            coeff_array_col_indices=mode_indices)
        proj_modes = vec_space.lin_combine(
            vecs[:, :-1], build_coeffs_proj,
            coeff_array_col_indices=mode_indices)
        adjoint_modes = vec_space.lin_combine(
            vecs[:, :-1], build_coeffs_adjoint,
            coeff_array_col_indices=mode_indices)
    # For non-sequential data, user must provide as many vecs as columns of
    # build_coeffs.
    elif vecs.shape[1] == build_coeffs_exact.shape[0]:
        exact_modes = vec_space.lin_combine(
            adv_vecs, build_coeffs_exact, coeff_array_col_indices=mode_indices)
        proj_modes = vec_space.lin_combine(
            vecs, build_coeffs_proj, coeff_array_col_indices=mode_indices)
        adjoint_modes = vec_space.lin_combine(
            vecs, build_coeffs_adjoint, coeff_array_col_indices=mode_indices)
    else:
        raise ValueError(
            'Number of cols in vecs does not match number of rows in '
            'build_coeffs array.')

    # Return a namedtuple
    DMD_results = namedtuple(
        'DMD_results', [
            'eigvals', 'spectral_coeffs',
            'exact_modes', 'proj_modes', 'adjoint_modes',
            'proj_coeffs', 'adv_proj_coeffs',
            'R_low_order_eigvecs', 'L_low_order_eigvecs',
            'correlation_array_eigvals', 'correlation_array_eigvecs',
            'correlation_array', 'cross_correlation_array'])
    return DMD_results(
        eigvals=eigvals, spectral_coeffs=spectral_coeffs,
        exact_modes=exact_modes,
        proj_modes=proj_modes, adjoint_modes=adjoint_modes,
        proj_coeffs=proj_coeffs, adv_proj_coeffs=adv_proj_coeffs,
        R_low_order_eigvecs=R_low_order_eigvecs,
        L_low_order_eigvecs=L_low_order_eigvecs,
        correlation_array_eigvals=correlation_array_eigvals,
        correlation_array_eigvecs=correlation_array_eigvecs,
        correlation_array=correlation_array,
        cross_correlation_array=cross_correlation_array)


def compute_DMD_arrays_direct_method(
    vecs, adv_vecs=None, mode_indices=None, inner_product_weights=None,
    atol=1e-13, rtol=None, max_num_eigvals=None):
    """Computes DMD modes using data stored in arrays, using direct method.

    Args:
        ``vecs``: Array whose columns are data vectors.

    Kwargs:
        ``adv_vecs``: Array whose columns are data vectors advanced in time.
        If not provided, then it is assumed that the vectors describe a
        sequential time-series. Thus ``vecs`` becomes ``vecs[:, :-1]`` and
        ``adv_vecs`` becomes ``vecs[:, 1:]``.

        ``mode_indices``: List of indices describing which modes to compute.
        Examples are ``range(10)`` or ``[3, 0, 6, 8]``.  If no mode indices are
        specified, then all modes will be computed.

        ``inner_product_weights``: 1D or 2D array of inner product weights.
        Corresponds to :math:`W` in inner product :math:`v_1^* W v_2`.

        ``atol``: Level below which eigenvalues of correlation array are
        truncated.

        ``rtol``: Maximum relative difference between largest and smallest
        eigenvalues of correlation array.  Smaller ones are truncated.

        ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
        computed.  This is enforced by truncating the basis onto which the
        approximating linear map is projected.  Computationally, this
        corresponds to truncating the eigendecomposition of the correlation
        array. If set to None, no truncation will be performed, and the
        maximum possible number of DMD eigenvalues will be computed.

    Returns:
        ``res``: Results of DMD computation, stored in a namedtuple with
        the following attributes:

        * ``eigvals``: 1D array of eigenvalues of approximating low-order linear
          map (DMD eigenvalues).

        * ``spectral_coeffs``: 1D array of DMD spectral coefficients, calculated
          as the magnitudes of the projection coefficients of first data vector.
          The projection is onto the span of the DMD modes using the
          (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``exact_modes``: Array whose columns are exact DMD modes.

        * ``proj_modes``: Array whose columns are projected DMD modes.

        * ``adjoint_modes``: Array whose columns are adjoint DMD modes.

        * ``proj_coeffs``: Array of projection coefficients for data vectors
          expressed as a linear combination of DMD modes.  Columns correspond to
          vector objects, rows correspond to DMD modes. The projection is onto
          the span of the DMD modes using the (biorthogonal) adjoint DMD modes.
          Note that this is the same as a least-squares projection onto the span
          of the DMD modes.

        * ``adv_proj_coeffs``: Array of projection coefficients for data vectors
          advanced in time, expressed as a linear combination of DMD modes.
          Columns correspond to vector objects, rows correspond to DMD
          modes. The projection is onto the span of the DMD modes using the
          (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``R_low_order_eigvecs``: Array of right eigenvectors of approximating
          low-order linear map.

        * ``L_low_order_eigvecs``: Array of left eigenvectors of approximating
          low-order linear map.

        * ``correlation_array_eigvals``: 1D array of eigenvalues of
          correlation array.

        * ``correlation_array_eigvecs``: Array of eigenvectors of
          correlation array.

        Attributes can be accessed using calls like ``res.exact_modes``.  To
        see all available attributes, use ``print(res)``.

    This method does not square the array of vectors as in the method of
    snapshots (:py:func:`compute_DMD_arrays_snaps_method`). It's slightly
    more accurate, but slower when the number of elements in a vector is more
    than the number of vectors, i.e.,  when ``vecs`` has more rows than
    columns.
    """
    if parallel.is_distributed():
        raise RuntimeError('Cannot run in parallel.')

    # Force data to be arrays (not matrices)
    vecs = np.array(vecs)
    if adv_vecs is not None:
        adv_vecs = np.array(adv_vecs)

    # Set up vector space (for inner products)
    vec_space = VectorSpaceArrays(weights=inner_product_weights)

    # Weight data as necessary
    if inner_product_weights is None:
        vecs_weighted = vecs
        if adv_vecs is not None:
            adv_vecs_weighted = adv_vecs
    elif inner_product_weights.ndim == 1:
        sqrt_weights = np.diag(inner_product_weights ** 0.5)
        vecs_weighted = sqrt_weights.dot(vecs)
        if adv_vecs is not None:
            adv_vecs_weighted = sqrt_weights.dot(adv_vecs)
    elif inner_product_weights.ndim == 2:
        if inner_product_weights.shape[0] > 500:
            print('Warning: Cholesky decomposition could be time consuming.')
        sqrt_weights = np.linalg.cholesky(inner_product_weights).conj().T
        vecs_weighted = sqrt_weights.dot(vecs)
        if adv_vecs is not None:
            adv_vecs_weighted = sqrt_weights.dot(adv_vecs)

    # Compute low-order linear map for sequential snapshot set.  This takes
    # advantage of the fact that for a sequential dataset, the unadvanced
    # and advanced vectors overlap.
    if adv_vecs is None:
        U, sing_vals, correlation_array_eigvecs = util.svd(
            vecs_weighted[:, :-1], atol=atol, rtol=rtol)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < sing_vals.size):
            U = U[:, :max_num_eigvals]
            sing_vals = sing_vals[:max_num_eigvals]
            correlation_array_eigvecs = correlation_array_eigvecs[
                :, :max_num_eigvals]

        # Compute correlation arrays.  Note that these include the truncation of
        # the singular vectors, so they are not equal to the full correlation
        # arrays.  However, they are useful for computations involving the
        # singular vectors.
        correlation_array_eigvals = sing_vals ** 2.
        correlation_array_eigvals_sqrt_inv = np.diag(sing_vals ** -1.)
        correlation_array = correlation_array_eigvecs.dot(
            np.diag(correlation_array_eigvals).dot(
                correlation_array_eigvecs.conj().T))
        cross_correlation_array = np.hstack((
            correlation_array[:, 1:],
            util.atleast_2d_col(
                vecs_weighted[:, :-1].conj().T.dot(vecs_weighted[:, -1]))))

    else:
        if vecs.shape != adv_vecs.shape:
            raise ValueError(('vecs and adv_vecs are not the same shape.'))
        U, sing_vals, correlation_array_eigvecs = util.svd(
            vecs_weighted, atol=atol, rtol=rtol)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < sing_vals.size):
            U = U[:, :max_num_eigvals]
            sing_vals = sing_vals[:max_num_eigvals]
            correlation_array_eigvecs = correlation_array_eigvecs[
                :, :max_num_eigvals]

        # Compute correlation arrays.  Note that these include the truncation of
        # the singular vectors, so they are not equal to the full correlation
        # arrays.  However, they are useful for computations involving the
        # singular vectors.
        correlation_array_eigvals = sing_vals ** 2.
        correlation_array_eigvals_sqrt_inv = np.diag(sing_vals ** -1.)
        correlation_array = correlation_array_eigvecs.dot(
            np.diag(correlation_array_eigvals).dot(
                correlation_array_eigvecs.conj().T))
        cross_correlation_array = vecs_weighted.conj().T.dot(adv_vecs_weighted)

    # Compute low-order lienar map
    low_order_linear_map = correlation_array_eigvals_sqrt_inv.dot(
        correlation_array_eigvecs.conj().T.dot(
            cross_correlation_array.dot(
                correlation_array_eigvecs.dot(
                    correlation_array_eigvals_sqrt_inv))))

    # Compute eigendecomposition of low-order linear map
    eigvals, R_low_order_eigvecs, L_low_order_eigvecs =\
        util.eig_biorthog(low_order_linear_map, scale_choice='left')

    # Compute build coefficients
    build_coeffs_proj = correlation_array_eigvecs.dot(
        correlation_array_eigvals_sqrt_inv.dot(
            R_low_order_eigvecs))
    build_coeffs_exact = build_coeffs_proj.dot(np.diag(eigvals ** -1.))
    build_coeffs_adjoint = correlation_array_eigvecs.dot(
        correlation_array_eigvals_sqrt_inv.dot(
            L_low_order_eigvecs))

    # Compute spectral coefficients
    spectral_coeffs = np.abs(L_low_order_eigvecs.conj().T.dot(
        np.diag(correlation_array_eigvals ** 0.5).dot(
            correlation_array_eigvecs[0, :].conj().T))).squeeze()

    # Compute projection coefficients
    proj_coeffs = L_low_order_eigvecs.conj().T.dot(
        np.diag(correlation_array_eigvals ** 0.5).dot(
            correlation_array_eigvecs.conj().T))
    adv_proj_coeffs = L_low_order_eigvecs.conj().T.dot(
        np.diag(correlation_array_eigvals ** -0.5).dot(
            correlation_array_eigvecs.conj().T.dot(
                cross_correlation_array)))

    # For sequential data, user must provide one more vec than columns of
    # build_coeffs.
    if vecs.shape[1] - build_coeffs_exact.shape[0] == 1:
        exact_modes = vec_space.lin_combine(
            vecs[:, 1:], build_coeffs_exact,
            coeff_array_col_indices=mode_indices)
        proj_modes = vec_space.lin_combine(
            vecs[:, :-1], build_coeffs_proj,
            coeff_array_col_indices=mode_indices)
        adjoint_modes = vec_space.lin_combine(
            vecs[:, :-1], build_coeffs_adjoint,
            coeff_array_col_indices=mode_indices)
    # For sequential data, user must provide as many vecs as columns of
    # build_coeffs.
    elif vecs.shape[1] == build_coeffs_exact.shape[0]:
        exact_modes = vec_space.lin_combine(
            adv_vecs, build_coeffs_exact, coeff_array_col_indices=mode_indices)
        proj_modes = vec_space.lin_combine(
            vecs, build_coeffs_proj, coeff_array_col_indices=mode_indices)
        adjoint_modes = vec_space.lin_combine(
            vecs, build_coeffs_adjoint, coeff_array_col_indices=mode_indices)
    else:
        raise ValueError(('Number of cols in vecs does not match '
            'number of rows in build_coeffs array.'))

    # Return a namedtuple
    DMD_results = namedtuple(
        'DMD_results', [
            'exact_modes', 'proj_modes', 'adjoint_modes',
            'spectral_coeffs', 'proj_coeffs', 'adv_proj_coeffs',
            'eigvals', 'R_low_order_eigvecs', 'L_low_order_eigvecs',
            'correlation_array_eigvals', 'correlation_array_eigvecs'])
    return DMD_results(
        exact_modes=exact_modes,
        proj_modes=proj_modes, adjoint_modes=adjoint_modes,
        spectral_coeffs=spectral_coeffs,
        proj_coeffs=proj_coeffs, adv_proj_coeffs=adv_proj_coeffs,
        eigvals=eigvals, R_low_order_eigvecs=R_low_order_eigvecs,
        L_low_order_eigvecs=L_low_order_eigvecs,
        correlation_array_eigvals=correlation_array_eigvals,
        correlation_array_eigvecs=correlation_array_eigvecs)


class DMDHandles(object):
    """Dynamic Mode Decomposition implemented for large datasets.

    Kwargs:
        ``inner_product``: Function that computes inner product of two vector
        objects.

        ``put_array``: Function to put an array out of modred, e.g., write it to
        file.

        ``get_array``: Function to get an array into modred, e.g., load it from
        file.

        ``max_vecs_per_node``: Maximum number of vectors that can be stored in
        memory, per node.

        ``verbosity``: 1 prints progress and warnings, 0 prints almost nothing.

    Computes DMD modes from vector objects (or handles).  It uses
    :py:class:`vectorspace.VectorSpaceHandles` for low level functions.

    Usage::

      myDMD = DMDHandles(inner_product=my_inner_product)
      myDMD.compute_decomp(vec_handles)
      myDMD.compute_exact_modes(range(50), mode_handles)

    See also :func:`compute_DMD_arrays_snaps_method`,
    :func:`compute_DMD_arrays_direct_method`, and :mod:`vectors`.
    """
    def __init__(
        self, inner_product=None, get_array=util.load_array_text,
        put_array=util.save_array_text, max_vecs_per_node=None, verbosity=1):
        """Constructor"""
        self.get_array = get_array
        self.put_array = put_array
        self.verbosity = verbosity
        self.eigvals = None
        self.correlation_array = None
        self.cross_correlation_array = None
        self.correlation_array_eigvals = None
        self.correlation_array_eigvecs = None
        self.low_order_linear_map = None
        self.L_low_order_eigvecs = None
        self.R_low_order_eigvecs = None
        self.spectral_coeffs = None
        self.proj_coeffs = None
        self.adv_proj_coeffs = None
        self.vec_space = VectorSpaceHandles(
            inner_product=inner_product, max_vecs_per_node=max_vecs_per_node,
            verbosity=verbosity)
        self.vec_handles = None
        self.adv_vec_handles = None


    def get_decomp(
        self, eigvals_src, R_low_order_eigvecs_src, L_low_order_eigvecs_src,
        correlation_array_eigvals_src, correlation_array_eigvecs_src):
        """Gets the decomposition arrays from sources (memory or file).

        Args:
            ``eigvals_src``: Source from which to retrieve eigenvalues of
            approximating low-order linear map (DMD eigenvalues).

            ``R_low_order_eigvecs_src``: Source from which to retrieve right
            eigenvectors of approximating low-order linear DMD map.

            ``L_low_order_eigvecs_src``: Source from which to retrieve left
            eigenvectors of approximating low-order linear DMD map.

            ``correlation_array_eigvals_src``: Source from which to retrieve
            eigenvalues of correlation array.

            ``correlation_array_eigvecs_src``: Source from which to retrieve
            eigenvectors of correlation array.
        """
        self.eigvals = np.squeeze(np.array(
            parallel.call_and_bcast(self.get_array, eigvals_src)))
        self.R_low_order_eigvecs = parallel.call_and_bcast(
            self.get_array, R_low_order_eigvecs_src)
        self.L_low_order_eigvecs = parallel.call_and_bcast(
            self.get_array, L_low_order_eigvecs_src)
        self.correlation_array_eigvals = np.squeeze(np.array(
            parallel.call_and_bcast(
            self.get_array, correlation_array_eigvals_src)))
        self.correlation_array_eigvecs = parallel.call_and_bcast(
            self.get_array, correlation_array_eigvecs_src)


    def get_correlation_array(self, src):
      """Gets the correlation array from source (memory or file).

        Args:
            ``src``: Source from which to retrieve correlation array.
        """
      self.correlation_array = parallel.call_and_bcast(self.get_array, src)


    def get_cross_correlation_array(self, src):
        """Gets the cross-correlation array from source (memory or file).

        Args:
            ``src``: Source from which to retrieve cross-correlation array.
        """
        self.cross_correlation_array = parallel.call_and_bcast(
            self.get_array, src)


    def get_spectral_coeffs(self, src):
        """Gets the spectral coefficients from source (memory or file).

        Args:
            ``src``: Source from which to retrieve spectral coefficients.
        """
        self.spectral_coeffs = parallel.call_and_bcast(self.get_array, src)


    def get_proj_coeffs(self, src, adv_src):
        """Gets the projection coefficients and advanced projection coefficients
        from sources (memory or file).

        Args:
            ``src``: Source from which to retrieve projection coefficients.

            ``adv_src``: Source from which to retrieve advanced projection
            coefficients.
        """
        self.proj_coeffs = parallel.call_and_bcast(self.get_array, src)
        self.adv_proj_coeffs = parallel.call_and_bcast(self.get_array, adv_src)


    def put_decomp(
        self, eigvals_dest, R_low_order_eigvecs_dest, L_low_order_eigvecs_dest,
        correlation_array_eigvals_dest, correlation_array_eigvecs_dest):
        """Puts the decomposition arrays in destinations (memory or file).

        Args:
            ``eigvals_dest``: Destination in which to put eigenvalues of
            approximating low-order linear map (DMD eigenvalues).

            ``R_low_order_eigvecs_dest``: Destination in which to put right
            eigenvectors of approximating low-order linear map.

            ``L_low_order_eigvecs_dest``: Destination in which to put left
            eigenvectors of approximating low-order linear map.

            ``correlation_array_eigvals_dest``: Destination in which to put
            eigenvalues of correlation array.

            ``correlation_array_eigvecs_dest``: Destination in which to put
            eigenvectors of correlation array.
        """
        # Don't check if rank is zero because the following methods do.
        self.put_eigvals(eigvals_dest)
        self.put_R_low_order_eigvecs(R_low_order_eigvecs_dest)
        self.put_L_low_order_eigvecs(L_low_order_eigvecs_dest)
        self.put_correlation_array_eigvals(correlation_array_eigvals_dest)
        self.put_correlation_array_eigvecs(correlation_array_eigvecs_dest)


    def put_eigvals(self, dest):
        """Puts eigenvalues of approximating low-order-linear map (DMD
        eigenvalues) to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.eigvals, dest)
        parallel.barrier()


    def put_R_low_order_eigvecs(self, dest):
        """Puts right eigenvectors of approximating low-order linear map to
        ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.R_low_order_eigvecs, dest)
        parallel.barrier()


    def put_L_low_order_eigvecs(self, dest):
        """Puts left eigenvectors of approximating low-order linear map to
        ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.L_low_order_eigvecs, dest)
        parallel.barrier()


    def put_correlation_array_eigvals(self, dest):
        """Puts eigenvalues of correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.correlation_array_eigvals, dest)
        parallel.barrier()


    def put_correlation_array_eigvecs(self, dest):
        """Puts eigenvectors of correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.correlation_array_eigvecs, dest)
        parallel.barrier()


    def put_correlation_array(self, dest):
        """Puts correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.correlation_array, dest)
        parallel.barrier()


    def put_cross_correlation_array(self, dest):
        """Puts cross-correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.cross_correlation_array, dest)
        parallel.barrier()


    def put_spectral_coeffs(self, dest):
        """Puts DMD spectral coefficients to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.spectral_coeffs, dest)
        parallel.barrier()


    def put_proj_coeffs(self, dest, adv_dest):
        """Puts projection coefficients to ``dest``, advanced projection
        coefficients to ``adv_dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.proj_coeffs, dest)
            self.put_array(self.adv_proj_coeffs, adv_dest)
        parallel.barrier()


    def sanity_check(self, test_vec_handle):
        """Checks that user-supplied vector handle and vector satisfy
        requirements.

        Args:
            ``test_vec_handle``: A vector handle to test.

        See :py:meth:`vectorspace.VectorSpaceHandles.sanity_check`.
        """
        self.vec_space.sanity_check(test_vec_handle)


    def compute_eigendecomp(self, atol=1e-13, rtol=None, max_num_eigvals=None):
        """Computes eigendecompositions of correlation array and approximating
        low-order linear map.

        Kwargs:
            ``atol``: Level below which eigenvalues of correlation array are
            truncated.

            ``rtol``: Maximum relative difference between largest and smallest
            eigenvalues of correlation array.  Smaller ones are truncated.

            ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
            computed.  This is enforced by truncating the basis onto which the
            approximating linear map is projected.  Computationally, this
            corresponds to truncating the eigendecomposition of the correlation
            array. If set to None, no truncation will be performed, and the
            maximum possible number of DMD eigenvalues will be computed.

        Useful if you already have the correlation array and cross-correlation
        array and want to avoid recomputing them.

        Usage::

          DMD.correlation_array = pre_existing_correlation_array
          DMD.cross_correlation_array = pre_existing_cross_correlation_array
          DMD.compute_eigendecomp()
          DMD.compute_exact_modes(
              mode_idx_list, mode_handles, adv_vec_handles=adv_vec_handles)

        Another way to use this is to compute a DMD using a truncated basis for
        the projection of the approximating linear map.  Start by either
        computing a full decomposition or by loading pre-computed correlation
        and cross-correlation arrays.

        Usage::

          # Start with a full decomposition
          DMD_eigvals, correlation_array_eigvals = DMD.compute_decomp(
              vec_handles)[0, 3]

          # Do some processing to determine the truncation level, maybe based
          # on the DMD eigenvalues and correlation array eigenvalues
          desired_num_eigvals = my_post_processing_func(
              DMD_eigvals, correlation_array_eigvals)

          # Do a truncated decomposition
          DMD_eigvals_trunc = DMD.compute_eigendecomp(
            max_num_eigvals=desired_num_eigvals)

          # Compute modes for truncated decomposition
          DMD.compute_exact_modes(
              mode_idx_list, mode_handles, adv_vec_handles=adv_vec_handles)

        Since it doesn't overwrite the correlation and cross-correlation
        arrays, ``compute_eigendecomp`` can be called many times in a row to
        do computations for different truncation levels.  However, the results
        of the decomposition (e.g., ``self.eigvals``) do get overwritten, so
        you may want to call a ``put`` method to save those results somehow.
        """
        # Compute eigendecomposition of correlation array
        self.correlation_array_eigvals, self.correlation_array_eigvecs =\
            parallel.call_and_bcast(
            util.eigh, self.correlation_array, atol=atol, rtol=None,
            is_positive_definite=True)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < self.correlation_array_eigvals.size):
            self.correlation_array_eigvals = self.correlation_array_eigvals[
                :max_num_eigvals]
            self.correlation_array_eigvecs = self.correlation_array_eigvecs[
                :, :max_num_eigvals]

        # Compute low-order linear map
        correlation_array_eigvals_sqrt_inv = np.diag(
            self.correlation_array_eigvals ** -0.5)
        self.low_order_linear_map = correlation_array_eigvals_sqrt_inv.dot(
            self.correlation_array_eigvecs.conj().T.dot(
                self.cross_correlation_array.dot(
                    self.correlation_array_eigvecs.dot(
                        correlation_array_eigvals_sqrt_inv))))

        # Compute eigendecomposition of low-order linear map
        self.eigvals, self.R_low_order_eigvecs, self.L_low_order_eigvecs =\
            parallel.call_and_bcast(
                util.eig_biorthog, self.low_order_linear_map,
                **{'scale_choice':'left'})


    def compute_decomp(
        self, vec_handles, adv_vec_handles=None, atol=1e-13, rtol=None,
        max_num_eigvals=None):
        """Computes eigendecomposition of low-order linear map approximating
        relationship between vector objects, returning various arrays
        necessary for computing and characterizing DMD modes.

        Args:
            ``vec_handles``: List of handles for vector objects.

        Kwargs:
            ``adv_vec_handles``: List of handles for vector objects advanced in
            time.  If not provided, it is assumed that the vector objects
            describe a sequential time-series. Thus ``vec_handles`` becomes
            ``vec_handles[:-1]`` and ``adv_vec_handles`` becomes
            ``vec_handles[1:]``.

            ``atol``: Level below which DMD eigenvalues are truncated.

            ``rtol``: Maximum relative difference between largest and smallest
            DMD eigenvalues.  Smaller ones are truncated.

            ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
            computed.  This is enforced by truncating the basis onto which the
            approximating linear map is projected.  Computationally, this
            corresponds to truncating the eigendecomposition of the correlation
            array. If set to None, no truncation will be performed, and the
            maximum possible number of DMD eigenvalues will be computed.

        Returns:
            ``eigvals``: 1D array of eigenvalues of low-order linear map, i.e.,
            the DMD eigenvalues.

            ``R_low_order_eigvecs``: Array whose columns are right
            eigenvectors of approximating low-order linear map.

            ``L_low_order_eigvecs``: Array whose columns are left eigenvectors
            of approximating low-order linear map.

            ``correlation_array_eigvals``: 1D array of eigenvalues of
            correlation array.

            ``correlation_array_eigvecs``: Array whose columns are eigenvectors
            of correlation array.
        """
        self.vec_handles = vec_handles
        if adv_vec_handles is not None:
            self.adv_vec_handles = adv_vec_handles
            if len(self.vec_handles) != len(self.adv_vec_handles):
                raise ValueError(('Number of vec_handles and adv_vec_handles'
                    ' is not equal.'))

        # For a sequential dataset, compute correlation array for all vectors.
        # This is more efficient because only one call is made to the inner
        # product routine, even though we don't need the last row/column yet.
        # Later we need all but the last element of the last column, so it is
        # faster to compute all of this now.  Only one extra element is
        # computed, since this is a symmetric inner product array.  Then
        # slice the expanded correlation array accordingly.
        if adv_vec_handles is None:
            self.expanded_correlation_array =\
                self.vec_space.compute_symm_inner_product_array(
                self.vec_handles)
            self.correlation_array = self.expanded_correlation_array[:-1, :-1]
            self.cross_correlation_array = self.expanded_correlation_array[
                :-1, 1:]
        # For non-sequential data, compute the correlation array from the
        # unadvanced snapshots only.  Compute the cross correlation array
        # involving the unadvanced and advanced snapshots separately.
        else:
            self.correlation_array =\
                self.vec_space.compute_symm_inner_product_array(
                    self.vec_handles)
            self.cross_correlation_array =\
                self.vec_space.compute_inner_product_array(
                self.vec_handles, self.adv_vec_handles)

        # Compute eigendecomposition of low-order linear map.
        self.compute_eigendecomp(
            atol=atol, rtol=rtol, max_num_eigvals=max_num_eigvals)

        # Return values
        return (
            self.eigvals,
            self.R_low_order_eigvecs,
            self.L_low_order_eigvecs,
            self.correlation_array_eigvals,
            self.correlation_array_eigvecs)


    def _compute_build_coeffs_exact(self):
        """Compute build coefficients for exact DMD modes."""
        return self.correlation_array_eigvecs.dot(
            np.diag(self.correlation_array_eigvals ** -0.5).dot(
                self.R_low_order_eigvecs.dot(
                    np.diag(self.eigvals ** -1.))))


    def _compute_build_coeffs_proj(self):
        """Compute build coefficients for projected DMD modes."""
        return self.correlation_array_eigvecs.dot(
            np.diag(self.correlation_array_eigvals ** -0.5).dot(
                self.R_low_order_eigvecs))


    def _compute_build_coeffs_adjoint(self):
        """Compute build coefficients for adjoint DMD modes."""
        return self.correlation_array_eigvecs.dot(
            np.diag(self.correlation_array_eigvals ** -0.5).dot(
                self.L_low_order_eigvecs))


    def compute_exact_modes(
        self, mode_indices, mode_handles, adv_vec_handles=None):
        """Computes exact DMD modes and calls ``put`` on them using mode
        handles.

        Args:
            ``mode_indices``: List of indices describing which exact modes to
            compute, e.g. ``range(10)`` or ``[3, 0, 5]``.

            ``mode_handles``: List of handles for exact modes to compute.

        Kwargs:
            ``vec_handles``: List of handles for vector objects. Optional if
            when calling :py:meth:`compute_decomp`.
        """
        # If advanced vec handles are passed in, set the internal attribute,
        if adv_vec_handles is not None:
            self.adv_vec_handles = adv_vec_handles

        # Compute build coefficient array
        build_coeffs_exact = self._compute_build_coeffs_exact()

        # If the internal attribute is set, then compute the modes
        if self.adv_vec_handles is not None:
            self.vec_space.lin_combine(
                mode_handles, self.adv_vec_handles, build_coeffs_exact,
                coeff_array_col_indices=mode_indices)
        # If the internal attribute is not set, then check to see if
        # vec_handles is set.  If so, assume a sequential dataset, in which
        # case adv_vec_handles can be taken from a slice of vec_handles.
        elif self.vec_handles is not None:
            if len(self.vec_handles) - build_coeffs_exact.shape[0] == 1:
                self.vec_space.lin_combine(
                    mode_handles, self.vec_handles[1:], build_coeffs_exact,
                    coeff_array_col_indices=mode_indices)
            else:
                raise(
                    ValueError, (
                    'Number of vec_handles is not correct for a sequential '
                    'dataset.'))
        else:
            raise(
                ValueError,
                'Neither vec_handles nor adv_vec_handles is defined.')


    def compute_proj_modes(self, mode_indices, mode_handles, vec_handles=None):
        """Computes projected DMD modes and calls ``put`` on them using mode
        handles.

        Args:
            ``mode_indices``: List of indices describing which projected modes
            to compute, e.g. ``range(10)`` or ``[3, 0, 5]``.

            ``mode_handles``: List of handles for projected modes to compute.

        Kwargs:
            ``vec_handles``: List of handles for vector objects. Optional if
            when calling :py:meth:`compute_decomp`.
        """
        if vec_handles is not None:
            self.vec_handles = vec_handles

        # Compute build coefficient array
        build_coeffs_proj = self._compute_build_coeffs_proj()

        # For sequential data, the user will provide a list vec_handles that
        # whose length is one larger than the number of rows of the
        # build_coeffs array.  This is to be expected, as vec_handles is
        # essentially partitioned into two sets of handles, each of length one
        # less than vec_handles.
        if len(self.vec_handles) - build_coeffs_proj.shape[0] == 1:
            self.vec_space.lin_combine(
                mode_handles, self.vec_handles[:-1], build_coeffs_proj,
                coeff_array_col_indices=mode_indices)
        # For a non-sequential dataset, the user will provide a list
        # vec_handles whose length is equal to the number of rows in the
        # build_coeffs array.
        elif len(self.vec_handles) == build_coeffs_proj.shape[0]:
            self.vec_space.lin_combine(
                mode_handles, self.vec_handles, build_coeffs_proj,
                coeff_array_col_indices=mode_indices)
        # Otherwise, raise an error, as the number of handles should fit one of
        # the two cases described above.
        else:
            raise ValueError((
                'Number of vec_handles does not match number of columns in '
                'build_coeffs_proj array.'))


    def compute_adjoint_modes(
        self, mode_indices, mode_handles, vec_handles=None):
        """Computes adjoint DMD modes and calls ``put`` on them using mode
        handles.

        Args:
            ``mode_indices``: List of indices describing which projected modes
            to compute, e.g. ``range(10)`` or ``[3, 0, 5]``.

            ``mode_handles``: List of handles for projected modes to compute.

        Kwargs:
            ``vec_handles``: List of handles for vector objects. Optional if
            when calling :py:meth:`compute_decomp`.
        """
        if vec_handles is not None:
            self.vec_handles = vec_handles

        # Compute build coefficient array
        build_coeffs_adjoint = self._compute_build_coeffs_adjoint()

        # For sequential data, the user will provide a list vec_handles that
        # whose length is one larger than the number of rows of the
        # build_coeffs array.  This is to be expected, as vec_handles is
        # essentially partitioned into two sets of handles, each of length one
        # less than vec_handles.
        if len(self.vec_handles) - build_coeffs_adjoint.shape[0] == 1:
            self.vec_space.lin_combine(
                mode_handles, self.vec_handles[:-1], build_coeffs_adjoint,
                coeff_array_col_indices=mode_indices)
        # For a non-sequential dataset, the user will provide a list
        # vec_handles whose length is equal to the number of rows in the
        # build_coeffs array.
        elif len(self.vec_handles) == build_coeffs_adjoint.shape[0]:
            self.vec_space.lin_combine(
                mode_handles, self.vec_handles, build_coeffs_adjoint,
                coeff_array_col_indices=mode_indices)
        # Otherwise, raise an error, as the number of handles should fit one of
        # the two cases described above.
        else:
            raise ValueError((
                'Number of vec_handles does not match number of columns in '
                'build_coeffs_proj array.'))


    def compute_spectrum(self):
        """Computes DMD spectral coefficients.  These coefficients come from a
        biorthogonal projection of the first vector object onto the exact DMD
        modes, which is analytically equivalent to doing a least-squares
        projection onto the projected DMD modes.

        Returns:
            ``spectral_coeffs``: 1D array of DMD spectral coefficients,
            calculated as the magnitudes of the projection coefficients of first
            data vector.  The projection is onto the span of the DMD modes using
            the (biorthogonal) adjoint DMD modes.  Note that this is the same as
            a least-squares projection onto the span of the DMD modes.

        """
        # TODO: maybe allow for user to choose which column to spectrum from?
        # ie first, last, or mean?
        self.spectral_coeffs = np.abs(
            self.L_low_order_eigvecs.conj().T.dot(
                np.diag(self.correlation_array_eigvals ** 0.5).dot(
                    self.correlation_array_eigvecs[0, :]).conj().T)).squeeze()
        return self.spectral_coeffs


    # Note that a biorthogonal projection onto the exact DMD modes is the same
    # as a least squares projection onto the projected DMD modes, so there is
    # only one method for computing the projection coefficients.
    def compute_proj_coeffs(self):
        """Computes projection of vector objects onto DMD modes.  Note that a
        biorthogonal projection onto exact DMD modes is analytically equivalent
        to a least-squares projection onto projected DMD modes.

        Returns:
            ``proj_coeffs``: Array of projection coefficients for vector
            objects, expressed as a linear combination of DMD modes.  Columns
            correspond to vector objects, rows correspond to DMD modes.

            ``adv_proj_coeffs``: Array of projection coefficients for vector
            objects advanced in time, expressed as a linear combination of DMD
            modes.  Columns correspond to vector objects, rows correspond to
            DMD modes.
        """
        self.proj_coeffs = self.L_low_order_eigvecs.conj().T.dot(
            np.diag(self.correlation_array_eigvals ** 0.5).dot(
                self.correlation_array_eigvecs.conj().T))
        self.adv_proj_coeffs = self.L_low_order_eigvecs.conj().T.dot(
            np.diag(self.correlation_array_eigvals ** -0.5).dot(
                self.correlation_array_eigvecs.conj().T.dot(
                    self.cross_correlation_array)))
        return self.proj_coeffs, self.adv_proj_coeffs


def compute_TLSqrDMD_arrays_snaps_method(
    vecs, adv_vecs=None, mode_indices=None, inner_product_weights=None,
    atol=1e-13, rtol=None, max_num_eigvals=None):
    """Computes Total Least Squares DMD modes using data stored in arrays,
    using method of snapshots.

    Args:
        ``vecs``: Array whose columns are data vectors.

    Kwargs:
        ``adv_vecs``: Array whose columns are data vectors advanced in time.
        If not provided, then it is assumed that the vectors describe a
        sequential time-series. Thus ``vecs`` becomes ``vecs[:, :-1]`` and
        ``adv_vecs`` becomes ``vecs[:, 1:]``.

        ``mode_indices``: List of indices describing which modes to compute.
        Examples are ``range(10)`` or ``[3, 0, 6, 8]``.  If no mode indices are
        specified, then all modes will be computed.

        ``inner_product_weights``: 1D or 2D array of inner product weights.
        Corresponds to :math:`W` in inner product :math:`v_1^* W v_2`.

        ``atol``: Level below which eigenvalues of correlation array are
        truncated.

        ``rtol``: Maximum relative difference between largest and smallest
        eigenvalues of correlation array.  Smaller ones are truncated.

        ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
        computed.  This is enforced by truncating the basis onto which the
        approximating linear map is projected.  Computationally, this
        corresponds to truncating the eigendecomposition of the correlation
        array. If set to None, no truncation will be performed, and the
        maximum possible number of DMD eigenvalues will be computed.

    Returns:
        ``res``: Results of DMD computation, stored in a namedtuple with
        the following attributes:

        * ``eigvals``: 1D array of eigenvalues of approximating low-order linear
          map (DMD eigenvalues).


        * ``spectral_coeffs``: 1D array of DMD spectral coefficients, calculated
          as the magnitudes of the projection coefficients of first (de-noised)
          data vector.  The projection is onto the span of the DMD modes using
          the (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``exact_modes``: Array whose columns are exact DMD modes.

        * ``proj_modes``: Array whose columns are projected DMD modes.

        * ``adjoint_modes``: Array whose columns are adjoint DMD modes.

        * ``proj_coeffs``: Array of projection coefficients for (de-noised) data
          vectors expressed as a linear combination of DMD modes.  Columns
          correspond to vector objects, rows correspond to DMD modes. The
          projection is onto the span of the DMD modes using the (biorthogonal)
          adjoint DMD modes.  Note that this is the same as a least-squares
          projection onto the span of the DMD modes.

        * ``adv_proj_coeffs``: Array of projection coefficients for (de-noised)
          data vectors advanced in time, expressed as a linear combination of
          DMD modes.  Columns correspond to vector objects, rows correspond to
          DMD modes. The projection is onto the span of the DMD modes using the
          (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``R_low_order_eigvecs``: Array of right eigenvectors of approximating
          low-order linear map.

        * ``L_low_order_eigvecs``: Array of left eigenvectors of approximating
          low-order linear map.

        * ``sum_correlation_array_eigvals``: 1D array of eigenvalues of
          sum correlation array.

        * ``sum_correlation_array_eigvecs``: Array whose columns are
          eigenvectors of sum correlation array.

        * ``proj_correlation_array_eigvals``: 1D array of eigenvalues of
          projected correlation array.

        * ``proj_correlation_array_eigvecs``: Array whose columns are
          eigenvectors of projected correlation array.

        * ``correlation_array``: Correlation array; elements are inner products
          of data vectors with each other.

        * ``adv_correlation_array``: Advanced correlation array; elements are
          inner products of advanced data vectors with each other.

        * ``cross_correlation_array``: Cross-correlation array; elements are
          inner products of data vectors with data vectors advanced in
          time. Going down rows, the data vector changes; going across columns
          the advanced data vector changes.

        Attributes can be accessed using calls like ``res.exact_modes``.  To
        see all available attributes, use ``print(res)``.

    This uses the method of snapshots, which is faster than the direct method
    (see :py:func:`compute_TLSqrDMD_arrays_direct_method`) when ``vecs`` has
    more rows than columns, i.e., when there are more elements in a vector than
    there are vectors. However, it "squares" this array and its singular
    values, making it slightly less accurate than the direct method.

    Note that max_num_eigvals must be set to a value smaller than the rank of
    the dataset.  In other words, if the projection basis for
    total-least-squares DMD is not truncated, then the algorithm reduces to
    standard DMD.  For over-constrained datasets (number of columns in data
    array is larger than the number of rows), this occurs naturally.  For
    under-constrained datasets, (number of vector objects is smaller than size
    of vector objects), this must be done explicitly by the user.  At this
    time, there is no standard method for choosing a truncation level.  One
    approach is to look at the roll-off of the correlation array eigenvalues,
    which contains information about the "energy" content of each projection
    basis vectors.

    """
    if parallel.is_distributed():
        raise RuntimeError('Cannot run in parallel.')

    # Force data to be arrays (not matrices)
    vecs = np.array(vecs)
    if adv_vecs is not None:
        adv_vecs = np.array(adv_vecs)

    # Set up vector space (for inner products)
    vec_space = VectorSpaceArrays(weights=inner_product_weights)

    # Sequential dataset
    if adv_vecs is None:
        # Compute correlation array for all vectors.
        # This is more efficient because only one call is made to the inner
        # product routine, even though we don't need the last row and column
        # yet.
        expanded_correlation_array = vec_space.compute_symm_inner_product_array(
            vecs)
        correlation_array = expanded_correlation_array[:-1, :-1]
        cross_correlation_array = expanded_correlation_array[:-1, 1:]
        adv_correlation_array = expanded_correlation_array[1:, 1:]
    # Non-sequential data
    else:
        if vecs.shape != adv_vecs.shape:
            raise ValueError(('vecs and adv_vecs are not the same shape.'))
        # Compute the correlation array from the unadvanced snapshots only.
        correlation_array = vec_space.compute_symm_inner_product_array(vecs)
        cross_correlation_array = vec_space.compute_inner_product_array(
            vecs, adv_vecs)
        adv_correlation_array = vec_space.compute_symm_inner_product_array(
            adv_vecs)

    sum_correlation_array_eigvals, sum_correlation_array_eigvecs =\
        util.eigh(
            correlation_array + adv_correlation_array,
            is_positive_definite=True,
            atol=atol, rtol=rtol)

    # Truncate if necessary
    if max_num_eigvals is not None and (
        max_num_eigvals < sum_correlation_array_eigvals.size):
        sum_correlation_array_eigvals = sum_correlation_array_eigvals[
            :max_num_eigvals]
        sum_correlation_array_eigvecs = sum_correlation_array_eigvecs[
            :, :max_num_eigvals]

    # Compute eigendecomposition of projected correlation array
    proj_correlation_array = sum_correlation_array_eigvecs.dot(
        sum_correlation_array_eigvecs.conj().T.dot(
            correlation_array.dot(
                sum_correlation_array_eigvecs.dot(
                    sum_correlation_array_eigvecs.conj().T))))
    proj_correlation_array_eigvals, proj_correlation_array_eigvecs = util.eigh(
        proj_correlation_array, atol=atol, rtol=None, is_positive_definite=True)

    # Truncate if necessary
    if max_num_eigvals is not None and (
        max_num_eigvals < proj_correlation_array_eigvals.size):
        proj_correlation_array_eigvals = proj_correlation_array_eigvals[
            :max_num_eigvals]
        proj_correlation_array_eigvecs = proj_correlation_array_eigvecs[
            :, :max_num_eigvals]

    # Compute low-order linear map
    proj_correlation_array_eigvals_sqrt_inv = np.diag(
        proj_correlation_array_eigvals ** -0.5)
    low_order_linear_map = proj_correlation_array_eigvals_sqrt_inv.dot(
        proj_correlation_array_eigvecs.conj().T.dot(
            sum_correlation_array_eigvecs.dot(
                sum_correlation_array_eigvecs.conj().T.dot(
                    cross_correlation_array.dot(
                        sum_correlation_array_eigvecs.dot(
                            sum_correlation_array_eigvecs.conj().T.dot(
                                proj_correlation_array_eigvecs.dot(
                                    proj_correlation_array_eigvals_sqrt_inv
                                ))))))))

    # Compute eigendecomposition of low-order linear map.
    eigvals, R_low_order_eigvecs, L_low_order_eigvecs = util.eig_biorthog(
        low_order_linear_map, scale_choice='left')

    # Compute build coefficients
    build_coeffs_proj = sum_correlation_array_eigvecs.dot(
        sum_correlation_array_eigvecs.conj().T.dot(
            proj_correlation_array_eigvecs.dot(
                np.diag(proj_correlation_array_eigvals ** -0.5).dot(
                    R_low_order_eigvecs))))
    build_coeffs_exact = build_coeffs_proj.dot(np.diag(eigvals ** -1.))
    build_coeffs_adjoint = sum_correlation_array_eigvecs.dot(
        sum_correlation_array_eigvecs.conj().T.dot(
            proj_correlation_array_eigvecs.dot(
                np.diag(proj_correlation_array_eigvals ** -0.5).dot(
                    L_low_order_eigvecs))))

    # Compute spectral coefficients
    spectral_coeffs = np.abs(
        L_low_order_eigvecs.conj().T.dot(
            np.diag(proj_correlation_array_eigvals ** 0.5).dot(
                proj_correlation_array_eigvecs[0, :].T))).squeeze()

    # Compute projection coefficients
    proj_coeffs = L_low_order_eigvecs.conj().T.dot(
        np.diag(proj_correlation_array_eigvals ** 0.5).dot(
            proj_correlation_array_eigvecs.conj().T))
    adv_proj_coeffs = L_low_order_eigvecs.conj().T.dot(
        np.diag(proj_correlation_array_eigvals ** -0.5).dot(
            proj_correlation_array_eigvecs.conj().T.dot(
                sum_correlation_array_eigvecs.dot(
                    sum_correlation_array_eigvecs.conj().T.dot(
                        cross_correlation_array.dot(
                            sum_correlation_array_eigvecs.dot(
                                sum_correlation_array_eigvecs.conj().T
                            )))))))

    # For sequential data, user must provide one more vec than columns of
    # build_coeffs.
    if vecs.shape[1] - build_coeffs_exact.shape[0] == 1:
        exact_modes = vec_space.lin_combine(
            vecs[:, 1:], build_coeffs_exact,
            coeff_array_col_indices=mode_indices)
        proj_modes = vec_space.lin_combine(
            vecs[:, :-1], build_coeffs_proj,
            coeff_array_col_indices=mode_indices)
        adjoint_modes = vec_space.lin_combine(
            vecs[:, :-1], build_coeffs_adjoint,
            coeff_array_col_indices=mode_indices)
    # For non-sequential data, user must provide as many vecs as columns of
    # build_coeffs.
    elif vecs.shape[1] == build_coeffs_exact.shape[0]:
        exact_modes = vec_space.lin_combine(
            adv_vecs, build_coeffs_exact, coeff_array_col_indices=mode_indices)
        proj_modes = vec_space.lin_combine(
            vecs, build_coeffs_proj, coeff_array_col_indices=mode_indices)
        adjoint_modes = vec_space.lin_combine(
            vecs, build_coeffs_adjoint, coeff_array_col_indices=mode_indices)
    else:
        raise ValueError(('Number of cols in vecs does not match '
            'number of rows in build_coeffs array.'))

    # Return a namedtuple
    TLSqrDMD_results = namedtuple(
        'TLSqrDMD_results', [
            'eigvals','spectral_coeffs',
            'exact_modes', 'proj_modes', 'adjoint_modes',
            'proj_coeffs', 'adv_proj_coeffs',
            'R_low_order_eigvecs', 'L_low_order_eigvecs',
            'sum_correlation_array_eigvals', 'sum_correlation_array_eigvecs',
            'proj_correlation_array_eigvals', 'proj_correlation_array_eigvecs',
            'correlation_array', 'adv_correlation_array',
            'cross_correlation_array'])
    return TLSqrDMD_results(
        eigvals=eigvals, spectral_coeffs=spectral_coeffs,
        exact_modes=exact_modes,
        proj_modes=proj_modes, adjoint_modes=adjoint_modes,
        proj_coeffs=proj_coeffs, adv_proj_coeffs=adv_proj_coeffs,
        R_low_order_eigvecs=R_low_order_eigvecs,
        L_low_order_eigvecs=L_low_order_eigvecs,
        sum_correlation_array_eigvals=sum_correlation_array_eigvals,
        sum_correlation_array_eigvecs=sum_correlation_array_eigvecs,
        proj_correlation_array_eigvals=proj_correlation_array_eigvals,
        proj_correlation_array_eigvecs=proj_correlation_array_eigvecs,
        correlation_array=correlation_array,
        adv_correlation_array=adv_correlation_array,
        cross_correlation_array=cross_correlation_array)


def compute_TLSqrDMD_arrays_direct_method(
    vecs, adv_vecs=None, mode_indices=None, inner_product_weights=None,
    atol=1e-13, rtol=None, max_num_eigvals=None):
    """Computes Total Least Squares DMD modes using data stored in arrays,
    using direct method.

    Args:
        ``vecs``: Array whose columns are data vectors.

    Kwargs:
        ``adv_vecs``: Array whose columns are data vectors advanced in time.
        If not provided, then it is assumed that the vectors describe a
        sequential time-series. Thus ``vecs`` becomes ``vecs[:, :-1]`` and
        ``adv_vecs`` becomes ``vecs[:, 1:]``.

        ``mode_indices``: List of indices describing which modes to compute.
        Examples are ``range(10)`` or ``[3, 0, 6, 8]``.  If no mode indices are
        specified, then all modes will be computed.

        ``inner_product_weights``: 1D or 2D array of inner product weights.
        Corresponds to :math:`W` in inner product :math:`v_1^* W v_2`.

        ``atol``: Level below which eigenvalues of correlation array are
        truncated.

        ``rtol``: Maximum relative difference between largest and smallest
        eigenvalues of correlation array.  Smaller ones are truncated.

        ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
        computed.  This is enforced by truncating the basis onto which the
        approximating linear map is projected.  Computationally, this
        corresponds to truncating the eigendecomposition of the correlation
        array. If set to None, no truncation will be performed, and the
        maximum possible number of DMD eigenvalues will be computed.

    Returns:
        ``res``: Results of DMD computation, stored in a namedtuple with
        the following attributes:

        * ``eigvals``: 1D array of eigenvalues of approximating low-order linear
          map (DMD eigenvalues).

        * ``spectral_coeffs``: 1D array of DMD spectral coefficients, calculated
          as the magnitudes of the projection coefficients of first (de-noised)
          data vector.  The projection is onto the span of the DMD modes using
          the (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``exact_modes``: Array whose columns are exact DMD modes.

        * ``proj_modes``: Array whose columns are projected DMD modes.

        * ``adjoint_modes``: Array whose columns are adjoint DMD modes.

        * ``proj_coeffs``: Array of projection coefficients for (de-noised) data
          vectors expressed as a linear combination of DMD modes.  Columns
          correspond to vector objects, rows correspond to DMD modes. The
          projection is onto the span of the DMD modes using the (biorthogonal)
          adjoint DMD modes.  Note that this is the same as a least-squares
          projection onto the span of the DMD modes.

        * ``adv_proj_coeffs``: Array of projection coefficients for (de-noised)
          data vectors advanced in time, expressed as a linear combination of
          DMD modes.  Columns correspond to vector objects, rows correspond to
          DMD modes. The projection is onto the span of the DMD modes using the
          (biorthogonal) adjoint DMD modes.  Note that this is the same as a
          least-squares projection onto the span of the DMD modes.

        * ``R_low_order_eigvecs``: Array of right eigenvectors of approximating
          low-order linear map.

        * ``L_low_order_eigvecs``: Array of left eigenvectors of approximating
          low-order linear map.

        * ``sum_correlation_array_eigvals``: 1D array of eigenvalues of
          sum correlation array.

        * ``sum_correlation_array_eigvecs``: Array whose columns are
          eigenvectors of sum correlation array.

        * ``proj_correlation_array_eigvals``: 1D array of eigenvalues of
          projected correlation array.

        * ``proj_correlation_array_eigvecs``: Array whose columns are
          eigenvectors of projected correlation array.

        Attributes can be accessed using calls like ``res.exact_modes``.  To
        see all available attributes, use ``print(res)``.

    This method does not square the array of vectors as in the method of
    snapshots (:py:func:`compute_DMD_arrays_snaps_method`). It's slightly
    more accurate, but slower when the number of elements in a vector is more
    than the number of vectors, i.e.,  when ``vecs`` has more rows than
    columns.

    Note that max_num_eigvals must be set to a value smaller than the rank of
    the dataset.  In other words, if the projection basis for
    total-least-squares DMD is not truncated, then the algorithm reduces to
    standard DMD.  For over-constrained datasets (number of columns in data
    array is larger than the number of rows), this occurs naturally.  For
    under-constrained datasets, (number of vector objects is smaller than size
    of vector objects), this must be done explicitly by the user.  At this
    time, there is no standard method for choosing a truncation level.  One
    approach is to look at the roll-off of the correlation array eigenvalues,
    which contains information about the "energy" content of each projection
    basis vectors.

    """
    if parallel.is_distributed():
        raise RuntimeError('Cannot run in parallel.')

    # Force data to be arrays (not matrices)
    vecs = np.array(vecs)
    if adv_vecs is not None:
        adv_vecs = np.array(adv_vecs)

    # Set up vector space (for inner products)
    vec_space = VectorSpaceArrays(weights=inner_product_weights)

    # Weight data as necessary
    if inner_product_weights is None:
        vecs_weighted = vecs
        if adv_vecs is not None:
            adv_vecs_weighted = adv_vecs
    elif inner_product_weights.ndim == 1:
        sqrt_weights = np.diag(inner_product_weights ** 0.5)
        vecs_weighted = sqrt_weights.dot(vecs)
        if adv_vecs is not None:
            adv_vecs_weighted = sqrt_weights.dot(adv_vecs)
    elif inner_product_weights.ndim == 2:
        if inner_product_weights.shape[0] > 500:
            print('Warning: Cholesky decomposition could be time consuming.')
        sqrt_weights = np.linalg.cholesky(inner_product_weights).conj().T
        vecs_weighted = sqrt_weights.dot(vecs)
        if adv_vecs is not None:
            adv_vecs_weighted = sqrt_weights.dot(adv_vecs)

    # Compute projections of original data (to de-noise).  First consider the
    # sequential data case.
    if adv_vecs is None:
        stacked_U, stacked_sing_vals, sum_correlation_array_eigvecs =util.svd(
            np.vstack((vecs_weighted[:, :-1], vecs_weighted[:, 1:])),
            atol=atol, rtol=rtol)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < stacked_sing_vals.size):
            stacked_U = stacked_U[:, :max_num_eigvals]
            stacked_sing_vals = stacked_sing_vals[:max_num_eigvals]
            sum_correlation_array_eigvecs = sum_correlation_array_eigvecs[
                :, :max_num_eigvals]

        # Project original data to de-noise
        vecs_proj = vecs[:, :-1].dot(
            sum_correlation_array_eigvecs.dot(
                sum_correlation_array_eigvecs.conj().T))
        adv_vecs_proj = vecs[:, 1:].dot(
            sum_correlation_array_eigvecs.dot(
                sum_correlation_array_eigvecs.conj().T))
    # Non-sequential data case
    else:
        if vecs.shape != adv_vecs.shape:
            raise ValueError(('vecs and adv_vecs are not the same shape.'))
        stacked_U, stacked_sing_vals, sum_correlation_array_eigvecs = util.svd(
            np.vstack((vecs_weighted, adv_vecs_weighted)),
            atol=atol, rtol=rtol)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < stacked_sing_vals.size):
            stacked_U = stacked_U[:, :max_num_eigvals]
            stacked_sing_vals = stacked_sing_vals[:max_num_eigvals]
            sum_correlation_array_eigvecs = sum_correlation_array_eigvecs[
                :, :max_num_eigvals]

        # Project original data to de-noise
        vecs_proj = vecs.dot(
            sum_correlation_array_eigvecs.dot(
                sum_correlation_array_eigvecs.conj().T))
        adv_vecs_proj = adv_vecs.dot(
            sum_correlation_array_eigvecs.dot(
                sum_correlation_array_eigvecs.conj().T))

    # Now proceed with DMD of projected data
    sum_correlation_array_eigvals = stacked_sing_vals ** 2
    DMD_res = compute_DMD_arrays_direct_method(
        vecs_proj, adv_vecs=adv_vecs_proj, mode_indices=mode_indices,
        inner_product_weights=inner_product_weights, atol=atol, rtol=rtol,
        max_num_eigvals=max_num_eigvals)

    # Return a namedtuple
    TLSqrDMD_results = namedtuple(
        'TLSqrDMD_results', [
            'exact_modes', 'proj_modes', 'adjoint_modes',
            'spectral_coeffs', 'proj_coeffs', 'adv_proj_coeffs',
            'eigvals', 'R_low_order_eigvecs', 'L_low_order_eigvecs',
            'sum_correlation_array_eigvals', 'sum_correlation_array_eigvecs',
            'proj_correlation_array_eigvals', 'proj_correlation_array_eigvecs'])
    return TLSqrDMD_results(
        exact_modes=DMD_res.exact_modes,
        proj_modes=DMD_res.proj_modes, adjoint_modes=DMD_res.adjoint_modes,
        spectral_coeffs=DMD_res.spectral_coeffs,
        proj_coeffs=DMD_res.proj_coeffs,
        adv_proj_coeffs=DMD_res.adv_proj_coeffs,
        eigvals=DMD_res.eigvals,
        R_low_order_eigvecs=DMD_res.R_low_order_eigvecs,
        L_low_order_eigvecs=DMD_res.L_low_order_eigvecs,
        sum_correlation_array_eigvals=sum_correlation_array_eigvals,
        sum_correlation_array_eigvecs=sum_correlation_array_eigvecs,
        proj_correlation_array_eigvals=DMD_res.correlation_array_eigvals,
        proj_correlation_array_eigvecs=DMD_res.correlation_array_eigvecs)


class TLSqrDMDHandles(DMDHandles):
    """Total Least Squares Dynamic Mode Decomposition implemented for large
    datasets.

    Kwargs:
        ``inner_product``: Function that computes inner product of two vector
        objects.

        ``put_array``: Function to put an array out of modred, e.g., write it to
        file.

        ``get_array``: Function to get an array into modred, e.g., load it from
        file.

        ``max_vecs_per_node``: Maximum number of vectors that can be stored in
        memory, per node.

        ``verbosity``: 1 prints progress and warnings, 0 prints almost nothing.

    Computes Total-Least-Squares DMD modes from vector objects (or handles).
    It uses :py:class:`vectorspace.VectorSpaceHandles` for low level functions.

    Usage::

      myDMD = TLSqrDMDHandles(inner_product=my_inner_product)
      myDMD.compute_decomp(vec_handles, max_num_eigvals=100)
      myDMD.compute_exact_modes(range(50), mode_handles)

    The total-least-squares variant of DMD accounts for an asymmetric treatment
    of the non-time-advanced and time-advanced data in standard DMD, in which
    the former is assumed to be perfect information, whereas all noise is
    limited to the latter.  In many cases, since the non-time-advanced and
    advanced data come from the same source, noise could be present in either.

    Note that max_num_eigvals must be set to a value smaller than the rank of
    the dataset.  In other words, if the projection basis for
    total-least-squares DMD is not truncated, then the algorithm reduces to
    standard DMD.  For over-constrained datasets (number of vector objects  is
    larger than the size of each vector object), this occurs naturally.  For
    under-constrained datasets, (number of vector objects is smaller than size
    of vector objects), this must be done explicitly by the user.  At this
    time, there is no standard method for choosing a truncation level.  One
    approach is to look at the roll-off of the correlation array eigenvalues,
    which contains information about the "energy" content of each projection
    basis vectors.

    Also, note that :class:`TLSqrDMDHandles` inherits from
    :class:`DMDHandles`, so
    certain methods are available, even though they are not
    implemented/documented here (namely several ``put`` functions).

    See also :func:`compute_TLSqrDMD_arrays_snaps_method`,
    :func:`compute_TLSqrDMD_arrays_direct_method`, and :mod:`vectors`.
    """
    def __init__(
        self, inner_product=None, get_array=util.load_array_text,
        put_array=util.save_array_text, max_vecs_per_node=None, verbosity=1):
        """Constructor"""
        self.get_array = get_array
        self.put_array = put_array
        self.verbosity = verbosity
        self.eigvals = None
        self.correlation_array = None
        self.cross_correlation_array = None
        self.adv_correlation_array = None
        self.sum_correlation_array = None
        self.sum_correlation_array_eigvals = None
        self.sum_correlation_array_eigvecs = None
        self.proj_correlation_array = None
        self.proj_correlation_array_eigvals = None
        self.proj_correlation_array_eigvecs = None
        self.low_order_linear_map = None
        self.L_low_order_eigvecs = None
        self.R_low_order_eigvecs = None
        self.spectral_coeffs = None
        self.proj_coeffs = None
        self.adv_proj_coeffs = None
        self.vec_space = VectorSpaceHandles(inner_product=inner_product,
            max_vecs_per_node=max_vecs_per_node, verbosity=verbosity)
        self.vec_handles = None
        self.adv_vec_handles = None


    def compute_eigendecomp(self, atol=1e-13, rtol=None, max_num_eigvals=None):
        """Computes eigendecompositions of correlation array and approximating
        low-order linear map.

        Kwargs:
            ``atol``: Level below which eigenvalues of correlation array are
            truncated.

            ``rtol``: Maximum relative difference between largest and smallest
            eigenvalues of correlation array.  Smaller ones are truncated.

            ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
            computed.  This is enforced by truncating the basis onto which the
            approximating linear map is projected.  Computationally, this
            corresponds to truncating the eigendecomposition of the correlation
            array. If set to None, no truncation will be performed, and the
            maximum possible number of DMD eigenvalues will be computed.

        Useful if you already have the correlation array and cross-correlation
        array and want to avoid recomputing them.

        Usage::

          TLSqrDMD.correlation_array = pre_existing_correlation_array
          TLSqrDMD.adv_correlation_array = pre_existing_adv_correlation_array
          TLSqrDMD.cross_correlation_array =
              pre_existing_cross_correlation_array
          TLSqrDMD.compute_eigendecomp()
          TLSqrDMD.compute_exact_modes(
              mode_idx_list, mode_handles, adv_vec_handles=adv_vec_handles)

        Another way to use this is to compute TLSqrDMD using different basis
        truncation levels for the projection of the approximating linear map.
        Start by either computing a full decomposition or by loading
        pre-computed correlation and cross-correlation arrays.

        Usage::

          # Start with a full decomposition
          DMD_eigvals, correlation_array_eigvals = TLSqrDMD.compute_decomp(
              vec_handles)[0, 3]

          # Do some processing to determine the truncation level, maybe based
          # on the DMD eigenvalues and correlation array eigenvalues
          desired_num_eigvals = my_post_processing_func(
              DMD_eigvals, correlation_array_eigvals)

          # Do a truncated decomposition
          DMD_eigvals_trunc = TLSqrDMD.compute_eigendecomp(
            max_num_eigvals=desired_num_eigvals)

          # Compute modes for truncated decomposition
          TLSqrDMD.compute_exact_modes(
              mode_idx_list, mode_handles, adv_vec_handles=adv_vec_handles)

        Since it doesn't overwrite the correlation and cross-correlation
        arrays, ``compute_eigendecomp`` can be called many times in a row to
        do computations for different truncation levels.  However, the results
        of the decomposition (e.g., ``self.eigvals``) do get overwritten, so
        you may want to call a ``put`` method to save those results somehow.

        Note that the truncation level (corresponding to ``max_num_eigvals``)
        must be set to a value smaller than the rank of the dataset, otherwise
        total-least-squares DMD reduces to standard DMD.  This occurs naturally
        for over-constrained datasets, but must be enforced by the user for
        under-constrined datasets.
        """
        # Compute eigendecomposition of stacked correlation array
        self.sum_correlation_array = (
            self.correlation_array + self.adv_correlation_array)
        (self.sum_correlation_array_eigvals,
        self.sum_correlation_array_eigvecs) = parallel.call_and_bcast(
            util.eigh, self.sum_correlation_array,
            atol=atol, rtol=None, is_positive_definite=True)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < self.sum_correlation_array_eigvals.size):
            self.sum_correlation_array_eigvals =\
                self.sum_correlation_array_eigvals[:max_num_eigvals]
            self.sum_correlation_array_eigvecs =\
                self.sum_correlation_array_eigvecs[:, :max_num_eigvals]

        # Compute eigendecomposition of projected correlation array
        self.proj_correlation_array = self.sum_correlation_array_eigvecs.dot(
            self.sum_correlation_array_eigvecs.conj().T.dot(
                self.correlation_array.dot(
                    self.sum_correlation_array_eigvecs.dot(
                        self.sum_correlation_array_eigvecs.conj().T))))
        (self.proj_correlation_array_eigvals,
        self.proj_correlation_array_eigvecs) = parallel.call_and_bcast(
            util.eigh, self.proj_correlation_array,
            atol=atol, rtol=None, is_positive_definite=True)

        # Truncate if necessary
        if max_num_eigvals is not None and (
            max_num_eigvals < self.proj_correlation_array_eigvals.size):
            self.proj_correlation_array_eigvals =\
                self.proj_correlation_array_eigvals[:max_num_eigvals]
            self.proj_correlation_array_eigvecs =\
                self.proj_correlation_array_eigvecs[:, :max_num_eigvals]

        # Compute low-order linear map
        proj_correlation_array_eigvals_sqrt_inv = np.diag(
            self.proj_correlation_array_eigvals ** -0.5)
        self.low_order_linear_map = proj_correlation_array_eigvals_sqrt_inv.dot(
            self.proj_correlation_array_eigvecs.conj().T.dot(
                self.sum_correlation_array_eigvecs.dot(
                    self.sum_correlation_array_eigvecs.conj().T.dot(
                        self.cross_correlation_array.dot(
                            self.sum_correlation_array_eigvecs.dot(
                                self.sum_correlation_array_eigvecs.conj().T.dot(
                                    self.proj_correlation_array_eigvecs.dot(
                                        proj_correlation_array_eigvals_sqrt_inv
                                    ))))))))

        # Compute eigendecomposition of low-order linear map
        self.eigvals, self.R_low_order_eigvecs, self.L_low_order_eigvecs =\
            parallel.call_and_bcast(
            util.eig_biorthog, self.low_order_linear_map,
            **{'scale_choice':'left'})


    def compute_decomp(
        self, vec_handles, adv_vec_handles=None, atol=1e-13, rtol=None,
        max_num_eigvals=None):
        """Computes eigendecomposition of low-order linear map approximating
        relationship between vector objects, returning various arrays
        necessary for computing and characterizing DMD modes.

        Args:
            ``vec_handles``: List of handles for vector objects.

        Kwargs:
            ``adv_vec_handles``: List of handles for vector objects advanced in
            time.  If not provided, it is assumed that the vector objects
            describe a sequential time-series. Thus ``vec_handles`` becomes
            ``vec_handles[:-1]`` and ``adv_vec_handles`` becomes
            ``vec_handles[1:]``.

            ``atol``: Level below which DMD eigenvalues are truncated.

            ``rtol``: Maximum relative difference between largest and smallest
            DMD eigenvalues.  Smaller ones are truncated.

            ``max_num_eigvals``: Maximum number of DMD eigenvalues that will be
            computed.  This is enforced by truncating the basis onto which the
            approximating linear map is projected.  Computationally, this
            corresponds to truncating the eigendecomposition of the correlation
            array. If set to None, no truncation will be performed, and the
            maximum possible number of DMD eigenvalues will be computed.

        Returns:
            ``eigvals``: 1D array of eigenvalues of low-order linear map, i.e.,
            the DMD eigenvalues.

            ``R_low_order_eigvecs``: Array whose columns are right
            eigenvectors of approximating low-order linear map.

            ``L_low_order_eigvecs``: Array whose columns are left eigenvectors
            of approximating low-order linear map.

            ``sum_correlation_array_eigvals``: 1D array of eigenvalues of
            sum correlation array.

            ``sum_correlation_array_eigvecs``: Array whose columns are
            eigenvectors of sum correlation array.

            ``proj_correlation_array_eigvals``: 1D array of eigenvalues of
            projected correlation array.

            ``proj_correlation_array_eigvecs``: Array whose columns are
            eigenvectors of projected correlation array.

        Note that the truncation level (corresponding to ``max_num_eigvals``)
        must be set to a value smaller than the rank of the dataset, otherwise
        total-least-squares DMD reduces to standard DMD.  This occurs naturally
        for over-constrained datasets, but must be enforced by the user for
        under-constrined datasets.
        """
        self.vec_handles = vec_handles
        if adv_vec_handles is not None:
            self.adv_vec_handles = adv_vec_handles
            if len(self.vec_handles) != len(self.adv_vec_handles):
                raise ValueError(('Number of vec_handles and adv_vec_handles'
                    ' is not equal.'))

        # For a sequential dataset, compute correlation array for all vectors.
        # This is more efficient because only one call is made to the inner
        # product routine, even though we don't need the last row/column yet.
        # Later we need all but the last element of the last column, so it is
        # faster to compute all of this now.  Only one extra element is
        # computed, since this is a symmetric inner product array.  Then
        # slice the expanded correlation array accordingly.
        if adv_vec_handles is None:
            self.expanded_correlation_array =\
                self.vec_space.compute_symm_inner_product_array(
                    self.vec_handles)
            self.correlation_array = self.expanded_correlation_array[:-1, :-1]
            self.cross_correlation_array = self.expanded_correlation_array[
                :-1, 1:]
            self.adv_correlation_array = self.expanded_correlation_array[1:, 1:]
        # For non-sequential data, compute the correlation array from the
        # unadvanced snapshots only.  Compute the cross correlation array
        # involving the unadvanced and advanced snapshots separately.
        else:
            self.correlation_array =\
                self.vec_space.compute_symm_inner_product_array(
                    self.vec_handles)
            self.cross_correlation_array =\
                self.vec_space.compute_inner_product_array(
                self.vec_handles, self.adv_vec_handles)
            self.adv_correlation_array =\
                self.vec_space.compute_symm_inner_product_array(
                    self.adv_vec_handles)

        # Compute eigendecomposition of low-order linear map.
        self.compute_eigendecomp(
            atol=atol, rtol=rtol, max_num_eigvals=max_num_eigvals)

        # Return values
        return (
            self.eigvals,
            self.R_low_order_eigvecs,
            self.L_low_order_eigvecs,
            self.sum_correlation_array_eigvals,
            self.sum_correlation_array_eigvecs,
            self.proj_correlation_array_eigvals,
            self.proj_correlation_array_eigvecs)


    def _compute_build_coeffs_exact(self):
        """Compute build coefficients for exact DMD modes."""
        return self.sum_correlation_array_eigvecs.dot(
            self.sum_correlation_array_eigvecs.conj().T.dot(
                self.proj_correlation_array_eigvecs.dot(
                    np.diag(self.proj_correlation_array_eigvals ** -0.5)).dot(
                        self.R_low_order_eigvecs.dot(
                            np.diag(self.eigvals ** -1.)))))


    def _compute_build_coeffs_proj(self):
        """Compute build coefficients for projected DMD modes."""
        return self.sum_correlation_array_eigvecs.dot(
            self.sum_correlation_array_eigvecs.conj().T.dot(
                self.proj_correlation_array_eigvecs.dot(
                    np.diag(self.proj_correlation_array_eigvals ** -0.5).dot(
                        self.R_low_order_eigvecs))))


    def _compute_build_coeffs_adjoint(self):
        """Compute build coefficients for adjoint DMD modes."""
        return self.sum_correlation_array_eigvecs.dot(
            self.sum_correlation_array_eigvecs.conj().T.dot(
                self.proj_correlation_array_eigvecs.dot(
                    np.diag(self.proj_correlation_array_eigvals ** -0.5).dot(
                        self.L_low_order_eigvecs))))


    def get_decomp(
        self, eigvals_src, R_low_order_eigvecs_src, L_low_order_eigvecs_src,
        sum_correlation_array_eigvals_src,
        sum_correlation_array_eigvecs_src, proj_correlation_array_eigvals_src,
        proj_correlation_array_eigvecs_src):
        """Gets the decomposition arrays from sources (memory or file).

        Args:
            ``eigvals_src``: Source from which to retrieve eigenvalues of
            approximating low-order linear map (DMD eigenvalues).

            ``R_low_order_eigvecs_src``: Source from which to retrieve right
            eigenvectors of approximating low-order linear DMD map.

            ``L_low_order_eigvecs_src``: Source from which to retrieve left
            eigenvectors of approximating low-order linear DMD map.

            ``sum_correlation_array_eigvals_src``: Source from which to
            retrieve eigenvalues of sum correlation array.

            ``sum_correlation_array_eigvecs_src``: Source from which to
            retrieve eigenvectors of sum correlation array.

            ``proj_correlation_array_eigvals_src``: Source from which to
            retrieve eigenvalues of projected correlation array.

            ``proj_correlation_array_eigvecs_src``: Source from which to
            retrieve eigenvectors of projected correlation array.
        """
        self.eigvals = np.squeeze(np.array(
            parallel.call_and_bcast(self.get_array, eigvals_src)))
        self.R_low_order_eigvecs = parallel.call_and_bcast(
            self.get_array, R_low_order_eigvecs_src)
        self.L_low_order_eigvecs = parallel.call_and_bcast(
            self.get_array, L_low_order_eigvecs_src)
        self.sum_correlation_array_eigvals = np.squeeze(np.array(
            parallel.call_and_bcast(
            self.get_array, sum_correlation_array_eigvals_src)))
        self.sum_correlation_array_eigvecs = parallel.call_and_bcast(
            self.get_array, sum_correlation_array_eigvecs_src)
        self.proj_correlation_array_eigvals = np.squeeze(np.array(
            parallel.call_and_bcast(
            self.get_array, proj_correlation_array_eigvals_src)))
        self.proj_correlation_array_eigvecs = parallel.call_and_bcast(
            self.get_array, proj_correlation_array_eigvecs_src)


    def get_adv_correlation_array(self, src):
        """Gets the advanced correlation array from source (memory or file).

        Args:
            ``src``: Source from which to retrieve advanced correlation array.
        """
        self.adv_correlation_array = parallel.call_and_bcast(
            self.get_array, src)


    def get_sum_correlation_array(self, src):
        """Gets the sum correlation array from source (memory or file).

        Args:
            ``src``: Source from which to retrieve sum correlation array.
        """
        self.sum_correlation_array = parallel.call_and_bcast(
            self.get_array, src)


    def get_proj_correlation_array(self, src):
        """Gets the projected correlation array from source (memory or file).

        Args:
            ``src``: Source from which to retrieve projected correlation
            array.
        """
        self.proj_correlation_array = parallel.call_and_bcast(
            self.get_array, src)


    def put_decomp(
        self, eigvals_dest, R_low_order_eigvecs_dest, L_low_order_eigvecs_dest,
        sum_correlation_array_eigvals_dest,
        sum_correlation_array_eigvecs_dest,
        proj_correlation_array_eigvals_dest,
        proj_correlation_array_eigvecs_dest):
        """Puts the decomposition arrays in destinations (file or memory).

        Args:
            ``eigvals_dest``: Destination in which to put eigenvalues of
            approximating low-order linear map (DMD eigenvalues).

            ``R_low_order_eigvecs_dest``: Destination in which to put right
            eigenvectors of approximating low-order linear map.

            ``L_low_order_eigvecs_dest``: Destination in which to put left
            eigenvectors of approximating low-order linear map.

            ``sum_correlation_array_eigvals_dest``: Destination in which to
            put eigenvalues of sum correlation array.

            ``sum_correlation_array_eigvecs_dest``: Destination in which to
            put eigenvectors of sum correlation array.

            ``proj_correlation_array_eigvals_dest``: Destination in which to put
            eigenvalues of projected correlation array.

            ``proj_correlation_array_eigvecs_dest``: Destination in which to put
            eigenvectors of projected correlation array.
        """
        # Don't check if rank is zero because the following methods do.
        self.put_eigvals(eigvals_dest)
        self.put_R_low_order_eigvecs(R_low_order_eigvecs_dest)
        self.put_L_low_order_eigvecs(L_low_order_eigvecs_dest)
        self.put_sum_correlation_array_eigvals(
            sum_correlation_array_eigvals_dest)
        self.put_sum_correlation_array_eigvecs(
            sum_correlation_array_eigvecs_dest)
        self.put_proj_correlation_array_eigvals(
            proj_correlation_array_eigvals_dest)
        self.put_proj_correlation_array_eigvecs(
            proj_correlation_array_eigvecs_dest)


    def put_correlation_array_eigvals(self, dest):
        """This method is not available for total least squares DMD"""
        raise NotImplementedError(
            'This method is not available.  Use '
            'put_sum_correlation_array_eigvals instead.')


    def put_correlation_array_eigvecs(self, dest):
        """This method is not available for total least squares DMD"""
        raise NotImplementedError(
            'This method is not available.  Use '
            'put_sum_correlation_array_eigvecs instead.')


    def put_sum_correlation_array_eigvals(self, dest):
        """Puts eigenvalues of sum correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.sum_correlation_array_eigvals, dest)
        parallel.barrier()


    def put_sum_correlation_array_eigvecs(self, dest):
        """Puts eigenvectors of sum correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.sum_correlation_array_eigvecs, dest)
        parallel.barrier()


    def put_proj_correlation_array_eigvals(self, dest):
        """Puts eigenvalues of projected correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.proj_correlation_array_eigvals, dest)
        parallel.barrier()


    def put_proj_correlation_array_eigvecs(self, dest):
        """Puts eigenvectors of projected correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.proj_correlation_array_eigvecs, dest)
        parallel.barrier()


    def put_adv_correlation_array(self, dest):
        """Puts advanced correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.adv_correlation_array, dest)
        parallel.barrier()


    def put_sum_correlation_array(self, dest):
        """Puts sum correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.sum_correlation_array, dest)
        parallel.barrier()


    def put_proj_correlation_array(self, dest):
        """Puts projected correlation array to ``dest``."""
        if parallel.is_rank_zero():
            self.put_array(self.proj_correlation_array, dest)
        parallel.barrier()


    def compute_spectrum(self):
        """Computes DMD spectral coefficients.  These coefficients come from a
        biorthogonal projection of the first (de-noised) vector object onto the
        exact DMD modes, which is analytically equivalent to doing a
        least-squares projection onto the projected DMD modes.

        Returns:
            ``spectral_coeffs``: 1D array of DMD spectral coefficients.
        """
        # TODO: maybe allow for user to choose which column to spectrum from?
        # ie first, last, or mean?
        self.spectral_coeffs = np.abs(self.L_low_order_eigvecs.conj().T.dot(
            np.diag(self.proj_correlation_array_eigvals ** 0.5).dot(
                self.proj_correlation_array_eigvecs[0, :]).conj().T)).squeeze()
        return self.spectral_coeffs


    # Note that a biorthogonal projection onto the exact DMD modes is the same
    # as a least squares projection onto the projected DMD modes, so there is
    # only one method for computing the projection coefficients.
    def compute_proj_coeffs(self):
        """Computes projection of (de-noised) vector objects onto DMD modes.
        Note that a biorthogonal projection onto exact DMD modes is
        analytically equivalent to a least-squares projection onto projected
        DMD modes.

        Returns:
            ``proj_coeffs``: Array of projection coefficients for
            (de-noised)vector objects, expressed as a linear combination of DMD
            modes.  Columns correspond to vector objects, rows correspond to
            DMD modes.

            ``adv_proj_coeffs``: Array of projection coefficients for
            (de-noised) vector objects advanced in time, expressed as a linear
            combination of DMD modes.  Columns correspond to vector objects,
            rows correspond to DMD modes.
        """
        self.proj_coeffs = self.L_low_order_eigvecs.conj().T.dot(
            np.diag(self.proj_correlation_array_eigvals ** 0.5).dot(
                self.proj_correlation_array_eigvecs.conj().T))
        self.adv_proj_coeffs = self.L_low_order_eigvecs.conj().T.dot(
            np.diag(self.proj_correlation_array_eigvals ** -0.5).dot(
                self.proj_correlation_array_eigvecs.conj().T.dot(
                    self.sum_correlation_array_eigvecs.dot(
                        self.sum_correlation_array_eigvecs.conj().T.dot(
                            self.cross_correlation_array.dot(
                                self.sum_correlation_array_eigvecs.dot(
                                    self.sum_correlation_array_eigvecs.conj().T
                                )))))))
        return self.proj_coeffs, self.adv_proj_coeffs
