ArnauMiro / pyLowOrder

High performance parallel Low Order Modelling library by BSC - UPM
MIT License
8 stars 3 forks source link

Unexpected behaviour in POD reconstruction #38

Closed Miguel-San closed 5 months ago

Miguel-San commented 5 months ago

Hola,

Estoy obteniendo un comportamiento inesperado al realizar reconstrucciones mediante pyLOM.POD.reconstruct(). En ocasiones la matriz está mal reconstruida. Sin embargo, si haces un flatten a la matriz a reconstruir, la devuelves a su forma original (de tal forma que la nueva matriz obtenida tiene exactamente los mismos valores y forma que la matriz original) y realizas el POD sobre esta nueva matriz, la reconstrucción obtenida es la correcta.

A continuación doy un ejemplo:

Minimal example

import numpy as np
import matplotlib.pyplot as plt
import pyLOM

def generate_X():
    x = np.linspace(0, 1, 100)
    y = np.linspace(0, 1, 200)
    z = np.linspace(0, 1, 300)
    XX, YY, ZZ = np.meshgrid(x, y, z)
    F = np.sin(2*np.pi*(XX + ZZ)) + np.cos(2*np.pi*(2*YY + ZZ))
    return F

# Generate synthetic database
X = generate_X()

########################################################################################################
# Using the whole database. This works fine, reconstructed matrix resembles original matrix
########################################################################################################

X_rsh = X.reshape([-1, X.shape[-1]])

PSI, S, V = pyLOM.POD.run(X_rsh, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
X_rsh_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)
print("RMSE for whole database =", np.sqrt(np.mean((X_rsh-X_rsh_rec)**2)))

########################################################################################################
# Using just the first "snapshot" (as if we are doing data compression). Reconstructed matrix is wrong
########################################################################################################

X_0 = X[:, :, 0]

PSI, S, V = pyLOM.POD.run(X_0, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
X_0_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)
print("RMSE for just first snapshot =", np.sqrt(np.mean((X_0-X_0_rec)**2)))

fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].pcolor(X_0)
axs[0].set_title("Original matrix")
axs[1].pcolor(X_0_rec)
axs[1].set_title("Reconstructed matrix")
plt.show()

########################################################################################################
# Workaround the issue by redefining the shape of the matrix
########################################################################################################

X_redef = X_0.flatten()
X_redef = X_redef.reshape(X_0.shape)
print("RMSE(X_0, X_redefined) =", np.sqrt(np.mean((X_0-X_redef)**2))) # Making sure that these two are the same matrix

PSI, S, V = pyLOM.POD.run(X_redef, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
X_redef_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)

print("RMSE for first snapshot with redefined shape =", np.sqrt(np.mean((X_0-X_redef_rec)**2))) # Now pyLOM works fine!

fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].pcolor(X_0)
axs[0].set_title("Original matrix")
axs[1].pcolor(X_redef_rec)
axs[1].set_title("Reconstructed (redefined) matrix")
plt.show()

Outputs

Shape of X:  (200, 100, 300)
RMSE for whole database = 8.876380209461154e-16
RMSE for just first snapshot = 1.532451748953142
RMSE(X_0, X_redefined) = 0.0
RMSE for first snapshot with redefined shape = 2.82832281954419e-16

fig1 fig2

Installation procedure

Para la instalación estoy siguiendo los pasos de la wiki. El único cambio que estoy haciendo es en la versión de MKL instalada. En el archivo Deps/oneAPI/install_mkl.sh estoy cambiando la línea SRC=https://registrationcenter-download.intel.com/akdlm/irc_nas/17977/l_BaseKit_p_${VERS}.sh por SRC=https://registrationcenter-download.intel.com/akdlm/IRC_NAS/fdc7a2bc-b7a8-47eb-8876-de6201297144/l_BaseKit_p_2024.1.0.596.sh para obtener la última versión (por lo visto Intel no deja descargar de forma gratuita versiones anteriores). Link obtenido aquí

Adjunto también el output del comando conda list:

# packages in environment at /home/miguel/miniconda3/envs/AIRML:
#
# Name                    Version                   Build  Channel
_libgcc_mutex             0.1                        main  
_openmp_mutex             5.1                       1_gnu  
ca-certificates           2024.3.11            h06a4308_0  
certifi                   2024.2.2                 pypi_0    pypi
charset-normalizer        3.3.2                    pypi_0    pypi
contourpy                 1.1.1                    pypi_0    pypi
cycler                    0.12.1                   pypi_0    pypi
cython                    3.0.10                   pypi_0    pypi
ensight-reader            0.9.0                    pypi_0    pypi
exceptiongroup            1.2.1                    pypi_0    pypi
fonttools                 4.52.4                   pypi_0    pypi
h5py                      3.11.0                   pypi_0    pypi
idna                      3.7                      pypi_0    pypi
importlib-resources       6.4.0                    pypi_0    pypi
iniconfig                 2.0.0                    pypi_0    pypi
kiwisolver                1.4.5                    pypi_0    pypi
ld_impl_linux-64          2.38                 h1181459_1  
libffi                    3.4.4                h6a678d5_1  
libgcc-ng                 11.2.0               h1234567_1  
libgomp                   11.2.0               h1234567_1  
libstdcxx-ng              11.2.0               h1234567_1  
matplotlib                3.7.5                    pypi_0    pypi
mpi4py                    3.1.6                    pypi_0    pypi
ncurses                   6.4                  h6a678d5_0  
nfft                      0.1                      pypi_0    pypi
numpy                     1.24.4                   pypi_0    pypi
openssl                   3.0.13               h7f8727e_2  
packaging                 24.0                     pypi_0    pypi
pillow                    10.3.0                   pypi_0    pypi
pip                       24.0             py38h06a4308_0  
platformdirs              4.2.2                    pypi_0    pypi
pluggy                    1.5.0                    pypi_0    pypi
pooch                     1.8.1                    pypi_0    pypi
pyloworder                1.3.5                    pypi_0    pypi
pyparsing                 3.1.2                    pypi_0    pypi
pytest                    8.2.1                    pypi_0    pypi
python                    3.8.19               h955ad1f_0  
python-dateutil           2.9.0.post0              pypi_0    pypi
pyvista                   0.43.8                   pypi_0    pypi
readline                  8.2                  h5eee18b_0  
requests                  2.32.2                   pypi_0    pypi
scipy                     1.10.1                   pypi_0    pypi
scooby                    0.10.0                   pypi_0    pypi
setuptools                69.5.1           py38h06a4308_0  
six                       1.16.0                   pypi_0    pypi
sqlite                    3.45.3               h5eee18b_0  
tk                        8.6.14               h39e8969_0  
tomli                     2.0.1                    pypi_0    pypi
urllib3                   2.2.1                    pypi_0    pypi
vtk                       9.3.0                    pypi_0    pypi
wheel                     0.43.0           py38h06a4308_0  
xz                        5.4.6                h5eee18b_1  
zipp                      3.19.0                   pypi_0    pypi
zlib                      1.2.13               h5eee18b_1 

Gracias de antemano!

ArnauMiro commented 5 months ago

Hola,

El uso de pyLOM es un poco diferente al que propones. Por ello el ejemplo2 que pones conceptualmente esta mal respecto a la que la tool espera y el "workaround" que defines es la forma correcta de hacerlo.

El hecho es que pyLOM está disenyado para aceptar una matriz tipo: Screenshot_20240614_093426 por lo que en la primera dimensión siempre se trata de la dimensión en espacio y la segunda en tiempo. Si hay más de una variable esas se apilan a continuación de la matriz X. Por ejemplo en https://github.com/ArnauMiro/pyLowOrder/blob/126b6da95805204afd4902555626e74cb3e2854c/Converters/alya2h5_channel.py#L37-L59 se muestra como se hace.

En el caso que propones, X_0 tiene dimensiones de (100,200) por lo que pyLOM entiende que 100 són los puntos en espacio y 200 en tiempo mientras que (20000,1) entiende que 20000 són los puntos en espacio y 1 en tiempo.

El ejemplo mínimo, adaptado a la forma de trabajar de pyLOM sería:

import numpy as np
import matplotlib.pyplot as plt
import pyLOM

def generate_X():
    x = np.linspace(0, 1, 100)
    y = np.linspace(0, 1, 200)
    z = np.linspace(0, 1, 300)
    XX, YY, ZZ = np.meshgrid(x, y, z)
    F = np.sin(2*np.pi*(XX + ZZ)) + np.cos(2*np.pi*(2*YY + ZZ))
    return F

# Generate synthetic database
X = generate_X()
print('Shape of synthetic database:',X.shape) # (200, 100, 300)

#########################################################################################################
## However, pyLOM takes a database of the shape (spatial_dims,time_dims)
## thus we need to reshape this matrix accordingly.
#########################################################################################################
X_rsh = X.reshape([-1, X.shape[-1]])
print('Shape of reshaped database:',X_rsh.shape) # (200*100 = 20000, 300)

#########################################################################################################
## POD using the whole database.
#########################################################################################################
PSI, S, V = pyLOM.POD.run(X_rsh, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
X_rsh_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)
print('Shape of reconstructed database:',X_rsh_rec.shape) # (200*100 = 20000, 300)
print("RMSE for whole database =", pyLOM.math.RMSE(X_rsh_rec,X_rsh))

fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].pcolor(X_rsh[:,0].reshape((200,100)))
axs[0].set_title("Original matrix")
axs[1].pcolor(X_rsh_rec[:,0].reshape((200,100)))
axs[1].set_title("Reconstructed matrix")
plt.show()

########################################################################################################
# POD using just the first "snapshot" (as if we are doing data compression).
########################################################################################################
X_rsh_2 = X_rsh[:,:1].copy()
PSI, S, V = pyLOM.POD.run(X_rsh_2, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
X_rsh_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)
print('Shape of reconstructed database:',X_rsh_rec.shape) # (200*100 = 20000, 0)
print("RMSE for whole database =", pyLOM.math.RMSE(X_rsh_rec,X_rsh_2))

fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].pcolor(X_rsh_2[:,0].reshape((200,100)))
axs[0].set_title("Original matrix")
axs[1].pcolor(X_rsh_rec[:,0].reshape((200,100)))
axs[1].set_title("Reconstructed matrix")
plt.show()

Aprobecho y añado ese ejemplo en los de pyLOM, ya que puede servir en un futuro.

Miguel-San commented 5 months ago

Hola @ArnauMiro,

Muchas gracias por tu respuesta.

Creo que no termino de entender por qué mi forma de plantearlo es incorrecta. Imagínate que tengo una base de datos X de una única dimensión espacial + la dimensión temporal, y pongamos que el mallado espacial tiene 200 puntos y el temporal 100, de tal forma que mi base de datos X es de 200x100. Si ejecuto POD sobre esta base de datos, usando X como matriz de snapshots, obtengo los resultados correctos.

Sin embargo, pongamos ahora que tengo otra base de datos Y que comprende cinco de las anteriores X, de tal forma que mi base de datos Y es de 200x100x5. Si construyo una matriz de snapshots Q a partir de la primera de las X (de tal forma que Q es de 200x100), y ejecuto POD sobre Q, obtengo resultados erróneos.

Además, en el primero de los casos (usando X como matriz de snapshots) el RMSE que obtiene pyLOM.math.RMSE y el que calculo yo usando directamente numpy coinciden. Sin embargo, en el segundo caso (usando Q como matriz de snapshots), obtengo valores completamente diferentes de RMSE.

Adjunto un mini-ejemplo con esto que comento y los resultados que obtengo.

Tampoco entiendo por qué el workaround que puse en el primer ejemplo es la forma correcta de hacerlo. Lo que hago es:

X_redef = X_0.flatten()
X_redef = X_redef.reshape(X_0.shape)

X_redef y X_0 en teoría representan la misma matriz (mismos valores, misma dimensión), pero ejecutando POD con X_redef los resultados son correctos, y con X_0 no.

Gracias de antemano!

Ejemplo 2

import numpy as np
import matplotlib.pyplot as plt
import pyLOM

def generate_X():
    x = np.linspace(0, 1, 200)
    t = np.linspace(0, 1, 100)
    TT, XX = np.meshgrid(t, x)
    F = np.sin(2*np.pi*(XX + TT)) + np.cos(2*np.pi*(5*XX + 2*TT))
    return F

# Generate synthetic database
X = generate_X()
print('Shape of synthetic database:',X.shape) # (200, 100)

#########################################################################################################
## POD using X as database.
#########################################################################################################
PSI, S, V = pyLOM.POD.run(X, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
X_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)
print('Shape of reconstructed database:',X_rec.shape) # (200, 100)
print("RMSE for whole database (pyLOM) =", pyLOM.math.RMSE(X_rec,X))
print("RMSE for whole database (numpy) =", np.sqrt(np.mean((X_rec-X)**2)))

fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].pcolor(X)
axs[0].set_title("Original matrix")
axs[1].pcolor(X_rec)
axs[1].set_title("Reconstructed matrix")
fig.suptitle("Using X database")
plt.show()

# ########################################################################################################
# # POD using Q (the first X of the Y database)
# ########################################################################################################
Y = np.stack([X.copy()]*5, axis=-1) # Y database consists of five X stacked
Q = Y[:, :, 0] # Snapshot matrix: First X of database Y
print("Shape of Y database: ", Y.shape) # (200, 100, 5)
print("Shape of Q (first X of Y database): ", Q.shape) # (200, 100)
PSI, S, V = pyLOM.POD.run(Q, remove_mean=False)
PSI_trunc, S_trunc, V_trunc = pyLOM.POD.truncate(PSI,S,V,r=1e-6)
Q_rec = pyLOM.POD.reconstruct(PSI_trunc, S_trunc, V_trunc)
print('Shape of reconstructed database:',Q_rec.shape) # (200, 100)
print("RMSE for whole database (pyLOM) =", pyLOM.math.RMSE(Q_rec,Q))
print("RMSE for whole database (numpy) =", np.sqrt(np.mean((Q_rec-Q)**2)))

fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].pcolor(Q)
axs[0].set_title("Original matrix")
axs[1].pcolor(Q_rec)
axs[1].set_title("Reconstructed matrix")
fig.suptitle("Using first X of Y database")
plt.show()

Output ejemplo 2

Shape of synthetic database: (200, 100)
Shape of reconstructed database: (200, 100)
RMSE for whole database (pyLOM) = 1.623651262933147e-15
RMSE for whole database (numpy) = 1.6236512629331486e-15

fig_1

Shape of Y database:  (200, 100, 5)
Shape of Q (first X of Y database):  (200, 100)
Shape of reconstructed database: (200, 100)
RMSE for whole database (pyLOM) = 9.798795481470964e-16
RMSE for whole database (numpy) = 1.413993616509635

fig2

ArnauMiro commented 5 months ago

Hola,

Lo que dices es cierto y funciona tal como lo describes. Por lo que entiendo usas la versión compilada de pyLOM. Eso es correcto, ya que da mucha más performace que la versión de python puro. Sin embargo, al usar la versión compilada hay algunas cosas que se deben tener en cuenta y hay un pequeño detalle que hace que pyLOM no funcione bien. Te cuento a continuación.

Vayamos primero a la versión sin compilar. Para acceder a ella debes instalar con

USE_COMPILED  = OFF

En el archivo options.cfg. Si has instalado con

make install_dev

simplemente puedes hacer

make cleanall

para acceder a la versión sin compilar. Porque te digo esto: pues porque verás que en la versión sin compilar tu ejemplo funciona a la perfección.

¿Entonces porque no te funciona el ejemplo? Funcionará si haces lo siguiente:

Q = Y[:, :, 0].copy() # Snapshot matrix: First X of database Y

ese .copy() es esencial. Para obtener perfomance en el compilado usamos los buffers de python directamente (y así no tener que gestionar la memoria), lo que nos ahorra algún dolor de cabeza siempre que seamos cuidadosos y entendemos que hace python por detrás. El comando

Q = Y[:, :, 0] # Snapshot matrix: First X of database Y

lo que hace es copiar el apuntador de Y a Q. Para python, el buffer de Q se define como (200,100) sin embargo sigue apuntando a Q que es de (200,100,5). Al llegar a nuestra layer de código C, el apuntador que vemos, aunque le pasemos Q es de (200,100,5), es decir, Y. Python copia apuntadores y cambia cierta metadata para no tener que copiar todo el espacio de memoria, lo que se llama una soft copy (MATLAB y Julia, por ejemplo, no hacen eso). Al poner el .copy() estamos creando el buffer de (200,100) que pyLOM necesita para trabajar y que no gestionamos internamente. Y luego tu ejemplo funciona.

Lo que no entendí de tu primer ejemplo es el porque de la tercera dimensión. Por lo que me has dicho ahora sería como el numero de variables del dataset, es decir, en tu Y tendrías (200,100,5) con 200 puntos espaiales, 100 puntos temporales y 5 variables. Entonces, para correr pyLOM en todo el dataset los tendrías que apilar de forma que te quede una base de datos tipo (200*5,100).

Espero haberme explicado bien!

Por cierto, repasando tu primer mensaje. Veo que hay problemas con el OneAPI. Nosotros también nos hemos encontrado con eso. Ahora no recuerdo si por defecto funciona o Intel ha cambiado todos los links, de todas formas le daremos una ojeada a ver.

Saludos,

Miguel-San commented 5 months ago

Hola Arnau,

Muchas gracias por tu respuesta. Ya entiendo dónde estaba el problema!

Un saludo y gracias de nuevo!