Source code for Muscat.LinAlg.SVD

# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
from typing import Optional

import numpy as np

from Muscat.Types import MuscatFloat, MuscatIndex
from Muscat.Helpers.Logger import Warning




def TruncatedSVDSymLower(matrix, epsilon: Optional[MuscatFloat] = None, nbModes: Optional[MuscatIndex] = None):
    """
    Computes a truncated singular value decomposition of a symmetric definite
    matrix in scipy.sparse.csr format. Only the lower triangular part needs
    to be defined

    Parameters
    ----------
    matrix : scipy.sparse.csr
        the input matrix
    epsilon : float
        the truncation tolerance, determining the number of keeps eigenvalues
    nbModes : int
        the number of keeps eigenvalues

    Returns
    -------
    np.ndarray
        kept eigenvalues, of size (numberOfEigenvalues)
    np.ndarray
        kept eigenvectors, of size (numberOfEigenvalues, numberOfSnapshots)
    """

    if epsilon != None and nbModes != None:  # pragma: no cover
        raise RuntimeError("cannot specify both epsilon and nbModes")

    eigenValues, eigenVectors = np.linalg.eigh(matrix, UPLO="L")

    idx = eigenValues.argsort()[::-1]
    eigenValues = eigenValues[idx]
    eigenVectors = eigenVectors[:, idx]

    if nbModes is None:
        if epsilon == None:
            nbModes = matrix.shape[0]
        else:
            nbModes = MuscatIndex(0)
            bound = (epsilon ** 2) * eigenValues[0]
            for e in eigenValues:
                if e > bound:
                    nbModes += 1
            id_max2: MuscatIndex = MuscatIndex(0)
            bound = (1 - epsilon ** 2) * np.sum(eigenValues)
            temp = 0
            for e in eigenValues:
                temp += e
                if temp < bound:
                    id_max2 += 1  # pragma: no cover

            nbModes = max(nbModes, id_max2)

    if nbModes > matrix.shape[0]:
        Warning("nbModes taken to max possible value of "+str(matrix.shape[0])+" instead of provided value "+str(nbModes))
        nbModes = matrix.shape[0]

    index = np.where(eigenValues < 0)
    if len(eigenValues[index]) > 0:
        if index[0][0] < nbModes:
            # print(nbModes, index[0][0])
            Warning("removing numerical noise from eigenvalues, nbModes is set to "+str(index[0][0])+" instead of "+str(nbModes))
            nbModes = index[0][0]

    return eigenValues[0:nbModes], eigenVectors[:, 0:nbModes]





def CheckIntegrity(GUI: bool = False):

    Mat = np.arange(9).reshape(3, 3)
    Mat[np.triu_indices(3, 1)] = 0.0

    ref = (np.array([16.01240515]), np.array([[-0.38252793], [-0.53464575], [-0.7535425]]))

    res = TruncatedSVDSymLower(Mat)
    np.testing.assert_almost_equal(res[0], ref[0])
    np.testing.assert_almost_equal(res[1], ref[1])

    res = TruncatedSVDSymLower(Mat, 1.0e-6)
    np.testing.assert_almost_equal(res[0], ref[0])
    np.testing.assert_almost_equal(res[1], ref[1])

    res = TruncatedSVDSymLower(Mat, nbModes=5)
    np.testing.assert_almost_equal(res[0], ref[0])
    np.testing.assert_almost_equal(res[1], ref[1])

    res = TruncatedSVDSymLower(Mat, nbModes=11)
    np.testing.assert_almost_equal(res[0], ref[0])
    np.testing.assert_almost_equal(res[1], ref[1])

    return "ok"



if __name__ == '__main__':
    print(CheckIntegrity())  # pragma: no cover