Atelier 4 - Impact du nombre de données sur l’estimation et l’ajustement du variogramme

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

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

def exponential_cov(h, range_):
    return np.exp(-3*h / range_)

def gaussian_cov(h, range_):
    return np.exp(-np.sqrt(3)*(h / range_)**2)

def nugget_model(h):
    return np.where(h == 0, 1, 0)

def get_cov_model(name, h, range_, sill):
    if name == "Spherique":
        return spherical_cov(h, range_) * sill
    elif name == "Exponentiel":
        return exponential_cov(h, range_) * sill
    elif name == "Gaussien":
        return gaussian_cov(h, range_) * sill
    elif name == "Pepite":
        return nugget_model(h) * sill
    return np.zeros_like(h)

# --- FFT-MA 1D ---
def compute_nested_covariance(h, models):
    cov = np.zeros_like(h, dtype=float)
    for model in models:
        type_ = model["type"]
        range_ = model["range"]
        sill = model.get("sill", 1.0)

        if type_ == "Spherique":
            cov += sill * spherical_cov(h, range_)
        elif type_ == "Exponentiel":
            cov += sill * exponential_cov(h, range_)
        elif type_ == "Gaussien":
            cov += sill * gaussian_cov(h, range_)
        elif type_ == "Pepite":
            cov += sill * nugget_model(h)
        else:
            raise ValueError(f"Modèle de covariance inconnu : {type_}")
    return cov

def fftma_1d(n, models, sigma=1.0, seed=544):
    if seed is not None:
        np.random.seed(seed)

    nfft = 2 * n  
    h = np.arange(nfft)
    cov = compute_nested_covariance(h, models)
    cov = np.concatenate([cov[:n], cov[:n][::-1]])

    spectrum = np.real(np.fft.fft(cov))
    spectrum[spectrum < 0] = 0
    amp = np.sqrt(np.maximum(spectrum, 1e-10))

    noise = np.random.normal(0, 1, nfft)
    noise_fft = np.fft.fft(noise)

    field_fft = amp * noise_fft
    field = np.fft.ifft(field_fft).real[:n]
    
    return field

# --- Variogramme expérimental ---
def empirical_variogram(z, max_lag, lag_step, x=None, n_sample=None, seed=5652):
    if x is None:
        x = np.arange(len(z))
    
    if seed is not None:
        np.random.seed(seed)

    if n_sample is not None and n_sample < len(z):
        indices_sample = np.random.choice(len(z), n_sample, replace=False)
        x = x[indices_sample]
        z = z[indices_sample]
    else:
        indices_sample = np.arange(len(z))

    n = len(z)
    distances = []
    semivars = []

    for i in range(n):
        for j in range(i + 1, n):
            h = abs(x[j] - x[i])
            sv = 0.5 * (z[j] - z[i])**2
            distances.append(h)
            semivars.append(sv)

    distances = np.array(distances)
    semivars = np.array(semivars)

    bins = np.arange(0, max_lag + lag_step * 1.5, lag_step)
    bin_indices = np.digitize(distances, bins, right=True) - 1

    gamma = []
    lags = []

    for i in range(len(bins) - 1):
        mask = bin_indices == i
        gamma.append(np.mean(semivars[mask]) if np.any(mask) else np.nan)
        lags.append((bins[i] + bins[i+1]) / 2)

    return np.array(lags), np.array(gamma), indices_sample

# --- Données simulées ---
n = 500
models_sim = [
    {"type": "Spherique", "range": 100, "sill": 0.8},
    {"type": "Pepite", "range": 0, "sill": 0.2}
]

models_adj = [
    {"type": "Spherique", "range": 161, "sill": 0.69},
    {"type": "Pepite", "range": 2, "sill": 0.20}
]

z_true = fftma_1d(n, models_sim)
lags_exp, gamma_exp, _ = empirical_variogram(z_true, max_lag=n//2, lag_step=5)

# --- Visualisation interactive ---
def plot_fitting(lag_step=5, n_sample=100):
    n_sample = int(n_sample)
    lags_exp, gamma_exp, sample_indices = empirical_variogram(z_true, max_lag=n//2, lag_step=lag_step, n_sample=n_sample)

    models_user = [
        {"type": "Spherique", "range": 161, "sill": 0.69},
        {"type": "Pepite", "range": 0, "sill": 0.20}
    ]
    gamma_model = [0.69 + 0.20 - compute_nested_covariance(np.array([lag]), models_user)[0] for lag in lags_exp]

    gamma_target = []
    sill = sum(model['sill'] for model in models_adj)
    gamma_target = [sill - compute_nested_covariance(np.array([lag]), models_adj)[0] for lag in lags_exp]

    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    axes[0].plot(lags_exp, gamma_exp, 'o', label="Variogramme expérimental")
    axes[0].plot(lags_exp, gamma_model, 'k-', label="Modèle ajusté")
    axes[0].set_xlabel("h")
    axes[0].set_ylabel("γ(h)")
    axes[0].set_title("Ajustement du variogramme")
    axes[0].grid(True)
    axes[0].legend()
    axes[0].set_xlim(0, n//2)
    axes[0].set_ylim(0, 2)

    x_full = np.arange(len(z_true))
    axes[1].plot(x_full, z_true, label="Champ simulé")
    axes[1].plot(sample_indices, z_true[sample_indices], 'ro', label="Échantillons")
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("z(x)")
    axes[1].set_title("Champ simulé et points utilisés")
    axes[1].grid(True)
    axes[1].legend()

    plt.tight_layout()
    plt.show()

interact(
    plot_fitting,
    lag_step=FloatSlider(min=1, max=50, step=1, value=5, description="Pas (lag)"),
    n_sample=FloatSlider(min=20, max=n, step=10, value=100, description="n échantillons"),
)