# -*- 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 numpy as np
from scipy import sparse
try:
from mpi4py import MPI
except ImportError:# pragma: no cover
print("MPI capabilities not available")
from BasicTools.Helpers.ProgressBar import printProgressBar
[docs]def ScalarProductXY(snapshots1, snapshots2, massMatrix = None):
"""
Constructs the scalar product matrix of two sets of snapshots
Parameters
----------
snapshots1 : np.ndarray
of size (numberOfSnapshots1, numberOfDofs)
snapshots2 : np.ndarray
of size (numberOfSnapshots2, numberOfDofs)
massMatrix : None or scipy.sparse.csr
the sparse matrix of the L2 scalar product, of size (numberOfDofs, numberOfDofs)
Returns
-------
np.ndarray
of size (numberOfSnapshots1, numberOfSnapshots2)
"""
numberOfSnapshots1 = snapshots1.shape[0]
numberOfSnapshots2 = snapshots2.shape[0]
if massMatrix is None:
massMatrix = sparse.eye(snapshots1.shape[1])
scalarProductMatrix = np.zeros((numberOfSnapshots1, numberOfSnapshots2))
for i, snapshot1 in enumerate(snapshots1):
matVecProduct = massMatrix.dot(snapshot1)
for j, snapshot2 in enumerate(snapshots2):
#if j <= i and j < numberOfSnapshots:
scalarProductMatrix[i, j] = np.dot(matVecProduct, snapshot2)
mpiScalarProductMatrix = np.zeros((numberOfSnapshots1, numberOfSnapshots2))
MPI.COMM_WORLD.Allreduce([scalarProductMatrix, MPI.DOUBLE], [mpiScalarProductMatrix, MPI.DOUBLE])
return mpiScalarProductMatrix
[docs]def ScalarProductXX(snapshots, massMatrix = None, verbose = False):
"""
Constructs the scalar product matrix of two sets of snapshots
Parameters
----------
snapshots : np.ndarray
of size (numberOfSnapshots, numberOfDofs)
massMatrix : None or scipy.sparse.csr
the sparse matrix of the L2 scalar product, of size (numberOfDofs, numberOfDofs)
Returns
-------
np.ndarray
of size (numberOfSnapshots1, numberOfSnapshots2)
"""
numberOfSnapshots = snapshots.shape[0]
if massMatrix is None:
massMatrix = sparse.eye(snapshots.shape[1])
scalarProductMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
if verbose:
printProgressBar(0, int(numberOfSnapshots*(numberOfSnapshots+1)/2), prefix = 'Progress ScalarProductXX:', suffix = 'Complete', length = 50)
count = 0
for i, snapshot1 in enumerate(snapshots):
matVecProduct = massMatrix.dot(snapshot1)
for j, snapshot2 in enumerate(snapshots):
if j <= i and j < numberOfSnapshots:
scalarProductMatrix[i, j] = np.dot(matVecProduct, snapshot2)
count += 1
if verbose:
printProgressBar(count, int(numberOfSnapshots*(numberOfSnapshots+1)/2), prefix = 'Progress ScalarProductXX:', suffix = 'Complete', length = 50)
scalarProductMatrix = scalarProductMatrix + scalarProductMatrix.T
mpiScalarProductMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
MPI.COMM_WORLD.Allreduce([scalarProductMatrix, MPI.DOUBLE], [mpiScalarProductMatrix, MPI.DOUBLE])
return mpiScalarProductMatrix
[docs]def Norm(snapshots, massMatrix = None):
"""
Constructs the norm vector of a set of snapshots
Parameters
----------
snapshots : np.ndarray
of size (numberOfSnapshots, numberOfDofs)
Returns
-------
np.ndarray
of size (numberOfSnapshots)
"""
numberOfSnapshots = snapshots.shape[0]
if massMatrix is None:
massMatrix = sparse.eye(snapshots.shape[1])
normVector = np.zeros(numberOfSnapshots)
for i, snapshot in enumerate(snapshots):
matVecProduct = massMatrix.dot(snapshot)
normVector[i] = np.dot(matVecProduct, snapshot)
return np.sqrt(normVector)
[docs]def SineDissimilarity(snapshots, massMatrix = None, verbose = False):
"""
Constructs the sine dissimilarity matrix of a set of snapshots
Parameters
----------
snapshots : np.ndarray
of size (numberOfSnapshots, numberOfDofs)
Returns
-------
np.ndarray
of size (numberOfSnapshots, numberOfSnapshots)
"""
#import time
#start = time.time()
normalizedSnapshots = snapshots / Norm(snapshots, massMatrix)[:, np.newaxis]
cosineDissimilarity = ScalarProductXX(normalizedSnapshots, massMatrix, verbose)
sineDissimilarity = np.sqrt(1. - np.minimum(1.,cosineDissimilarity**2))
np.fill_diagonal(sineDissimilarity, 0)
#print("d1 =", time.time()-start)
"""from sklearn.metrics.pairwise import cosine_similarity
#start = time.time()
cosineDissimilarity2 = cosine_similarity(snapshots)
sineDissimilarity2 = np.sqrt(1. - np.minimum(1.,cosineDissimilarity2**2))
np.fill_diagonal(sineDissimilarity2, 0)
#print("d2 =", time.time()-start)
print("rel diff =", np.linalg.norm(sineDissimilarity2-sineDissimilarity)/np.linalg.norm(sineDissimilarity2))
print("np.linalg.norm(sineDissimilarity) =", np.linalg.norm(sineDissimilarity))
print("np.linalg.norm(sineDissimilarity2) =", np.linalg.norm(sineDissimilarity2))"""
return sineDissimilarity
[docs]def L2Dissimilarity(snapshots, massMatrix = None, faster = True, verbose = False):
"""
Constructs the L2 dissimilarity matrix of a set of snapshots
Parameters
----------
snapshots : np.ndarray
of size (numberOfSnapshots, numberOfDofs)
Returns
-------
np.ndarray
of size (numberOfSnapshots, numberOfSnapshots)
"""
numberOfSnapshots = snapshots.shape[0]
if faster == True:
scalarProductXX = ScalarProductXX(snapshots, massMatrix, verbose)
l2DissimilarityMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
for i in range(numberOfSnapshots):
for j in range(numberOfSnapshots):
if j <= i and j < numberOfSnapshots:
l2DissimilarityMatrix[i,j] = np.sqrt(scalarProductXX[i,i] + scalarProductXX[j,j] -2*scalarProductXX[i,j])
else:
if massMatrix is None:
massMatrix = sparse.eye(snapshots.shape[1])
if verbose:
printProgressBar(0, int(numberOfSnapshots*(numberOfSnapshots+1)/2), prefix = 'Progress ScalarProduct:', suffix = 'Complete', length = 50)
count = 0
l2DissimilarityMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
for i, snapshot1 in enumerate(snapshots):
for j, snapshot2 in enumerate(snapshots):
if j <= i and j < numberOfSnapshots:
delta = snapshot1-snapshot2
matVecProduct = massMatrix.dot(delta)
l2DissimilarityMatrix[i, j] = np.sqrt(np.dot(matVecProduct, delta))
count += 1
if verbose:
printProgressBar(count, int(numberOfSnapshots*(numberOfSnapshots+1)/2), prefix = 'Progress ScalarProduct:', suffix = 'Complete', length = 50)
l2DissimilarityMatrix = l2DissimilarityMatrix+l2DissimilarityMatrix.T
np.fill_diagonal(l2DissimilarityMatrix, 0)
mpiL2DissimilarityMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
MPI.COMM_WORLD.Allreduce([l2DissimilarityMatrix, MPI.DOUBLE], [mpiL2DissimilarityMatrix , MPI.DOUBLE])
"""from scipy.spatial import distance_matrix
mpiL2DissimilarityMatrix2 = distance_matrix(snapshots,snapshots)
np.fill_diagonal(mpiL2DissimilarityMatrix2, 0)
print("rel diff =", np.linalg.norm(mpiL2DissimilarityMatrix-mpiL2DissimilarityMatrix2)/np.linalg.norm(mpiL2DissimilarityMatrix2))
print("np.linalg.norm(mpiL2DissimilarityMatrix) =", np.linalg.norm(mpiL2DissimilarityMatrix))
print("np.linalg.norm(mpiL2DissimilarityMatrix2) =", np.linalg.norm(mpiL2DissimilarityMatrix2))"""
return mpiL2DissimilarityMatrix
if __name__ == "__main__":# pragma: no cover
from genericROM import RunTestFile
RunTestFile(__file__)