Atelier 5 - Ajustement d’un modèle théorique anisotrope (modèle 2D)

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider, IntSlider, Dropdown
import matplotlib.cm as cm

# --- FONCTIONS DE BASE ---

def anisotropic_covariance(hx, hy,
                           ranges_angles_models_weights,
                           c0=0.0):
    """
    Calcule une covariance imbriquée avec plusieurs structures anisotropes.

    ranges_angles_models_weights : liste de tuples (range_major, range_minor, angle_deg, model, c)
    c0 : effet nugget (valeur à h=0 uniquement)
    """
    cov_total = np.zeros_like(hx, dtype=float)

    for range_major, range_minor, angle_deg, model, c in ranges_angles_models_weights:
        angle_rad = np.deg2rad(angle_deg)
        cos_a, sin_a = np.cos(angle_rad), np.sin(angle_rad)
        h_rot = cos_a * hx + sin_a * hy
        v_rot = -sin_a * hx + cos_a * hy
        h_scaled = h_rot / range_major
        v_scaled = v_rot / range_minor
        h = np.sqrt(h_scaled**2 + v_scaled**2)

        if model == 'spherical':
            base_cov = np.where(h < 1, 1 - 1.5 * h + 0.5 * h**3, 0)
        elif model == 'exponential':
            base_cov = np.exp(-3 * h)
        elif model == 'gaussian':
            base_cov = np.exp(-3 * h**2)
        elif model == 'nugget':
            base_cov = np.where(h == 0, 1, 0)
        else:
            raise ValueError("Modèle inconnu : {}".format(model))

        cov_total += c * base_cov

    # Nugget uniquement au point h=0
    cov_total += c0 * (hx == 0) * (hy == 0)

    return cov_total


def fftma_2d_nested(nx, ny, nested_params, c0=0.0, seed=None):
    if seed is not None:
        np.random.seed(seed)
    nfft_x, nfft_y = 2 * nx, 2 * ny
    x = np.arange(-nx, nx)
    y = np.arange(-ny, ny)
    hx, hy = np.meshgrid(x, y, indexing='ij')

    cov = anisotropic_covariance(hx, hy, nested_params, c0)

    spectrum = np.fft.fft2(np.fft.ifftshift(cov))
    spectrum = np.real(spectrum)
    spectrum[spectrum < 0] = 0
    amp = np.sqrt(spectrum)
    noise = np.random.normal(0, 1, (nfft_x, nfft_y))
    noise_fft = np.fft.fft2(noise)
    field_fft = amp * noise_fft
    field = np.fft.ifft2(field_fft).real
    field = field[:nx, :ny]
    return field


def directional_variogram(field, n_samples=500, max_lag=20, lag_step=1, directions_deg=None, tol=10):
    nx, ny = field.shape
    X, Y = np.meshgrid(np.arange(nx), np.arange(ny), indexing='ij')
    indices = np.random.choice(nx * ny, size=n_samples, replace=False)
    coords = np.column_stack((X.ravel()[indices], Y.ravel()[indices]))
    values = field.ravel()[indices]
    h_vectors = coords[:, None, :] - coords[None, :, :]
    h_dist = np.linalg.norm(h_vectors, axis=2)
    h_angle = np.arctan2(h_vectors[..., 1], h_vectors[..., 0]) * 180 / np.pi
    h_angle = (h_angle + 180) % 180
    gamma_by_dir = {}
    if directions_deg is None:
        directions_deg = np.arange(0, 180, 22.5)
    lag_bins = np.arange(0, max_lag + lag_step, lag_step)
    lag_centers = 0.5 * (lag_bins[:-1] + lag_bins[1:])
    for theta in directions_deg:
        mask_dir = (np.abs(h_angle - theta) <= tol)
        gammas = []
        for k in range(len(lag_bins) - 1):
            mask_lag = (h_dist >= lag_bins[k]) & (h_dist < lag_bins[k + 1])
            mask = mask_dir & mask_lag
            if not np.any(mask):
                gammas.append(np.nan)
                continue
            i, j = np.where(mask)
            gamma = 0.5 * np.mean((values[i] - values[j])**2)
            gammas.append(gamma)
        gamma_by_dir[theta] = gammas
    return gamma_by_dir, lag_centers, coords, values, X, Y, field

# --- SIMULATION FIXE ---

nx, ny = 250, 250
seed = 42

# Deux structures imbriquées + nugget
nested_params = [
    (80, 20, 22.5, 'spherical', 0.7),   # Structure large orientée N-S
    (1, 1, 1, 'nugget', 0.3)     # Structure fine orientée E-W
]

field = fftma_2d_nested(nx, ny, nested_params, c0=0.1, seed=42)

# --- WIDGET D'AJUSTEMENT ---
def update_model(n_samples, lag_step, c0, c1, model, range_major, range_minor, angle_deg):

    gamma_exp, lags, coords, values, X, Y, _ = directional_variogram(field, n_samples=n_samples, max_lag=50, lag_step=lag_step, tol=10)

    directions = list(gamma_exp.keys())
    fig, axs = plt.subplots(3, 3, figsize=(10, 8))
    for i in range(8):
        angle = directions[i]
        ax = axs[i // 3, i % 3]
        ax.plot(lags, gamma_exp[angle], 'o-', label='Expérimental', color='tab:blue')
        # Courbe théorique
        h = lags / range_major  # Approximation isotrope sur direction principale
        if model == 'spherical':
            gamma_th = np.where(h < 1, c0 + c1 * (1.5 * h - 0.5 * h**3), c0 + c1)
        elif model == 'exponential':
            gamma_th = c0 + c1 * (1 - np.exp(-3 * h))
        elif model == 'gaussian':
            gamma_th = c0 + c1 * (1 - np.exp(-3 * h**2))
        ax.plot(lags, gamma_th, '-', color='black', label='Théorique')
        ax.set_title(f'Direction {angle}°')
        ax.set_xlabel('Lag')
        ax.set_ylabel('γ(h)')
        ax.set_ylim(0, 1.8)
        ax.set_xlim(0, np.max(lags))
        ax.grid(True)
    axs[2, 2].axis('off')
    plt.tight_layout()
    plt.show()
    
    # Figure 2 : champ simulé + points échantillonnés
    fig2, ax = plt.subplots(figsize=(6, 4))
    im = ax.imshow(field, cmap='jet', origin='lower', alpha=0.4)
    ax.set_title('Champ simulé avec points échantillonnés')

    nx_grid, ny_grid = field.shape
    X_grid, Y_grid = np.meshgrid(np.arange(nx_grid), np.arange(ny_grid), indexing='ij')

    sampled_indices = coords[:,0]*ny_grid + coords[:,1]
    
    mask_all = np.ones(nx_grid * ny_grid, dtype=bool)
    mask_all[sampled_indices] = False

    X_all = X_grid.ravel()[mask_all]
    Y_all = Y_grid.ravel()[mask_all]
    vals_all = field.ravel()[mask_all]

    X_samp = coords[:,0]
    Y_samp = coords[:,1]
    vals_samp = values

    norm = plt.Normalize(vmin=np.min(field), vmax=np.max(field))
    cmap = cm.jet
    colors = cmap(norm(vals_samp))
    ax.scatter(Y_samp, X_samp, s=80, facecolors='none', edgecolors='black', linewidth=1.2, label='Échantillonnés')
    sc = ax.scatter(Y_samp, X_samp, s=60, c=vals_samp, cmap='jet', norm=norm, label='Valeur')

    plt.colorbar(sc, ax=ax, fraction=0.046, pad=0.04)
    ax.set_xlim(0, nx-1)
    ax.set_ylim(0, ny-1)

    plt.tight_layout()
    plt.show()

# Créer le widget interactif
interact(update_model,
         range_major=FloatSlider(min=5, max=50, step=1, value=20, description='Range major'),
         range_minor=FloatSlider(min=2, max=50, step=1, value=10, description='Range minor'),
         angle_deg=IntSlider(min=0, max=180, step=5, value=45, description='Angle (°)'),
         model=Dropdown(options=['spherical', 'exponential', 'gaussian'], value='spherical', description='Model'),
         c0=FloatSlider(min=0, max=1, step=0.01, value=0, description='Nugget c0'),
         c1=FloatSlider(min=0, max=2, step=0.01, value=1, description='Sill c1'),
         n_samples=IntSlider(min=50, max=2000, step=25, value=100, description='n_samples'),
         lag_step=IntSlider(min=1, max=20, step=1, value=5, description='Lag'))