#Install package
import sys
!{sys.executable} -m pip install numpy
!{sys.executable} -m pip install matplotlib

# import packages
import numpy as np
import matplotlib.pyplot as plt

def generate_translation_snapshots(Ne=400,x0=-100,x1=100,delta=0.001,n_mu=100):
    """
    Ne+1 mesh elements
    mesh starting point x0
    mesh last point x1
    jump width: delta
    number of tested parameters n_mu<Ne
    """
    
    # fixed mesh
    x = ...

    # uniform mu
    # n_mu values between x0 and x1
    mu_values = ...

    snapshots = []

    for mu in mu_values:
        u = ...
        ...
    
    return np.array(snapshots), x, mu_values # should return snapshots, mesh, mu_values


## Test

snapshots,x,_=generate_translation_snapshots(Ne=50)
for u in snapshots:
    plt.plot(x,u)
plt.xlabel("x")
plt.ylabel("u(x)")
plt.tight_layout()
plt.show()

""" OFFLINE """
""" POD """

def Construct_RB(NumberOfSnapshots=100,Ne=100,NumberOfModes=80):
    
    """
    Number of training parameters: NumberOfSnapshots
    Mesh size : Ne=100
    Number of basis functions: NumberOfModes
    """
    
    Snapshots,x,_=generate_translation_snapshots(Ne=Ne)
    NumberOfSnapshots=len(Snapshots)
    print("Number of snapshots:", NumberOfSnapshots)
    
    ## SVD ##
    # you can use (.,.)_l2 or the scaled scalar product 1./|Omega| * (.,.)_l2
    # indication : you can use the previous finite volume TP 
    
    # We first compute the correlation matrix C_ij = (u_i,u_j)_(l2)
    # Then, we compute the eigenvalues/eigenvectors of C 
    #SVD: C eigenVectors=eigenValues eigenVectors
    # sort the eigenvalues
    # retrieve N=NumberOfModes first eigenvalues
    ...
    
    print("eigenvalues: ",EigenValues)

    RIC=... #must be close to 0
    print("Relativ Information Content (must be close to 0): ",RIC)

    ReducedBasis = ...

    # orthogonality test
    """
    for i in range(NumberOfModes):
        MatVecProduct = L2.dot(ReducedBasis[i])
        for j in range(NumberOfModes):    
            test = np.dot(MatVecProduct, ReducedBasis[j])
            print("orthogonal:",test)
    """
    return ReducedBasis,Snapshots,EigenValues

## L2 projection 
ReducedBasis,Snapshots,eigvals=Construct_RB(NumberOfSnapshots=80,Ne=100,NumberOfModes=80)

# POD eigenvalues
plt.figure(figsize=(7, 5))
plt.semilogy(np.arange(1, len(eigvals) + 1), eigvals, 'o-')
plt.xlabel("N")
plt.ylabel("POD eigenvalues")
plt.grid(True)
plt.show()

### test

ReducedBasis,Snapshots,eigvals=Construct_RB(NumberOfSnapshots=80,Ne=100,NumberOfModes=80)
x = np.linspace(-100, 100, 100 + 1)
plt.plot(x,ReducedBasis[0,:])
plt.plot(x,ReducedBasis[1,:])
plt.xlabel("x")
plt.ylabel("Phi[x]")
plt.tight_layout()
plt.show()

def project_snapshot(snapshot,x, ReducedBasis):
    
    """
    l2 projection of a snapshot on V_N spanned by (Phi)_i
    Phi: shape (Ndofs, n)
    should return \sum_i alpha_i Phi_i where alpha_i = (u,Phi_i)_l2
    """
    ...
    return ...

def compute_kolmogorov_decay(Snapshots, x,ReducedBasis):
 
    """
    Snapshots : array of snapshots
    x: mesh
    ReducedBasis: basis functions
    """
    NumberOfSnapshots, Ndofs = Snapshots.shape
    NumberOfModes = ReducedBasis.shape[0]
    #print(np.shape(ReducedBasis))
    
    dmax = np.zeros(NumberOfModes)

    for n in range(1, NumberOfModes + 1): #forloop over n= 1 ... N
        ReducedBasis_n = ...  # shape (Ndofs, n)

        errs = []
        for i in range(NumberOfSnapshots): #for each snapshots
            ...
            # compute the error ||u-P_Nu||_l2
            err_l2 = ...
            errs.append(err_l2)

        errs = np.array(errs)
        
        dmax[n - 1] = np.max(errs)

    return dmax


# -------------------------
# Test
# -------------------------
NumberOfSnapshots = 80
Ne = 100
NumberOfModes = 40

ReducedBasis,Snapshots,_= Construct_RB(NumberOfSnapshots=NumberOfSnapshots,Ne=Ne,NumberOfModes=NumberOfModes)

x0=-100
x1=100
x = np.linspace(x0, x1, Ne + 1)
dmax = compute_kolmogorov_decay(Snapshots,x, ReducedBasis)

# log-log estimation of the slope
nvals = np.arange(1, len(dmax) + 1)

#Use np.polyfit(x,y,deg) to fit a polynomial p[0] * x + p[1] to points (nvals, dmax ) in log scale
#
coef = ...
alpha = ...
print(f"Slope ~ n^(alpha) : alpha = {alpha:.3f}")

# Plot
plt.figure(figsize=(7, 5))
plt.loglog(nvals, dmax, 'o-', label=r'$d_n^{\max}$')

plt.loglog(nvals,np.exp(coef[1]) * nvals**coef[0],  '--',label=fr'fit $\sim n^{{{alpha:.2f}}}$')
plt.xlabel("n")
plt.ylabel("Projection error")
plt.grid(True, which="both")
plt.legend()
plt.show()