fbientrigo / grav_lensing

A Python toy model for gravitational lensing with educative means
0 stars 0 forks source link

Separar el modelo en 2 estimadores #6

Closed fbientrigo closed 1 month ago

fbientrigo commented 1 month ago

Se tienen dos metodos bastante prometedores, una es el uso de Principal Component Analysis (PCA) y la otra el uso de Gaussian Mixtures (GM).

PCA

El PCA es bastante bueno cuando tenemos muchos datos y poca variabilidad en ellos, permite encontrar los patrones y generar bases que expliquen de mejor manera los datos observados. Para funcionar bien requiere de datos escalados correctamente. Un modelo que utilice PCA sería el que generé los coeficientes de cada componente para construir la imagen final.

requiere

Un pipeline

Un modelo que tenga

Gaussian Mixtures Model

Utilizando la distribución como una distribución probabilistica (necesario el escalado) se pueden generar un muestro de puntos, de manera (x,y), esta colección de puntos se puede pasar a un modelo de GM para que estime donde posicionar las gaussianas para mejor adaptarse a los datos, generando los coeficientes del modelo, las medias y las desviaciones estandar $GMM(c_\pi, \vec \mu, \vec \sigma) = \sumj c\pi(j) \exp(- (\vec r - \vec \mu(j))^2/2\sigma(j)^2)$ Un modelo que estima GM sería entonces una red que entregue 5 parametros por cada Gaussiana.

requiere

Un pipeline

Objetivo del Issue

Generar estudios de ambos modelos, y estimar como aprovecharlos para crear eventualmente una red que sea capaz de usar uno o ambos.

Posibles problematicas

Un enfoque naive sería el de directamente darle las bases del PCA, que la red estime los coeficientes y luego además darle la posición de los peaks de la misma imagen, provocaría que el modelo intente generar la misma información dos veces, consiguiendo una correlación indeseada.

Posibles arreglos: Separar la información

Si se puede separar la imagen en 2, una que sea más facil de explicar con PCA y una que sea más facil con GMM, entonces el modelo puede aprender a generar las dos partes por separado y así conseguir un mejor resultado

fbientrigo commented 1 month ago

Avances con la separación de información

Usando la transformada de Fourier, se implementaron hard-filters de manera

def apply_lowpass_filter(f_transform, cutoff):
    """
    Aplica un filtro pasa bajas (low-pass) a la Transformada de Fourier.

    Parámetros:
        f_transform (numpy array): Transformada de Fourier desplazada.
        cutoff (float): Frecuencia de corte para el filtro (entre 0 y 1, relativa al tamaño de la imagen).

    Retorna:
        lowpass_filtered (numpy array): Transformada de Fourier filtrada con pasa bajas.
    """
    rows, cols = f_transform.shape
    crow, ccol = rows // 2, cols // 2  # Centro de la imagen

    # Crear una máscara pasa bajas (con 1 en el centro y 0 en los bordes)
    mask = np.zeros((rows, cols), dtype=np.float32)
    radius = cutoff * min(crow, ccol)  # Radio de corte basado en la frecuencia relativa
    for i in range(rows):
        for j in range(cols):
            if np.sqrt((i - crow) ** 2 + (j - ccol) ** 2) <= radius:
                mask[i, j] = 1

    # Aplicar la máscara a la Transformada de Fourier
    lowpass_filtered = f_transform * mask
    return lowpass_filtered

def apply_highpass_filter(f_transform, cutoff):
    """
    Aplica un filtro pasa altas (high-pass) a la Transformada de Fourier.

    Parámetros:
        f_transform (numpy array): Transformada de Fourier desplazada.
        cutoff (float): Frecuencia de corte para el filtro (entre 0 y 1, relativa al tamaño de la imagen).

    Retorna:
        highpass_filtered (numpy array): Transformada de Fourier filtrada con pasa altas.
    """
    rows, cols = f_transform.shape
    crow, ccol = rows // 2, cols // 2  # Centro de la imagen

    # Crear una máscara pasa altas (con 0 en el centro y 1 en los bordes)
    mask = np.ones((rows, cols), dtype=np.float32)
    radius = cutoff * min(crow, ccol)  # Radio de corte basado en la frecuencia relativa
    for i in range(rows):
        for j in range(cols):
            if np.sqrt((i - crow) ** 2 + (j - ccol) ** 2) <= radius:
                mask[i, j] = 0

    # Aplicar la máscara a la Transformada de Fourier
    highpass_filtered = f_transform * mask
    return highpass_filtered

De manera que fuera posible separar imagenes en 2, una que pase el filtro y nada pase al otro y viceversa. Esto provoca que podamos sumar los dos resultados en la imagen original y asi tener información separable.

image

Posibles errores

Si bien existen correlaciones entre la densidad, una zona con un peak concentrado (alta frecuencia) tambien presentará una mayor densidad a gran escala, lo que podría provocar que la información no sea totalmente ortogonal e independiente. Por ello se requiere comprobar este approach con experimentación

fbientrigo commented 1 month ago

Aplicando PCA

Primero se implementó el IncrementalPCA de sklearn, el cual permite cargar la información de a poco y así no tener que alocar los gigas de memoria.

A lo cual con unos 32 componentes se explica apenas un 25% de la varianza y se requieren más de 100 componentes para explicar el 70% de la varianza image

De este resultado se pueden ver interesantes bases aprendidas image

El modelo es capaz de construir imagenes aproximadas asumiendo el uso de los 256 componentes que se entrenarón con todo el dataset image estas imagenes generadas asumen que la red aprendería a la perfección los componentes, por lo que estamos observando la información maxima contenida en este metodo.

Usar ambos componentes por separado

Sobre ambos componentes, los de baja y alta frecuencia se aplicó el PCA, lo cual tomo alrededor de 90 minutos.

# Crear y aplicar IncrementalPCA para baja frecuencia
ipca_low = IncrementalPCA(n_components=n_components, batch_size=batch_size)

# Crear y aplicar IncrementalPCA para alta frecuencia
ipca_high = IncrementalPCA(n_components=n_components, batch_size=batch_size)

# Extraer imágenes del dataset en lotes y ajustar el IncrementalPCA
for X_batch, y_batch in train_dataset:
    # Obtener las imágenes de baja y alta frecuencia
    low_freq_batch, high_freq_batch = process_batch_filters(y_batch.numpy(), cutoff=0.05)

    # Aplanar las imágenes y apilarlas para baja frecuencia
    low_freq_stack = np.vstack([img.reshape(-1, 128*128) for img in low_freq_batch])

    # Aplanar las imágenes y apilarlas para alta frecuencia
    high_freq_stack = np.vstack([img.reshape(-1, 128*128) for img in high_freq_batch])

    # Aplicar el ajuste incremental al IPCA para baja frecuencia
    ipca_low.partial_fit(low_freq_stack)

    # Aplicar el ajuste incremental al IPCA para alta frecuencia
    ipca_high.partial_fit(high_freq_stack)

A lo cual se consiguierón los resultados en el siguiente grafico image

con unos 32 componentes que fue el maximo utilizado en este caso, el modelo puede explicar un 68% de la varianza para la baja frecuencia, loq ue sugiere que podemos usar PCA para estimar el fondo poco variable. Y los componentes de alta frecuencia que son gaussianas más concentradas estimarlas con el modelo de gaussian mixtures

fbientrigo commented 1 month ago

Los componentes de baja frecuencia image

Y de alta frecuencia image

Es posible generar además el vector de coeficientes de forma facil, primero se tiene que cargar el modelo que esta guardado en 05_parametricNet/ipca_low.pkl

# Cargar el modelo PCA para baja frecuencia
with open('ipca_low.pkl', 'rb') as f:
    ipca_low_loaded = pickle.load(f)

# generar los vectores
coefficients = ipca_low_loaded.transform(image.reshape(1,-1))

# bases de imagenes
eigen_images_low = ipca_low_loaded.components_[:n_components].reshape((n_components, 128, 128))

En 2_freq_separation.ipynb se desarrolla con batchs, a lo que vemos la capacidad de reconstrucción para el background

image

fbientrigo commented 1 month ago

Se ha logrado implementar el modelo de dos salidas los cambios se pueden encontrar en 06_model_parametric en particular "a_parametrize_tutorial.ipynb" para una introducción de como trabajar con el modelo y sus salidas