Atelier 1 - Introduction à la variance de blocs

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import VBox, Layout, IntSlider, FloatSlider, Dropdown, Output
from numpy.fft import fft2, ifft2, fftshift
from scipy.stats import norm
from IPython.display import display

# --- 1. FONCTIONS OPTIMISÉES ET CORRIGÉES ---

def theoretical_block_variance_fast(range_x, range_y, sill, nugget, block_size, pixel_size, angle_deg, model_name, n_points=20):
    """
    Approximation rapide de la variance de bloc par moyenne discrète.
    Cette version est optimisée pour éviter la création de très grands tableaux intermédiaires,
    ce qui la rend beaucoup plus rapide.
    """
    # Si le bloc est un simple pixel, la variance est la variance ponctuelle totale.
    if block_size <= 1:
        return sill + nugget

    theta = np.deg2rad(angle_deg)
    
    # Grille régulière de points pour discrétiser le bloc
    coords_1d = np.linspace(-block_size * pixel_size / 2, block_size * pixel_size / 2, n_points)
    x_pts, y_pts = np.meshgrid(coords_1d, coords_1d)
    x_pts_flat = x_pts.ravel()
    y_pts_flat = y_pts.ravel()

    # Calcul optimisé des différences de coordonnées (lags) entre toutes les paires de points
    dx = np.subtract.outer(x_pts_flat, x_pts_flat)
    dy = np.subtract.outer(y_pts_flat, y_pts_flat)

    # Rotation inverse des vecteurs de lag pour aligner avec l'anisotropie
    cos_t = np.cos(theta)
    sin_t = np.sin(theta)
    h_rot_x = dx * cos_t + dy * sin_t
    h_rot_y = -dx * sin_t + dy * cos_t
    
    # Calcul de la distance anisotropique équivalente
    h_aniso = np.sqrt((h_rot_x / range_x)**2 + (h_rot_y / range_y)**2)
    
    # Calcul de la matrice de covariance structurée (sans pépite)
    if model_name == 'Sphérique':
        cov_mat = np.where(h_aniso < 1, sill * (1 - 1.5 * h_aniso + 0.5 * h_aniso**3), 0.0)
    elif model_name == 'Exponentiel':
        cov_mat = sill * np.exp(-3 * h_aniso)
    elif model_name == 'Gaussien':
        cov_mat = sill * np.exp(-3 * h_aniso**2)
    else:
        cov_mat = np.zeros_like(h_aniso)
        
    # La formule utilisée ici approxime la variance de bloc.
    # C'est la variance ponctuelle totale moins la covariance structurée moyenne dans le bloc.
    cov_mean = np.mean(cov_mat)
    return cov_mean + nugget / ((block_size * pixel_size)**2)

def variogram_model(hx, hy, range_x, range_y, sill, nugget, angle_deg, model_name):
    """Calcule le modèle de variogramme théorique (sans l'effet de pépite)."""
    theta = np.deg2rad(angle_deg)
    hx_rot = hx * np.cos(theta) + hy * np.sin(theta)
    hy_rot = -hx * np.sin(theta) + hy * np.cos(theta)

    h = np.sqrt((hx_rot / range_x)**2 + (hy_rot / range_y)**2)

    if model_name == 'Sphérique':
        gamma = np.where(h < 1, sill * (1.5 * h - 0.5 * h**3), sill)
    elif model_name == 'Exponentiel':
        gamma = sill * (1 - np.exp(-3 * h))
    elif model_name == 'Gaussien':
        gamma = sill * (1 - np.exp(-3 * h**2))
    else:
        raise ValueError("Modèle de variogramme non reconnu")
    # L'effet de pépite est géré lors du calcul de la covariance (C(0) = c0 + c1)
    return gamma

# --- 2. AUTRES FONCTIONS (INCHANGÉES MAIS NÉCESSAIRES) ---

def normal_score_transform(field):
    """Transformation en scores normaux."""
    flat = field.flatten()
    ranks = flat.argsort().argsort()
    uniform_scores = (ranks + 1) / (len(flat) + 1)
    normal_scores = norm.ppf(uniform_scores)
    return normal_scores.reshape(field.shape)

def aggregate(field: np.ndarray, block_size: int) -> np.ndarray:
    """Agrége un champ par moyenne sur des blocs de taille block_size x block_size."""
    s = field.shape[0]
    # Assure que la taille du champ est un multiple de la taille du bloc
    trimmed_size = s - (s % block_size)
    trimmed = field[:trimmed_size, :trimmed_size]
    
    # Redimensionnement efficace pour calculer la moyenne par bloc
    reshaped = trimmed.reshape(
        trimmed_size // block_size, block_size,
        trimmed_size // block_size, block_size
    )
    return reshaped.mean(axis=(1, 3))

def variance_blocks(field, max_support):
    """Calcul de la variance expérimentale pour différentes tailles de blocs."""
    sizes = list(range(1, max_support + 1))
    variances = []

    # Calcul de la variance totale pour s = 1
    var_total = np.var(field, ddof=1)

    for s in sizes:
        if s == 1:
            # Pour s = 1, on renvoie la variance globale du champ
            variances.append(var_total)
        else:
            # Agréger le champ avec un support (taille de bloc) s
            aggregated_field = aggregate(field, s)
            
            # Calcul de la variance des valeurs agrégées
            block_variance = np.var(aggregated_field, ddof=1)
            variances.append(block_variance)

    return sizes[:len(variances)], variances


def fftma_simulation(size, range_x, range_y, sill, nugget, angle_deg, model_name, seed=None):
    """
    Simulation d'un champ gaussien par la méthode FFT-MA avec gestion de la périodicité.
    La covariance est calculée sur une grille de taille 2*size pour éliminer les effets de bord.
    """
    if seed is not None:
        np.random.seed(seed)

    size_big = 2 * size  # Taille étendue pour gérer la périodicité

    # Grille centrée autour de 0
    x = np.arange(-size_big // 2, size_big // 2)
    y = np.arange(-size_big // 2, size_big // 2)
    hx, hy = np.meshgrid(x, y, indexing='ij')

    # Variogramme structuré sans pépite
    gamma_structured = variogram_model(hx, hy, range_x, range_y, sill, 0, angle_deg, model_name)
    
    # Covariance : C(h) = c0 + c1 - γ(h)
    cov = sill - gamma_structured
    cov[size_big // 2, size_big // 2] += nugget  # Ajouter pépite à C(0)

    cov_fft = fft2(fftshift(cov))
    
    white_noise = np.random.normal(0, 1, (size_big, size_big))
    wn_fft = fft2(white_noise)
    
    field_fft = wn_fft * np.sqrt(np.abs(cov_fft))
    field_big = np.real(ifft2(field_fft))

    # Extraire la partie centrale de taille size x size
    center = size_big // 2
    field = field_big[center - size // 2:center + size // 2,
                      center - size // 2:center + size // 2]

    return field


# --- 3. INTERFACE INTERACTIVE (UTILISANT LE CODE OPTIMISÉ) ---

# Définition des widgets
slider_style = {'description_width': '160px'}
slider_layout = Layout(width='450px')
support_slider = IntSlider(min=1, max=50, step=1, value=1, description='Support (pixels)', style=slider_style, layout=slider_layout)
range_y_slider = FloatSlider(min=1, max=50, step=1, value=15, description='Portée petite ($a_p$)', style=slider_style, layout=slider_layout)
range_x_slider = FloatSlider(min=1, max=50, step=1, value=15, description='Portée grande ($a_g$)', style=slider_style, layout=slider_layout)
sill_slider = FloatSlider(min=0.1, max=10, step=0.1, value=1.0, description='Palier structure ($c_1$)', style=slider_style, layout=slider_layout)
nugget_slider = FloatSlider(min=0, max=1, step=0.01, value=0, description='Effet de pépite ($c_0$)', style=slider_style, layout=slider_layout)
angle_slider = FloatSlider(min=0, max=180, step=1, value=30, description='Angle ($\\theta$)', style=slider_style, layout=slider_layout)
model_dropdown = Dropdown(options=['Sphérique', 'Exponentiel', 'Gaussien'], value='Sphérique', description='Modèle', style=slider_style, layout=slider_layout)
button = widgets.Button(description="Lancer la simulation", button_style='success', layout=Layout(width='200px', margin='10px 0 0 165px'))
out = Output()

def interactive_variance_plot(support, range_x, range_y, sill, nugget, angle_deg, model_name):
    """Fonction principale pour générer les graphiques."""
    field = fftma_simulation(
        size=1024, range_x=range_x, range_y=range_y, sill=sill, nugget=nugget,
        angle_deg=angle_deg, model_name=model_name, seed=4263
    )
    # La variance des scores normaux est 1. On remet à l'échelle à la variance désirée.
    field = normal_score_transform(field) * np.sqrt(sill + nugget)

    plt.figure(figsize=(14, 6))

    # Agrégation du champ simulé
    agg = aggregate(field, support)

    # Définir les limites de la colorbar pour la cohérence visuelle
    vmin, vmax = -2.5 * np.sqrt(sill), 2.5 * np.sqrt(sill)
    
    plt.subplot(1, 2, 1)
    plt.imshow(agg, cmap='viridis', origin='lower', vmin=vmin, vmax=vmax)
    plt.title(f'Champ agrégé (support = {support}x{support})')
    plt.colorbar(label='Valeur')
    plt.axis('off')

    # Calcul des variances expérimentale et théorique
    max_support = 50
    bsizes, var_exp = variance_blocks(field, max_support)
    
    # Appel à la fonction RAPIDE
    var_theo = [
        theoretical_block_variance_fast(
            range_x=range_x, range_y=range_y, sill=sill, nugget=nugget, 
            block_size=b, pixel_size=1.0, # pixel_size=1.0 si on travaille en unités pixel
            angle_deg=angle_deg, model_name=model_name, n_points=100
        ) for b in bsizes
    ]

    plt.subplot(1, 2, 2)
    plt.plot(bsizes, var_exp, 'o-', label='Variance expérimentale', zorder=2)
    plt.plot(bsizes, var_theo, 's--', label='Variance théorique', zorder=1)
    plt.axvline(support, color='red', linestyle=':', label=f'Support sélectionné = {support}')
    plt.xlabel('Taille du support (pixels)')
    plt.ylabel('Variance de bloc')
    plt.title('Variance des blocs vs. Taille de support')
    plt.legend(loc='upper right')
    plt.grid(True, linestyle='--', alpha=0.6)
    plt.xlim(0, max_support)
    plt.ylim(bottom=0, top=(sill + nugget) * 1.05)
    plt.tight_layout()
    plt.show()

def on_button_clicked(b):
    """Action déclenchée par le bouton."""
    with out:
        out.clear_output(wait=True)
        interactive_variance_plot(
            support=support_slider.value, range_x=range_x_slider.value, range_y=range_y_slider.value,
            sill=sill_slider.value, nugget=nugget_slider.value, angle_deg=angle_slider.value,
            model_name=model_dropdown.value
        )

button.on_click(on_button_clicked)

# Affichage de l'interface utilisateur
ui = VBox([
    support_slider, range_x_slider, range_y_slider, sill_slider,
    nugget_slider, angle_slider, model_dropdown, button
])

print("Configurez les paramètres et cliquez sur 'Lancer la simulation'.")
display(ui, out)
print("NOTE - La discrétisation est simplifiée pour accélérer les calculs.")