Atelier 1 - Définition variogramme et covariance spatiale

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, Dropdown, BoundedFloatText, BoundedIntText

# --- Modèles de covariance ---
def spherical_model(h, range_, sill):
    h = np.abs(h)
    return np.where(h <= range_, sill * (1 - 1.5 * h / range_ + 0.5 * (h / range_)**3), 0)

def exponential_model(h, range_, sill):
    h = np.abs(h)
    return sill * np.exp(-3 * h / range_)

def gaussian_model(h, range_, sill):
    h = np.abs(h)
    return sill * np.exp(-np.sqrt(3) * (h / range_)**2)

def nugget_model(h, sill):
    return sill * (h == 0)

def compute_covariance(model, range_, sill, n):
    h = np.fft.fftshift(np.arange(-n//2, n//2))
    if model == "Sphérique":
        cov = spherical_model(h, range_, sill)
    elif model == "Exponentiel":
        cov = exponential_model(h, range_, sill)
    elif model == "Gaussien":
        cov = gaussian_model(h, range_, sill)
    elif model == "Effet de pépite":
        cov = nugget_model(h, sill)
    return cov

# --- FFT-MA 1D ---
def fftma_1d(model="Sphérique", range_=20.0, sill=1.0, n=1024, seed=1542):
    if seed is not None:
        np.random.seed(seed)
    
    cov = compute_covariance(model, range_, sill, 2*n)
    spectrum = np.fft.fft(cov)
    spectrum = np.where(spectrum.real < 0, 0, spectrum)
    amp = np.sqrt(spectrum)
    noise = np.random.normal(0, 1, 2*n)
    field_fft = amp * np.fft.fft(noise)
    field = np.fft.ifft(field_fft).real[:n]
    field = (field - field.mean()) / field.std()
    field *= np.sqrt(sill)
    return field

# --- Affichage interactif avec deux figures ---
def plot_variogram_with_stats(model, range_, sill, h_shift):
    n = 10000
    z = fftma_1d(model, range_, sill, n)
    
    h = int(h_shift)
    z_x = z[:-h]
    z_xh = z[h:]

    cov_emp = np.cov(z_x, z_xh)[0, 1]
    corr_emp = np.corrcoef(z_x, z_xh)[0, 1]
    gamma_emp = 0.5 * np.mean((z_x - z_xh)**2)
    length = 2 * gamma_emp

    h_vals = np.arange(1, 100)
    gamma_vals = []
    cov_vals = []
    corr_vals = []
    for hh in h_vals:
        zx = z[:-hh]
        zxh = z[hh:]
        gamma_vals.append(0.5 * np.mean((zx - zxh)**2))
        cov_vals.append(np.cov(zx, zxh)[0, 1])
        corr_vals.append(np.corrcoef(zx, zxh)[0, 1])

    fig, axs = plt.subplots(1, 2, figsize=(14, 6))

    axs[0].scatter(z_x, z_xh, alpha=0.4)
    axs[0].plot([-4, 4], [-4, 4], 'k--', lw=1, label="z(x) = z(x+h)")
    mean_z = 0
    dx, dy = 1 / np.sqrt(2), -1 / np.sqrt(2)
    x0, y0 = mean_z - dx * length, mean_z - dy * length
    x1, y1 = mean_z + dx * length, mean_z + dy * length
    axs[0].plot([x0, x1], [y0, y1], 'r-', lw=4, label=r"$2\gamma(h)$")
    axs[0].set_xlabel("z(x)")
    axs[0].set_ylabel("z(x + h)")
    axs[0].set_title(f"{model} | h = {h} px\ncov = {cov_emp:.2f}, ρ = {corr_emp:.2f}, γ(h) = {gamma_emp:.2f}")
    axs[0].legend(loc='upper left')
    axs[0].grid(True)
    axs[0].set_xlim(-4, 4)
    axs[0].set_ylim(-4, 4)
    axs[0].set_aspect('equal', adjustable='box')

    axs[1].plot(h_vals, gamma_vals, 'r*', label=r"$\gamma(h)$ - Variogramme")
    axs[1].plot(h_vals, cov_vals, 'b*', label="Covariance")
    axs[1].plot(h_vals, corr_vals, 'g*', label="Corrélation")
    axs[1].axvline(h, color='k', linestyle='--', lw=1)
    axs[1].set_xlabel("Distance h (pixels)")
    axs[1].set_title("Évolution du variogramme, covariance et corrélation")
    axs[1].legend(loc='upper right')
    axs[1].grid(True)
    axs[1].set_xlim(0, 100)
    axs[1].set_ylim(-0.2*sill, max(1.5*sill, 1))

    plt.tight_layout()
    plt.show()

# --- Interface interactive avec champs bornés ---
interact(plot_variogram_with_stats,
         model=Dropdown(options=["Sphérique", "Exponentiel", "Gaussien", "Effet de pépite"], value="Sphérique", description="Modèle"),
         range_=BoundedFloatText(value=20.0, min=5.0, max=100.0, step=1.0, description="Portée"),
         sill=BoundedFloatText(value=1.0, min=0.1, max=2.0, step=0.1, description="Variance"),
         h_shift=BoundedIntText(value=5, min=1, max=100, step=1, description="Décalage h"));