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

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

# ------------------------
# Fonctions
# ------------------------

def anisotropic_covariance(hx, hy, ranges_angles_models_weights, c0=0.0):
    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 = np.sqrt((h_rot / range_major)**2 + (v_rot / range_minor)**2)

        if model == 'Sphérique':
            base_cov = np.where(h_scaled < 1, 1 - 1.5 * h_scaled + 0.5 * h_scaled**3, 0)
        elif model == 'Exponentiel':
            base_cov = np.exp(-3 * h_scaled)
        elif model == 'Gaussien':
            base_cov = np.exp(-3 * h_scaled**2)
        elif model == 'Pépite':
            base_cov = np.where(h_scaled == 0, 1, 0)
        else:
            raise ValueError(f"Modèle inconnu : {model}")

        cov_total += c * base_cov
    cov_total += c0 * (hx == 0) * (hy == 0)
    return cov_total

def fftma_2d(nx, ny, 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, 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]

    sill_sum = c0 + sum([c for _, _, _, _, c in params])
    var_field = np.var(field)
    if var_field > 0:
        norm_factor = np.sqrt(sill_sum / var_field)
        field = field * norm_factor

    return field


def directional_variogram(field, n_samples=2000, 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

# ------------------------
# Initialisation
# ------------------------

nx, ny = 250, 250
output = Output()
field = None
gamma_exp = None
lags = None
directions_used = np.arange(0, 180, 22.5)
show_solution = False
true_params = {}

def compute_experimental_variogram():
    global gamma_exp, lags
    gamma_exp, lags = directional_variogram(
        field, n_samples=500,
        max_lag=min(nx, ny) // 2,
        lag_step=w_lag_step.value,
        directions_deg=directions_used,
        tol=w_tol.value
    )

def generate_new_field():
    global field, true_params
    angle = np.random.randint(0, 180)
    major = np.random.randint(30, 90)
    minor = np.random.randint(10, major)
    model = np.random.choice(["Sphérique", "Exponentiel", "Gaussien"])
    sill = np.round(np.random.uniform(0.5, 1.0), 2)
    nugget = np.round(np.random.uniform(0.0, 0.3), 2)
    seed = np.random.randint(0, 1_000_000)
    sim_params = [(major, minor, angle, model, sill)]
    global field
    field = fftma_2d(nx, ny, sim_params, c0=nugget, seed=seed)

    true_params.update({
        'angle': angle,
        'portée majeure': major,
        'portée mineure': minor,
        'modèle': model,
        'sill': sill,
        'nugget': nugget
    })

    compute_experimental_variogram()

# ------------------------
# Affichage
# ------------------------

def update_plot(*args):
    global show_solution
    with output:
        clear_output(wait=True)
        fig, axs = plt.subplots(3, 3, figsize=(9, 7))
        for i, angle in enumerate(directions_used[:8]):
            ax = axs[i // 3, i % 3]
            ax.plot(lags, gamma_exp[angle], 'o', label="Expérimental")

            # Modèle théorique avec anisotropie (ajusté par sliders)
            theta_rad = np.deg2rad(w_angle.value)
            cos_t, sin_t = np.cos(theta_rad), np.sin(theta_rad)

            hx = lags * np.cos(np.deg2rad(angle - w_angle.value))
            hy = lags * np.sin(np.deg2rad(angle - w_angle.value))
            h_rot = cos_t * hx + sin_t * hy
            v_rot = -sin_t * hx + cos_t * hy
            h_scaled = np.sqrt((h_rot / w_range_major.value)**2 + (v_rot / w_range_minor.value)**2)

            if w_model.value == "Sphérique":
                gamma_th = np.where(h_scaled < 1,
                                    w_c0.value + w_c1.value * (1.5 * h_scaled - 0.5 * h_scaled**3),
                                    w_c0.value + w_c1.value)
            elif w_model.value == "Exponentiel":
                gamma_th = w_c0.value + w_c1.value * (1 - np.exp(-3 * h_scaled))
            elif w_model.value == "Gaussien":
                gamma_th = w_c0.value + w_c1.value * (1 - np.exp(-3 * h_scaled**2))

            lags_with_zero = np.insert(lags, 0, 0)
            gamma_th_with_zero = np.insert(gamma_th, 0, w_c0.value)

            ax.plot(lags_with_zero, gamma_th_with_zero, '-', color='black', label='Théorique')

            if show_solution:
                angle_sol = true_params['angle']
                maj_sol = true_params['portée majeure']
                min_sol = true_params['portée mineure']
                c0_sol = true_params['nugget']
                c1_sol = true_params['sill']
                model_sol = true_params['modèle']

                theta_sol = np.deg2rad(angle_sol)
                cos_sol, sin_sol = np.cos(theta_sol), np.sin(theta_sol)
                hx_sol = lags * np.cos(np.deg2rad(angle - angle_sol))
                hy_sol = lags * np.sin(np.deg2rad(angle - angle_sol))
                h_rot_sol = cos_sol * hx_sol + sin_sol * hy_sol
                v_rot_sol = -sin_sol * hx_sol + cos_sol * hy_sol
                h_scaled_sol = np.sqrt((h_rot_sol / maj_sol)**2 + (v_rot_sol / min_sol)**2)

                if model_sol == "Sphérique":
                    gamma_true = np.where(h_scaled_sol < 1,
                                          c0_sol + c1_sol * (1.5 * h_scaled_sol - 0.5 * h_scaled_sol**3),
                                          c0_sol + c1_sol)
                elif model_sol == "Exponentiel":
                    gamma_true = c0_sol + c1_sol * (1 - np.exp(-3 * h_scaled_sol))
                elif model_sol == "Gaussien":
                    gamma_true = c0_sol + c1_sol * (1 - np.exp(-3 * h_scaled_sol**2))

                gamma_true_with_zero = np.insert(gamma_true, 0, c0_sol)

                ax.plot(lags_with_zero, gamma_true_with_zero, '--', color='red', label='Solution (réelle)')

            ax.set_title(f"{angle}°")
            ax.set_ylim(0, w_c0.value + w_c1.value + 0.5)
            ax.grid(True)
            if i == 0:
                ax.legend()

        axs[2, 2].axis('off')
        plt.tight_layout()
        plt.show()

        if show_solution:
            print("Paramètres de la solution réelle :")
            for k, v in true_params.items():
                print(f"- {k} : {v}")


# ------------------------
# Widgets
# ------------------------

w_range_major = FloatSlider(min=5, max=100, step=1, value=40, description='$a_g$', layout=Layout(width='250px'))
w_range_minor = FloatSlider(min=5, max=100, step=1, value=20, description='$a_p$', layout=Layout(width='250px'))
w_angle = IntSlider(min=0, max=180, step=1, value=45, description=r'$\theta_{a_g}$', layout=Layout(width='250px'))
w_model = Dropdown(options=["Sphérique", "Exponentiel", "Gaussien"], value="Sphérique", description="Modèle", layout=Layout(width='250px'))
w_c0 = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.1, description='$c_0$', layout=Layout(width='250px'))
w_c1 = FloatSlider(min=0.0, max=2.0, step=0.01, value=1.0, description='$c_1$', layout=Layout(width='250px'))
w_lag_step = IntSlider(min=1, max=20, step=1, value=5, description='Pas de lag', layout=Layout(width='250px'))
w_tol = IntSlider(min=1, max=30, step=1, value=10, description='Tolérance', layout=Layout(width='250px'))
w_show_solution = Checkbox(value=False, description='Afficher la solution')  # Checkbox n'a pas besoin

btn = Button(description="🎲 Nouvelle simulation", button_style="info")

def on_button_clicked(b):
    global show_solution
    generate_new_field()
    show_solution = False
    w_show_solution.value = False
    update_plot()

btn.on_click(on_button_clicked)

def on_checkbox_change(change):
    global show_solution
    show_solution = change['new']
    update_plot()

w_show_solution.observe(on_checkbox_change, names='value')

def on_variogram_param_change(change):
    compute_experimental_variogram()
    update_plot()

for w in [w_range_major, w_range_minor, w_angle, w_model, w_c0, w_c1]:
    w.observe(update_plot, names='value')

w_lag_step.observe(on_variogram_param_change, names='value')
w_tol.observe(on_variogram_param_change, names='value')

# ------------------------
# Affichage interface
# ------------------------

generate_new_field()
update_plot()

display(VBox([
    HBox([w_range_major, w_range_minor, w_angle]),
    HBox([w_model, w_c0, w_c1]),
    HBox([w_lag_step, w_tol, w_show_solution]),
    btn,
    output
]))