# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE', which is part of this source code package.
#
#
import os
try:
from mpi4py import MPI
except ImportError:# pragma: no cover
print("MPI capabilities not available")
import numpy as np
from scipy import sparse
from Mordicus.BasicAlgorithms import SVD as SVD
#REMARK: do not work very well for low values of tolerance
[docs]def CompressData(
collectionProblemData, solutionName, tolerance = None, snapshotCorrelationOperator = None, snapshots = None, compressSolutions = False, nbModes = None
):
"""
Computes a reducedOrderBasis using the SnapshotPOD algorithm, from the
snapshots contained in the iterator snapshotsIterator, which a correlation
operator between the snapshots defined by the matrix
snapshotCorrelationOperator, with tolerance as target accuracy of the data
compression. If a reducedOrderBasis prexists, it is updated using the
snapshots from the solutions.
Parameters
----------
collectionProblemData : CollectionProblemData
input collectionProblemData containing the data
solutionName : str
name of the solutions from which snapshots are taken
tolerance : float, cannot be provided with nbModes
target accuracy of the data compression
snapshotCorrelationOperator : scipy.sparse.csr_matrix, optional
correlation operator between the snapshots
snapshots : np.ndarray, optional
of size (nbSnapshots, numberOfDofs): snapshots of the solutions to
compress
compressSolutions : bool, optional
True to compresse solutions using the constructed reducedOrderBasis
nbModes : int, cannot be provided with tolerance
the number of keps snapshot POD modes
"""
assert isinstance(solutionName, str)
if tolerance == None and nbModes == None:# pragma: no cover
raise("must specify epsilon or nbModes")
if tolerance != None and nbModes != None:# pragma: no cover
raise("cannot specify both epsilon and nbModes")
if snapshots is None:
snapshots = collectionProblemData.GetSnapshots(solutionName)
if snapshotCorrelationOperator is None:
snapshotCorrelationOperator = sparse.eye(snapshots.shape[1])
numberOfSnapshots = snapshots.shape[0]
previousReducedOrderBasis = collectionProblemData.GetReducedOrderBasis(solutionName)
if previousReducedOrderBasis is None:
correlationMatrix = np.zeros((numberOfSnapshots,numberOfSnapshots))
for i, snapshot1 in enumerate(snapshots):
matVecProduct = snapshotCorrelationOperator.dot(snapshot1)
for j, snapshot2 in enumerate(snapshots):
if j <= i and j < numberOfSnapshots:
correlationMatrix[i, j] = np.dot(matVecProduct, snapshot2)
mpiReducedCorrelationMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
MPI.COMM_WORLD.Allreduce([correlationMatrix, MPI.DOUBLE], [mpiReducedCorrelationMatrix, MPI.DOUBLE])
eigenValuesRed, eigenVectorsRed = SVD.TruncatedSVDSymLower(mpiReducedCorrelationMatrix, tolerance)
nbePODModes = eigenValuesRed.shape[0]
changeOfBasisMatrix = np.zeros((nbePODModes,numberOfSnapshots))
for j in range(nbePODModes):
changeOfBasisMatrix[j,:] = eigenVectorsRed[:,j]/np.sqrt(eigenValuesRed[j])
reducedOrderBasis = np.dot(changeOfBasisMatrix,snapshots)
collectionProblemData.AddReducedOrderBasis(solutionName, reducedOrderBasis)
localGamma = np.dot(reducedOrderBasis, snapshotCorrelationOperator.dot(snapshots.T))
globalGamma = np.zeros(localGamma.shape)
MPI.COMM_WORLD.Allreduce([localGamma, MPI.DOUBLE], [globalGamma, MPI.DOUBLE])
collectionProblemData.AddDataCompressionData(solutionName, globalGamma)
else:
for snap in snapshots:
snap.shape = (1,snap.shape[0])
gammaNPrevInd = collectionProblemData.GetDataCompressionData(solutionName)
localGammaNNewInd = np.dot(previousReducedOrderBasis, snapshotCorrelationOperator.dot(snap.T))
gammaNNewInd = np.zeros(localGammaNNewInd.shape)
MPI.COMM_WORLD.Allreduce([localGammaNNewInd, MPI.DOUBLE], [gammaNNewInd, MPI.DOUBLE])
snap = snap - np.einsum('kl,km->ml', previousReducedOrderBasis, gammaNNewInd, optimize = True)
localNorms = np.einsum('kl,lk->k', snap, snapshotCorrelationOperator.dot(snap.T), optimize = True)
globalNorms = np.zeros(localNorms.shape)
MPI.COMM_WORLD.Allreduce([localNorms, MPI.DOUBLE], [globalNorms, MPI.DOUBLE])
globalNorms = np.sqrt(globalNorms)
snap = np.divide(snap, globalNorms[:,np.newaxis])
indices = globalNorms > tolerance
snap = snap[indices,:]
globalNorms = globalNorms[indices]
gammaNNewInd = gammaNNewInd[:,indices]
Vinter = np.append(previousReducedOrderBasis, snap, axis=0)
numberOfAddedSnapshots = snap.shape[0]
gammaInterPrevInd1 = np.append(gammaNPrevInd, np.zeros((numberOfAddedSnapshots,gammaNPrevInd.shape[1])), axis = 0)
gammaInterPrevInd2 = np.append(gammaNNewInd, np.diag(globalNorms), axis = 0)
gammaInterLocal = np.append(gammaInterPrevInd1, gammaInterPrevInd2, axis = 1)
gammaInterLocal2 = np.dot(gammaInterLocal, gammaInterLocal.T)
gammaInterGlobal2 = np.zeros(gammaInterLocal2.shape)
MPI.COMM_WORLD.Allreduce([gammaInterLocal2, MPI.DOUBLE], [gammaInterGlobal2, MPI.DOUBLE])
eigenValuesRed, eigenVectorsRed = SVD.TruncatedSVDSymLower(gammaInterGlobal2, tolerance)
previousReducedOrderBasis = np.dot(Vinter.T, eigenVectorsRed).T
gammaNPrevInd = np.dot(eigenVectorsRed.T, gammaInterLocal)
collectionProblemData.AddReducedOrderBasis(solutionName, previousReducedOrderBasis)
collectionProblemData.AddDataCompressionData(solutionName, gammaNPrevInd)
if compressSolutions == True:
collectionProblemData.CompressSolutions(solutionName, snapshotCorrelationOperator)
if __name__ == "__main__":# pragma: no cover
from genericROM import RunTestFile
RunTestFile(__file__)