🧪 Atelier 2 – Exploration de la nuée variographique

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft2, ifft2
from scipy.spatial.distance import pdist
from ipywidgets import IntText, FloatText, VBox, Button, Output, Tab
import matplotlib.gridspec as gridspec
from IPython.display import display, clear_output

class VariogramWorkshop:
    def __init__(self):
        # Widgets paramètres (recalcul)
        self.n_samples_widget = IntText(value=100, description="Nb échantillons")
        self.range_widget = FloatText(value=30.0, description="Portée")
        self.sill_widget = FloatText(value=1.0, description="Sill")
        self.nugget_widget = FloatText(value=0.1, description="Nugget")
        self.bin_width_widget = FloatText(value=10.0, description="Bin (m)")
        self.calc_button = Button(description="Calculer", button_style='success')

        # Bouton suivant
        self.next_button = Button(description="Suivant", button_style='info')
        self.next_button.disabled = True

        # Outputs
        self.out_plot = Output()
        self.out_text = Output()

        # Onglets
        self.tabs = Tab(children=[self.out_plot, self.out_text])
        self.tabs.set_title(0, "Visualisation")
        self.tabs.set_title(1, "Calculs")

        # Data placeholders
        self.sample_coords = None
        self.sample_values = None
        self.dists = None
        self.gamma = None
        self.n_samples = 0
        self.field = None
        self.pair_index = 0

        # Connexions
        self.calc_button.on_click(self.on_calculate)
        self.next_button.on_click(self.on_next_pair)

        # Affichage interface
        param_box = VBox([
            self.n_samples_widget, self.range_widget, self.sill_widget,
            self.nugget_widget, self.bin_width_widget, self.calc_button])
        full_ui = VBox([param_box, self.next_button, self.tabs])
        display(full_ui)

    def generate_gaussian_field_fftma(self, size, range_, sill, nugget, seed=0):
        np.random.seed(seed)
        extended_size = (2 * size[0], 2 * size[1])
        x = np.arange(extended_size[0])
        y = np.arange(extended_size[1])
        X, Y = np.meshgrid(x, y, indexing='ij')
        h = np.sqrt((X - size[0])**2 + (Y - size[1])**2)
        cov = sill * np.exp(-h / range_)
        cov[0, 0] += nugget
        spectrum = np.abs(fft2(cov))**0.5
        white_noise = np.random.normal(size=extended_size)
        field = np.real(ifft2(fft2(white_noise) * spectrum))
        field_cropped = field[:size[0], :size[1]]
        std = np.std(field_cropped)
        return field_cropped / std * np.sqrt(sill+nugget) if std > 0 else field_cropped

    def pdist_to_indices(self, k, n):
        i = n - 2 - int(np.floor(np.sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5))
        j = k + i + 1 - n*(n-1)//2 + (n - i)*((n - i) -1)//2
        return i, j

    def on_calculate(self, b):
        # Récupération des paramètres de l'interface
        n_samples = self.n_samples_widget.value
        range_ = self.range_widget.value
        sill = self.sill_widget.value
        nugget = self.nugget_widget.value
        bin_width = self.bin_width_widget.value

        # Taille de la grille (200x200 ici)
        size = (200, 200)

        # Génération d'un nouveau champ gaussien
        self.field = self.generate_gaussian_field_fftma(size, range_, sill, nugget, seed=12345)

        # Sélection aléatoire des échantillons à partir de la grille
        sample_indices = np.random.choice(size[0]*size[1], n_samples, replace=False)
        self.sample_coords = np.array(np.unravel_index(sample_indices, size)).T
        self.sample_values = self.field[self.sample_coords[:, 0], self.sample_coords[:, 1]]
        self.n_samples = n_samples

        # Calcul des distances entre les échantillons
        self.dists = pdist(self.sample_coords)
    
        # Calcul des différences et de la demi-variance (γ)
        diffs = pdist(self.sample_values[:, None])
        self.gamma = 0.5 * diffs**2

        # Calcul du nombre maximal de paires d'échantillons
        self.max_pairs = n_samples * (n_samples - 1) // 2
        self.pair_index = 0

        # Déverrouiller le bouton suivant
        self.next_button.disabled = False

        # Mettre à jour les variables internes
        self.range_ = range_
        self.sill = sill
        self.nugget = nugget
        self.bin_width = bin_width
        self.size = size

        # Mettre à jour l'affichage du graphique
        self.update_plot()


    def update_plot(self):
        pair_index = self.pair_index
        n_samples = self.n_samples

        i, j = self.pdist_to_indices(pair_index, n_samples)
        p1 = self.sample_coords[i]
        p2 = self.sample_coords[j]
        z1 = self.sample_values[i]
        z2 = self.sample_values[j]
        dist_p = self.dists[pair_index]
        gamma_p = self.gamma[pair_index]

        # Binning variogramme expérimental
        max_dist = 200
        bins = np.arange(0, max_dist + self.bin_width, self.bin_width)
        bin_centers = 0.5 * (bins[:-1] + bins[1:])
        bin_indices = np.digitize(self.dists, bins) - 1
        experimental = [np.mean(self.gamma[bin_indices == k]) if np.any(bin_indices == k) else np.nan for k in range(len(bin_centers))]

        # Modèle théorique
        h_vals = np.linspace(0, max_dist, 200)
        theoretical = self.nugget + self.sill * (1 - np.exp(-h_vals / self.range_))

        # Affichage graphique
        with self.out_plot:
            clear_output(wait=True)
            fig = plt.figure(figsize=(14, 5))
            gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1.2])

            ax0 = fig.add_subplot(gs[0])
            ax0.imshow(self.field, cmap='viridis', origin='lower')
            ax0.scatter(self.sample_coords[:, 1], self.sample_coords[:, 0], c='r', s=10, label='Échantillons')
            ax0.scatter([p1[1], p2[1]], [p1[0], p2[0]], c='yellow', s=50, edgecolors='black', label='Points choisis')
            ax0.set_title("Champ Gaussien 2D")
            ax0.legend(loc='upper right')

            ax1 = fig.add_subplot(gs[1])
            ax1.scatter(self.dists, self.gamma, alpha=0.3, label='Nuée variographique', color='cornflowerblue')
            ax1.plot(h_vals, theoretical, 'r-', label='Modèle théorique')
            ax1.plot(bin_centers, experimental, 'ko-', label='Variogramme expérimental')
            for b in bins:
                ax1.axvline(b, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
            ax1.scatter(dist_p, gamma_p, c='red', s=100, label='Couple choisi')
            ax1.set_xlim(0, 200)
            ax1.set_xlabel("Distance h (m)")
            ax1.set_ylabel("½(zᵢ - zⱼ)²")
            ax1.set_title("Nuée variographique")
            ax1.legend(loc='upper right')

            plt.tight_layout()
            plt.show()

        # Affichage texte dans l'onglet dédié
        with self.out_text:
            clear_output(wait=True)
            txt = (
                f"Point 1 (y,x): {p1}, valeur z₁ = {z1:.3f}\n"
                f"Point 2 (y,x): {p2}, valeur z₂ = {z2:.3f}\n"
                f"Distance = ||p₁ - p₂|| = {dist_p:.2f}\n"
                f"Calcul demi-variance : γ(h) = ½ (z₁ - z₂)² = ½ * ({z1:.3f} - {z2:.3f})² = {gamma_p:.3f}"
            )
            print(txt)

    def on_next_pair(self, b):
        self.pair_index += 1
        if self.pair_index >= self.max_pairs:
            self.pair_index = 0
        self.update_plot()

# Lancer l'interface
VariogramWorkshop()