🛠️ Atelier 3 – Ajustement du variogramme (1D)

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider, Dropdown, Checkbox, Button, VBox, HBox, Output
from IPython.display import display, clear_output
import random

# -----------------------------
# 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, range_=None):
    return np.where(h == 0, 1, 0)

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)

        fn = {
            "Spherique": spherical_cov,
            "Exponentiel": exponential_cov,
            "Gaussien": gaussian_cov,
            "Pepite": nugget_model
        }.get(type_)

        if fn is not None:
            cov += sill * fn(h, range_)
    return cov

# -----------------------------
# Simulation 1D par FFT-MA
# -----------------------------
def fftma_1d(n, models, seed=1234):
    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)
    field_fft = amp * np.fft.fft(noise)
    field = np.fft.ifft(field_fft).real[:n]
    
    # Normalisation de la variance du champ à la somme des sill
    target_variance = sum(model.get("sill", 1.0) for model in models)
    current_variance = np.var(field)
    if current_variance > 0:
        field = field / np.sqrt(current_variance) * np.sqrt(target_variance)

    return field

# -----------------------------
# Variogramme expérimental
# -----------------------------
def empirical_variogram(z, max_lag, lag_step, x=None, n_sample=None, seed=4512):
    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))

    distances = []
    semivars = []

    n = len(z)
    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

# -----------------------------
# Interface interactive complète
# -----------------------------
n = 5000
output = Output()
solution_output = Output()

current_field = None
current_target = None

type_struct = Dropdown(options=["Spherique", "Exponentiel", "Gaussien"], value="Spherique", description="Type")
var_struct = FloatSlider(min=0.0, max=1.5, step=0.01, value=0.9, description="Variance ($c_1$)")
range_struct = FloatSlider(min=1, max=200, step=1, value=40, description="Portée (a)")
var_nugget = FloatSlider(min=0.0, max=0.6, step=0.01, value=0.1, description="Pépite ($c_0$)")
show_target = Checkbox(value=False, description="Afficher modèle cible")
btn_random = Button(description="🎲 Changer de scénario", button_style='info')
btn_show_solution = Button(description="Afficher solution", button_style='success')

def plot_fitting(var_nugget_val, type_struct_val, var_struct_val, range_struct_val, show_target_val):
    lag_step = 4
    max_lag = 250
    n_sample = 1000
    lags_exp, gamma_exp, sample_indices = empirical_variogram(current_field, max_lag=n//2, lag_step=lag_step, n_sample=n_sample)

    models_user = [
        {"type": type_struct_val, "range": range_struct_val, "sill": var_struct_val},
        {"type": "Pepite", "range": 0, "sill": var_nugget_val}
    ]

    gamma_model = [var_struct_val + var_nugget_val - compute_nested_covariance(np.array([lag]), models_user)[0] for lag in lags_exp]

    gamma_target = []
    if show_target_val:
        gamma_target = [current_target["sill"] + current_target["nugget"] - compute_nested_covariance(np.array([lag]), [
            {"type": current_target["type"], "range": current_target["range"], "sill": current_target["sill"]},
            {"type": "Pepite", "range": 0, "sill": current_target["nugget"]}
        ])[0] for lag in lags_exp]

    fig, axes = plt.subplots(2, 1, figsize=(8, 6))

    axes[0].plot(lags_exp, gamma_exp, 'o', label="Variogramme expérimental")
    axes[0].plot(lags_exp, gamma_model, 'k-', label="Modèle ajusté")
    if show_target_val:
        axes[0].plot(lags_exp, gamma_target, 'r--', label="Modèle cible", linewidth=2)
    axes[0].set_xlabel("h")
    axes[0].set_ylabel("γ(h)")
    axes[0].legend()
    axes[0].set_xlim(0, max_lag)
    axes[0].set_ylim(0, 2.5)
    axes[0].grid(True)
    axes[0].set_title("Ajustement du variogramme")

    axes[1].plot(np.arange(n), current_field, label="Champ simulé")
    axes[1].plot(sample_indices, current_field[sample_indices], 'ro', label="Échantillons")
    axes[1].legend()
    axes[1].set_title("Champ simulé et points utilisés")
    axes[1].grid(True)
    axes[1].set_ylim(-3, 3)
    axes[1].set_xlim(0, 500)

    plt.tight_layout()
    plt.show()

def generate_random_target():
    return {
        "type": random.choice(["Spherique", "Exponentiel", "Gaussien"]),
        "range": random.randint(25, 150),
        "sill": round(random.uniform(0.5, 1.0), 2),
        "nugget": round(random.uniform(0.0, 0.3), 2)
    }
    
def update_scenario():
    global current_field, current_target
    current_target = generate_random_target()
    models_sim = [
        {"type": current_target["type"], "range": current_target["range"], "sill": current_target["sill"]},
        {"type": "Pepite", "range": 0, "sill": current_target["nugget"]}
    ]
    current_field = fftma_1d(n, models_sim, seed=np.random.randint(0, 1_000_000))


def on_random_clicked(b):
    update_scenario()
    show_target.value = False
    with output:
        clear_output(wait=True)
        plot_fitting(var_nugget.value, type_struct.value, var_struct.value, range_struct.value, False)
    with solution_output:
        clear_output(wait=True)

def on_show_solution_clicked(b):
    show_target.value = True
    with output:
        clear_output(wait=True)
        plot_fitting(var_nugget.value, type_struct.value, var_struct.value, range_struct.value, True)
    with solution_output:
        clear_output(wait=True)
        print(f"🎯 Modèle cible : {current_target['type']}, a = {current_target['range']} m, c₁ = {current_target['sill']}, c₀ = {current_target['nugget']}")

btn_random.on_click(on_random_clicked)
btn_show_solution.on_click(on_show_solution_clicked)

def update_plot(change=None):
    with output:
        clear_output(wait=True)
        plot_fitting(var_nugget.value, type_struct.value, var_struct.value, range_struct.value, show_target.value)

for widget in [type_struct, var_struct, range_struct, var_nugget]:
    widget.observe(update_plot, names='value')

update_scenario()

display(VBox([
    HBox([type_struct, var_struct]),
    HBox([range_struct, var_nugget]),
    HBox([btn_random, btn_show_solution]),
    output,
    solution_output
]))
update_plot()