import sys
#!{sys.executable} -m pip install numpy
#!{sys.executable} -m pip install matplotlib
#!{sys.executable} -m pip install scikit-fem
#!{sys.executable} -m pip install tensorflow

# import packages
import skfem  # for Finite Element Method
import numpy as np
import matplotlib.pyplot as plt 
from sklearn.utils.extmath import randomized_svd
import random
import tensorflow as tf
import itertools

random.seed(10)

""" First we define a mesh for the unit square with the boundary decomposition """
mesh= skfem.MeshTri().refined(6).with_boundaries(
        {                                                                                                                                
            "left": lambda x: x[0] == 0,                                                                                            
            "right": lambda x: x[0] == 1,            
            "down": lambda x: x[1] == 0,                                                                                            
            "up": lambda x: x[1] == 1,     
            
        }                                                                                                                               
)

print(mesh)
mesh

""" Then we assemble the matrices of our model problem """

# Assembling matrices

from skfem.assembly import BilinearForm, LinearForm
from skfem.helpers import grad, dot,div,ddot

@BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))


@BilinearForm
def b1(u, v, w):
    telt=w.telt
    #print(telt)
    mu2=w.mu2
    b=[np.cos(np.pi/mu2*telt),np.sin(np.pi/mu2*telt)]
    return b[0]*grad(u)[0]*v+b[1]*grad(u)[1]*v
    
@BilinearForm
def mass(u, v, _):
    return u * v

@LinearForm
def linearterm(v, w):
    mu3=w.mu3
    mu4=w.mu4
    f=10*np.exp(-((w.x[0]-mu3)**2+(w.x[1]-mu4)**2)/0.07**2)
    return f*v

c=1
""" Solving the problem """
from skfem import *

element = {'u': ElementTriP1()} #[u]= Solution with P1 FEM elements
def visualize_solution(Mesh,solution):
    ## Function that allows us to visualize the solution
    from skfem.visuals.matplotlib import plot
    return plot(Mesh, solution, shading='gouraud', colorbar=True)


def SolveProblem_BDF2(mu, Mesh):
    """ Backward Differentiation Formula (BDF) of order 2 """
    mu1,mu2,mu3,mu4 = mu[0]
    T = 10*np.pi #final time
    dt =  5*np.pi / 10 #time step
    basis = {variable: Basis(Mesh, e) for variable, e in element.items()}  # FEM space
    u_prev = basis['u'].zeros()  # u_n (initial condition)
    t_list=[0] #time steps
    u_list=[u_prev]
    u_curr = u_prev.copy()  # u_{n+1}

    M = mass.assemble(basis['u'])
    laplaceMat = mu1 * laplace.assemble(basis['u'])
    lineartermV = linearterm.assemble(basis['u'], mu3=mu3, mu4=mu4)

    # implicit scheme - order 1 for init
    t = dt
    B1 = b1.assemble(basis['u'], mu2=mu2, telt=t)
    K1 = M / dt + laplaceMat + B1 + c*M #c=1
    lineartermT = lineartermV + M @ u_prev / dt
    u_curr = solve(K1, lineartermT)  # u_{n+1}
    t_list.append(dt)
    u_list.append(u_curr)
    
    # BDF2 iteration
    for t in np.arange(2 * dt, T , dt):
        B2 = b1.assemble(basis['u'], mu2=mu2, telt=t)
        K2 = (3 / (2 * dt)) * M + laplaceMat + B2 + c * M
        lineartermT = lineartermV + (2 / dt) * M @ u_curr - (1 / (2 * dt)) * M @ u_prev

        u_next = solve(K2, lineartermT)  # u_{n+2}
        u_prev, u_curr = u_curr, u_next  # Update
        t_list.append(t)
        u_list.append(u_curr)
    
    return u_list,t_list  #[u_curr], [t] #
velocity,t_list=SolveProblem_BDF2([[0.0045,50,0.2, 0.475]],mesh)
visualize_solution(mesh,velocity[-1]).show()

velocity,t_list=SolveProblem_BDF2([[0.006,35,0.475, 0.475]],mesh)
visualize_solution(mesh,velocity[-1]).show()
print("Number of time stepps", len(velocity))

## FINE MESH
FineMesh = skfem.MeshTri().refined(6)

FineBasis = {variable: Basis(FineMesh, e, intorder=4)
            for variable, e in element.items()} # FEM space

NumberOfNodesFineMesh = FineMesh.p.shape[1]
print("Number of nodes: ",NumberOfNodesFineMesh)
#num_dofs_uFineMesh = basis['u'].doflocs.shape[1] # or np.shape(velocity)[0] for DOFs

# """ TRAINING set """

print("-----------------------------------")
print("        Offline                    ")
print("-----------------------------------")
# The parameter set is given by $[0.002, 0.005] × [30, 70] × [0.4, 0.6]^2$
# and we take for the training set 4*5*5*5 uniformly distributed in each parametric direction
#Nt = 100 time instances

mu1_values = 0.004 # np.linspace(0.002, 0.005, 4)  # 4 valeurs dans [0.002, 0.005]
mu2_values = np.linspace(10, 100, 16)        # 5 valeurs dans [30, 70]
mu3_values = np.linspace(0.4, 0.6, 6)      # 5 valeurs dans [0.4, 0.6]
mu4_values = 0.5 #np.linspace(0.4, 0.6, 5)
mu5_values = t_list # additional parameter = time

print("mu1: ",mu1_values)
print("mu2: ",mu2_values)
print("mu3: ",mu3_values)
print("mu4: ",mu4_values)

#param_combinations = list(itertools.product(mu1_values, mu2_values, mu3_values, mu4_values)) 
param_combinations = list(itertools.product( mu2_values, mu3_values)) 

FineSnapshots=[]
muParam=[]

for mu in param_combinations: 
    #print(mu[0])
    
    mutmp=[mu1_values,mu[0],mu[1],mu4_values]

    solution,_=SolveProblem_BDF2([mutmp],FineMesh)
    mu_full = np.hstack([[mu] * len(velocity) , np.array(mu5_values).reshape((len(velocity),1))])  # we add time as a parameter mu5
    mu_full = mu_full.tolist()
    muParam += mu_full 
    FineSnapshots.append(solution)
    
print("Number of training parameters: ",len(muParam))
print("First parameter:", muParam[0])
print("Shape snapshots matrix: ",np.shape(np.array(FineSnapshots)))

visualize_solution(FineMesh,solution[-1]).show()

## randomSVD on the snapshots
NumberOfModes=64
import time

start_time=time.time()
FineSnapshotsMat=np.concatenate(FineSnapshots)
print("shape snapshots: ", FineSnapshotsMat.shape) #NS, NumberOfDegreesOfFreedrom (NS=Ntrain*NT)

U, Sigma, VT = randomized_svd(np.array(FineSnapshotsMat).transpose(), 
                              n_components=NumberOfModes, 
                              n_iter=5, 
                              random_state=42)  #m=10 (oversampling default value)

ReducedBasis = U[:, :NumberOfModes]  # (Ndof x NumberOfModes)
EigenValues = Sigma[:NumberOfModes] ** 2  
print("shape reduced basis: " , np.shape(ReducedBasis))

end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")

# In view of the DL-ROM strategy, we decompose Mtrain in [Mtrain, Mval] where Mval is of size alpha*NS 
#(and NS=Ntrain*NT, NT beeing the number of time steps) 
# but first a random shuffle

FineSnapshotsMat=np.concatenate(FineSnapshots)
print("shape snapshots: ", FineSnapshotsMat.shape) #NS, NumberOfDegreesOfFreedrom (NS=Ntrain*NT)

#randomly shuffle to decompose in training/testing parts
indices = list(range(len(FineSnapshotsMat))) 
random.shuffle(indices)

FineSnapshots_training = [FineSnapshotsMat[i,:] for i in indices] #R^NS=R^Ntrain*NT

SVector=FineSnapshots_training.copy() # full snapshots
MVector = np.array([muParam[i] for i in indices]) #full parameters

alpha = 0.1  # e.g. 10% for testing

n = MVector.shape[0]
split_point = int((1 - alpha) * n)

#parameters
M1 = MVector[:split_point]   # 90%
M2 = MVector[split_point:]   # 10%

#values
S1 = SVector[:split_point]
S2 = SVector[split_point:]

print("S1 (snapshots for training) shape:", np.shape(S1))
print("M2 (parameters for testing) shape:", np.shape(M2))

print("------------------------------------------------")
print("  Compress training/testing snapshots           ")
print("------------------------------------------------")

VTS1=np.zeros((np.shape(S1)[0],NumberOfModes))
VTS2=np.zeros((np.shape(S2)[0],NumberOfModes))
print(np.shape(S1))
for i,snap in enumerate(S1):#FineSnapshotsMat:
    ExactSolution =snap #size Nh
    CompressedSolutionU= ReducedBasis.transpose()@ExactSolution #snapshot compression
    VTS1[i]=ReducedBasis.transpose()@ExactSolution #snapshot compression
    #ReconstructedCompressedSolution = np.dot(VTS1[i], ReducedBasis.transpose())

for i,snap in enumerate(S2):#FineSnapshotsMat:
    ExactSolution =snap #size Nh
    VTS2[i]=ReducedBasis.transpose()@ExactSolution #snapshot compression
    #ReconstructedCompressedSolution = np.dot(VTS2[i], ReducedBasis.transpose())

print("Compressed size",VTS1.shape)
ReconstructedCompressedSolution = np.dot(VTS2[-1], ReducedBasis.transpose())
visualize_solution(FineMesh,ReconstructedCompressedSolution).show()

## Normalization
"""
## parameters training M1 and test M2 matrices ##
##
"""
print(M1.shape)
Mmax=np.zeros(M1.shape[1])
Mmin=np.zeros(M1.shape[1])

for j in range(M1.shape[1]):
    Mmax[j]=np.max(M1[:,j])
    Mmin[j]=np.min(M1[:,j])

# Training parameters
M1train=np.zeros(np.shape(M1))

for j in range(M1.shape[1]):
    M1train[:,j]=(M1[:,j]-Mmin[j])/(Mmax[j]-Mmin[j])
    
#M1train[:,4]=mu5_values[0] ##if only one time step
M1train=M1train.transpose()


# Testing part
M2test=np.zeros(np.shape(M2))
for j in range(M2.shape[1]):
    M2test[:,j]=(M2[:,j]-Mmin[j])/(Mmax[j]-Mmin[j])

#M2test[:,4]=mu5_values[0]


"""
## snapshots training S1 and test S2 matrices ##
##
"""
Smax=np.max(VTS1)
Smin=np.min(VTS1)
S1train=(VTS1-Smin)/(Smax-Smin)
S1train=S1train.transpose()

S2test=(VTS2-Smin)/(Smax-Smin)

NS1=S1train.shape[1]
print("number of training parameters: ",NS1)
#print(np.shape(S2))

from tensorflow.keras import layers, models, initializers
# 1/ encode the solution with a CAE that goes from R^N -> R^n with n as close as possible to the number of parameters (n<=N)
# 2/ learn to approximate the latent space of the parameters with a DFNN 
# 3/ decode this approximation to approach the original solution

print("shape S1 training:",np.shape(S1train))
print("shape S2 validation:",np.shape(S2test))
print("shape M1 parameters:",np.shape(M1train))
print("shape M2 parameters:",np.shape(M2test))

#Training inputs:
# Optimizers (learning rate = eta)
eta=5e-4
#NS1=S1train.shape[1]
print("Total number of training parameters:", NS1)
N_epochs=2500 # Number of iterations
print("Maximum iterations:", N_epochs)
BatchSize=108 #batch size (allows to save CPU memory) 
print("Number of minibatches:",NS1/BatchSize)
assert NS1%BatchSize==0
NbMiniBatches=int(NS1/BatchSize) # Number of Minibatches
# We can add an initialization of theta_0 from a coarser model

n=M1train.shape[0]
print("latent space dim:",n)

"""================"""
""" Initialization """
"""================"""


def build_encoder(latent_dim=n):
    # Input: need a latent space of dimension n < NumberOfModes 
    # Output: Encoder (CAE)

    # Three conv filters:
    # 1. 32 output pictures where each filter is a 5*5 matrix that scans the picture
    #    It yields 5*5*1=25 weights for one filter and one bias
    #    elu activation function as in the paper "A comprehensive DL-based approach to ROM of nonlinear time-dependent parametrized PDEs"
    # 2. 64 output pictures stride=2 -> jump one pixel at a time
    
    inputs = tf.keras.Input(shape=(8, 8, 1)) # picture of size 8*8 = sqrt(NumberOfModes)*sqrt(NumberOfModes)
    x = layers.Conv2D(32, kernel_size=5, strides=1, padding='same', activation='elu')(inputs)  # 8x8x32
    x = layers.Conv2D(64, kernel_size=5, strides=2, padding='same', activation='elu')(x)       # 4x4x64
    x = layers.Conv2D(128, kernel_size=5, strides=2, padding='same', activation='elu')(x)      # 2x2x128
    x = layers.Flatten()(x)
    x = layers.Dense(128, activation='elu')(x)
    latent = layers.Dense(latent_dim, activation='elu', name='latent_vector')(x)

    return tf.keras.Model(inputs, latent, name='encoder')
    

def build_dfnn(latent_dim=n):
    inputs = tf.keras.Input(shape=(n,)) 
    print(inputs.shape)
    x = tf.keras.layers.Dense(NumberOfModes, activation='elu')(inputs)
    x = tf.keras.layers.Dense(NumberOfModes, activation='elu')(x)
    x = tf.keras.layers.Dense(latent_dim)(x)
    print(x.shape)
    return tf.keras.Model(inputs, x)


def build_decoder(latent_dim=n):
    inputs = tf.keras.Input(shape=(latent_dim,))
    x = layers.Dense(128, activation='elu')(inputs)
    x = layers.Dense(2 * 2 * 128, activation='elu')(inputs)
    x = layers.Reshape((2, 2, 128))(x)  # last encoded part 
    x = layers.Conv2DTranspose(64, kernel_size=5, strides=2, padding='same', activation='elu')(x)   # 4x4x64
    x = layers.Conv2DTranspose(32, kernel_size=5, strides=2,padding='same', activation='elu')(x)   # 8x8x32
    outputs = layers.Conv2DTranspose(1, kernel_size=5, padding='same', activation='linear')(x)# Last layer to obtain a picture 8x8x1

    return tf.keras.Model(inputs, outputs, name='decoder')


encoder = build_encoder(latent_dim=n)
loss_enc = tf.keras.losses.MeanSquaredError()
theta_E_stock = []

dfnn = build_dfnn(latent_dim=n)
loss_dfnn = tf.keras.losses.MeanSquaredError()
theta_DFNN_stock = []

decoder = build_decoder(latent_dim=n)
loss_dec = tf.keras.losses.MeanSquaredError()
#loss_dec2 = lambda y_true, y_pred: 1 - tf.keras.losses.CosineSimilarity()(y_true, y_pred) #can try other errors
theta_D_stock = []

optimizer = tf.keras.optimizers.Adam(learning_rate=eta)

"""========================"""
""" Parameter optimization """
"""========================"""

# Step 2: reshape in picture frame for CAE
print(S2test.shape)
S2_reshaped=tf.reshape(S2test,((-1,int(np.sqrt(NumberOfModes)),int(np.sqrt(NumberOfModes)),1))) # for encoder
S2_batch=tf.reshape(S2test, (-1, NumberOfModes))  #for loss

ne=0
while ne<=N_epochs:
    #if ne==500 and total_loss>=total_loss_prec: #early-stopping criterion
    #    print(f"early stopping...ne=={ne}==500 and {total_loss}>={total_loss_prec}" )
    #   break
    #total_loss_prec=total_loss
    """==========="""
    ## Training part 
    """==========="""
    for k in range(NbMiniBatches):
        # Iteration (ne, k)
        iteration = NbMiniBatches * ne + k
        
        # Step 1: sample minibatch [M_batch,S_batch] from M1, S1
        start = k * BatchSize
        end = start + BatchSize
        M1_batch = M1train.T[start:end,:]  # shape (n, 120)
        S1_batch = S1train.T[start:end,:]  # shape (N, 120)

        # Step 2: reshape in picture frame for CAE
        S1_batch_reshaped = tf.reshape(tf.cast(S1_batch, tf.float32),(BatchSize, int(np.sqrt(NumberOfModes)), int(np.sqrt(NumberOfModes)), 1))          
            
        # Step 3: encoder from S1_batch [Convolutional autoencoder (CAE)] + DFNN + decoder
        # with tensorflow and an activation function ELU
         
        with tf.GradientTape() as tape: # for automatic differentiation
            # step a/
            S1tilde_n_batch = encoder(S1_batch_reshaped,training=True)   #latent space that we want to approximate of shape (BatchSize, latent_dim)

            # step b/
            latent_pred = dfnn(M1_batch, training=True)  # function to approximate the latent space. shape: (BatchSize, latent_dim)       

            # step c/
            S_reconstructed = decoder(latent_pred, training=True) 
            St1_batch = tf.reshape(S_reconstructed, (BatchSize, NumberOfModes))
            
            # Step 4: Losses
            loss_int = loss_dfnn(S1tilde_n_batch, latent_pred)  # latent space loss (S1tilde_n_batch= true latent space)
            loss_rec = loss_dec(S1_batch, St1_batch)      # reconstruction loss 
            total_loss = loss_int+loss_rec #alpha * loss_int + beta * loss_rec


        # Weights and biases
        weights_biases = (encoder.trainable_weights + dfnn.trainable_weights + decoder.trainable_weights)

        # Total loss gradient w.r.t. parameters (weights/biases): tape=automatic diff
        grads = tape.gradient(total_loss, weights_biases)

        # Step 5: Update weights/biases with one optimizer 
        optimizer.apply_gradients(zip(grads, weights_biases))

       
    """==========="""
    ## Validation part 
    """==========="""
   

    # we store the weights_biases 
    if ne % 1000 == 0:
        theta_E_stock.append(encoder.get_weights())
        theta_DFNN_stock.append(dfnn.get_weights())
        theta_D_stock.append(decoder.get_weights())


    # Redo step 3: with encoder from S2_batch [Convolutional autoencoder (CAE)]
    # with tensorflow and an activation function ELU

    # step a/
    S2tilde_n = encoder(S2_reshaped,training=True)   #latent space that we want to approximate of shape (BatchSize, latent_dim)
    
    # step b/
    latent_pred = dfnn(M2test, training=True)  # function to approximate the latent space. shape: (BatchSize, latent_dim)

    # step c/
    S_reconstructed = decoder(latent_pred, training=True) 
    St2 = tf.reshape(S_reconstructed, (-1, NumberOfModes))
     
    # Step 4: Losses #for early criterion
    loss_int = loss_dfnn(S2tilde_n, latent_pred)  # latent space loss (S2tilde_n= true latent space)
    loss_rec = loss_dec(S2_batch, St2)      # reconstruction loss 
    total_loss =   loss_int + loss_rec
         
    if ne % 500 == 0:
        print(f"loss_recStep {ne}/{N_epochs}. Loss testing : {total_loss.numpy():.6f}")

    ne+=1

        
ExactSolution=S2_batch[10]
ExactSolution=ExactSolution*(Smax-Smin)+Smin
ExactSolution = np.dot(ExactSolution, ReducedBasis.transpose())
        
CompressedSolutionU =St2[10] 
CompressedSolutionU=CompressedSolutionU*(Smax-Smin)+Smin
ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
   
#print(np.linalg.norm((S2_batch[10].numpy()-St2[10].numpy())))
visualize_solution(FineMesh,ExactSolution).show()
visualize_solution(FineMesh,ReconstructedCompressedSolution).show()

# validation visualization
print(S2test.shape)
S2_batch=tf.reshape(S2test, (-1, NumberOfModes))      
latent_pred = dfnn(M2test, training=True)  # function to approximate the latent space. shape: (BatchSize, latent_dim)       
S_reconstructed = decoder(latent_pred, training=True) 
St2 = tf.reshape(S_reconstructed, (-1, NumberOfModes))
      
ExactSolution=S2_batch[25] 
ExactSolution=ExactSolution*(Smax-Smin)+Smin
ExactSolution = np.dot(ExactSolution, ReducedBasis.transpose())

CompressedSolutionU =St2[25]#St[i] #size Nh
CompressedSolutionU=CompressedSolutionU*(Smax-Smin)+Smin
ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())

visualize_solution(FineMesh,ExactSolution).show()
visualize_solution(FineMesh,ReconstructedCompressedSolution).show()

""" TESTING set """

print("-----------------------------------")
print("        Online                     ")
print("-----------------------------------")
# The parameter set is given by $[0.002, 0.005] × [30, 70] × [0.4, 0.6]^2$
# Nt = 20 time instances
# testing values (different from the training ones)

mu1_values = 0.004 # np.linspace(0.002, 0.005, 4)  # 4 valeurs dans [0.002, 0.005]
mu2_values = np.linspace(35, 65, 4)        # 5 valeurs dans [30, 70]
mu3_values = np.linspace(0.425, 0.575, 4)       # 5 valeurs dans [0.4, 0.6]
mu4_values = 0.5 #np.linspace(0.4, 0.6, 5)
mu5_values = t_list # additional parameter = time

print("mu1: ",mu1_values)
print("mu2: ",mu2_values)
print("mu3: ",mu3_values)
print("mu4: ",mu4_values)

param_combinations = list(itertools.product( mu2_values, mu3_values)) 

FineSnapshots_testing=[]
muParam=[]

for mu in param_combinations: 
    mutmp=[mu1_values,mu[0],mu[1],mu4_values]
    
    solution,_=SolveProblem_BDF2([mutmp],FineMesh)
    mu_full = np.hstack([[mu] * len(solution) , np.array(mu5_values).reshape((len(solution),1))])  # we add time as a parameter mu5
    mu_full = mu_full.tolist()
    muParam += mu_full 
    FineSnapshots_testing.append(solution)


visualize_solution(FineMesh,solution[-1]).show()

FineSnapshotsMat=np.concatenate(FineSnapshots_testing)
print("shape snapshots: ", FineSnapshotsMat.shape) #NS, NumberOfDegreesOfFreedrom (NS=Ntest*NT)

indices = list(range(len(FineSnapshotsMat))) 
FineSnapshots_testing = [FineSnapshotsMat[i,:] for i in indices] #R^NS=R^Ntest*NT

S2 = FineSnapshots_testing.copy() # full snapshots
M2 = np.array([muParam[i] for i in indices]) #full parameters

print("S2 (snapshots for testing) shape:", np.shape(S2))
print("M2 (parameters for testing) shape:", np.shape(M2))

print("-----------------------------------")
print("        Data  Compression          ")
print("-----------------------------------")

VTS2=np.zeros((np.shape(S2)[0],NumberOfModes))

for i,snap in enumerate(S2):#FineSnapshotsMat:
    ExactSolution =snap #size Nh
    VTS2[i]=ReducedBasis.transpose()@ExactSolution #snapshot compression
    #ReconstructedCompressedSolution = VTS2[i]

CompressedSolutionU=VTS2[1]
ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
print("Compressed size",VTS2.shape)
visualize_solution(FineMesh,ReconstructedCompressedSolution).show()

## Normalization
"""
## parameters test M2 matrix ##
##
"""
# Testing part
M2test=np.zeros(np.shape(M2))
for j in range(M2.shape[1]):
    M2test[:,j]=(M2[:,j]-Mmin[j])/(Mmax[j]-Mmin[j])

#M2test[:,4]=mu5_values[0] #if one time step

"""
## snapshots test S2 matrix ##
##
"""
S2test=(VTS2-Smin)/(Smax-Smin)

"""==========="""
## Testing part 
"""==========="""

print(S2test.shape)
S2_batch=tf.reshape(S2test, (-1, NumberOfModes))      
latent_pred = dfnn(M2test, training=True)  # function to approximate the latent space. shape: (BatchSize, latent_dim)       
S_reconstructed = decoder(latent_pred, training=True) 
St2 = tf.reshape(S_reconstructed, (-1, NumberOfModes))
     

ExactSolution=S2_batch[10] 
ExactSolution=ExactSolution*(Smax-Smin)+Smin
ExactSolution = np.dot(ExactSolution, ReducedBasis.transpose()) #true compressed

CompressedSolutionU =St2[10]#St[i] #size Nh
CompressedSolutionU=CompressedSolutionU*(Smax-Smin)+Smin
ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
   
visualize_solution(FineMesh,FineSnapshots_testing[10]).show() #true solution FineSnapshots_testing[i]
visualize_solution(FineMesh,ExactSolution).show() #true compressed
visualize_solution(FineMesh,ReconstructedCompressedSolution).show() #POD-DL-ROM approx

"""==================="""
## Testing part accuracy
"""==================="""

L2compressionErrors=[]
L2=mass.assemble(FineBasis['u'],nu=1)

print(np.shape(VTS2))
for i,snap in enumerate(St2[:10]):
    ExactSolution=FineSnapshots_testing[i]
    
    CompressedSolutionU=snap*(Smax-Smin)+Smin #POD_DL_ROM compressed data
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)

    L2compressionErrors.append(relL2Error)
  
print("L2 averaged compression error =", np.mean(L2compressionErrors))

visualize_solution(FineMesh,ExactSolution).show()
visualize_solution(FineMesh,ReconstructedCompressedSolution).show()