Atelier 2 - Effet d’information

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft2, ifft2, fftshift
from numpy.random import Generator, PCG64
import ipywidgets as widgets
from ipywidgets import interact, Layout
from functools import partial

# --- Configuration Centralisée ---
# Regrouper les paramètres facilite la lecture et les modifications.
CONFIG = {
    "size": 200,
    "range": 0.2,
    "sill": 1.0,
    "seed": 42,
    "dtype": np.float32,  # Utiliser float32 réduit l'usage mémoire de 50%
    "clip_min": 0,
    "clip_max": 10,
}

# --- Fonctions de Simulation (Améliorées) ---
# Ces fonctions sont maintenant plus lisibles et utilisent les pratiques modernes.

def spherical_covariance_fft(size: int, range_val: float, sill: float = 1.0) -> np.ndarray:
    """Calcule la covariance sphérique sur une grille étendue pour la simulation FFT."""
    extended_size = 2 * size
    x = np.arange(-extended_size // 2, extended_size // 2)
    X, Y = np.meshgrid(x, x)
    h = np.sqrt(X**2 + Y**2).astype(CONFIG["dtype"])
    
    # Évite la division par zéro si range_val est nul
    effective_range = range_val * size
    if effective_range == 0:
        return np.zeros((extended_size, extended_size), dtype=CONFIG["dtype"])

    h_norm = h / effective_range
    cov = sill * (1 - 1.5 * h_norm + 0.5 * h_norm**3)
    cov[h > effective_range] = 0
    return fftshift(cov)

def fftma_simulation(size: int, range_val: float, rng: Generator, sill: float = 1.0) -> np.ndarray:
    """Génère un champ gaussien 2D en utilisant la méthode FFT-MA."""
    extended_size = 2 * size
    cov_model = spherical_covariance_fft(size, range_val, sill)
    cov_fft = np.sqrt(np.abs(fft2(cov_model)))
    
    white_noise = rng.normal(size=(extended_size, extended_size)).astype(CONFIG["dtype"])
    white_fft = fft2(white_noise)
    
    z_fft = cov_fft * white_fft
    z_ext = np.real(ifft2(z_fft))
    
    start, end = extended_size // 4, extended_size // 4 + size
    return z_ext[start:end, start:end]

def gaussian_to_lognormal(field: np.ndarray) -> np.ndarray:
    """Convertit un champ gaussien en champ log-normal."""
    return np.exp(field)

# --- ÉTAPE 1: Calcul unique et coûteux ---
# On génère le champ de base une seule fois ici.
rng = np.random.default_rng(CONFIG["seed"])
base_field_gauss = fftma_simulation(CONFIG["size"], CONFIG["range"], rng, CONFIG["sill"])
REAL_FIELD = gaussian_to_lognormal(base_field_gauss)



# --- ÉTAPE 2: Fonction d'affichage rapide ---
# Cette fonction ne fait que des calculs rapides et prend le champ réel en argument.
def plot_real_vs_estimated_bias(
    real_field: np.ndarray,
    bias_percent: float,
    noise_std: float,
    cutoff: float
):
    """Affiche les champs et le nuage de points en fonction des paramètres interactifs."""
    
    # Ajout du biais et du bruit (opérations rapides)
    biased_real = real_field * (1 + bias_percent / 100)
    # On utilise un nouveau générateur pour que le bruit change à chaque interaction
    noise_rng = np.random.default_rng()
    epsilon = noise_rng.normal(loc=0, scale=noise_std, size=real_field.shape).astype(CONFIG["dtype"])
    estimated = biased_real + epsilon
    
    # Clipping des valeurs pour l'affichage
    vmin, vmax = CONFIG["clip_min"], CONFIG["clip_max"]
    real_clipped = np.clip(real_field, vmin, vmax)
    estimated_clipped = np.clip(estimated, vmin, vmax)
    
    # --- Tracé ---
    fig, axes = plt.subplots(2, 2, figsize=(10, 8), constrained_layout=True)
    axes = axes.ravel()  # Transforme en tableau 1D de 4 éléments
    
    # Cartes des champs réel et estimé
    im0 = axes[0].imshow(real_clipped, cmap='viridis', vmin=vmin, vmax=vmax)
    axes[0].set_title("Champ réel")
    fig.colorbar(im0, ax=axes[0])
    
    im1 = axes[1].imshow(estimated_clipped, cmap='viridis', vmin=vmin, vmax=vmax)
    axes[1].set_title("Champ estimé")
    fig.colorbar(im1, ax=axes[1])

    for ax in [axes[0], axes[1]]:
        ax.set_xlabel("X (m)")
        ax.set_ylabel("Y (m)")

    # Nuage de points
    ax2 = axes[2]
    real_flat = real_clipped.flatten()
    estimated_flat = estimated_clipped.flatten()

    # Masques pour colorer les points
    is_ore_real = real_flat >= cutoff
    is_ore_estimated = estimated_flat >= cutoff
    mask_red = is_ore_real & ~is_ore_estimated   # Minerai ignoré
    mask_blue = ~is_ore_real & is_ore_estimated # Stérile traité
    mask_other = ~(mask_red | mask_blue)

    # --- Calcul des pourcentages ---
    total_points = len(real_flat)
    percent_red = 100 * np.sum(mask_red) / total_points
    percent_blue = 100 * np.sum(mask_blue) / total_points

    ax2.scatter(estimated_flat[mask_other], real_flat[mask_other], alpha=0.3, s=10, color="gray", label="Correctement classé")
    ax2.scatter(estimated_flat[mask_blue], real_flat[mask_blue], alpha=0.7, s=15, color="blue", ec="k", lw=0.2,
                label=f"Stérile traité ({percent_blue:.1f}%)")
    ax2.scatter(estimated_flat[mask_red], real_flat[mask_red], alpha=0.7, s=15, color="red", ec="k", lw=0.2,
                label=f"Minerai ignoré ({percent_red:.1f}%)")

    
    ax2.plot([vmin, vmax], [vmin, vmax], 'k-', lw=2, label="Ligne 1:1")
    ax2.axhline(cutoff, color='k', linestyle='--', lw=1)
    ax2.axvline(cutoff, color='k', linestyle='--', lw=1)

    ax2.set(xlabel="Teneur estimée (ppm)", ylabel="Teneur réelle (ppm)",
            title="Teneur réelle vs estimée",
            xlim=(vmin, vmax), ylim=(vmin, vmax))
    ax2.set_aspect('equal', adjustable='box')
    ax2.legend(loc='upper right')
    ax2.grid(True, linestyle=':', alpha=0.5)

    # --- Ligne de régression ---
    coeffs = np.polyfit(estimated_flat, real_flat, deg=1)
    reg_x = np.linspace(vmin, vmax, 100)
    reg_y = coeffs[0] * reg_x + coeffs[1]
    ax2.plot(reg_x, reg_y, 'g--', lw=2, label=f"Régression linéaire: y = {coeffs[0]:.2f}x + {coeffs[1]:.2f}")

    axes[3].axis("off")  # Masque la case vide (position bas droite)

    plt.show()

# --- ÉTAPE 3: Widgets et Interaction ---
# On crée les widgets comme avant.

# Widgets
common_layout = Layout(width='400px')  # Largeur totale du widget

bias_slider = widgets.FloatSlider(
    value=0, min=-50, max=50, step=5,
    description='Biais (%)',
    layout=common_layout,
    style={'description_width': '150px'}  # <- largeur de l'étiquette
)

noise_slider = widgets.FloatSlider(
    value=0.5, min=0, max=2.0, step=0.05,
    description='Écart-type bruit',
    layout=common_layout,
    style={'description_width': '150px'}
)

cutoff_slider = widgets.FloatSlider(
    value=2, min=1, max=8, step=0.1,
    description='Teneur coupure',
    layout=common_layout,
    style={'description_width': '150px'}
)
# Au lieu d'utiliser partial, on utilise une fonction lambda.
# La lambda accepte les arguments des curseurs et appelle notre fonction
# principale en lui passant ces arguments ainsi que le champ pré-calculé.
interact(
    lambda bias_percent, noise_std, cutoff: plot_real_vs_estimated_bias(REAL_FIELD, bias_percent, noise_std, cutoff),
    bias_percent=bias_slider,
    noise_std=noise_slider,
    cutoff=cutoff_slider
);