Atelier 3 - Calculateur de variance de blocs - Note de cours et ateliers interactifs

🎯 But pédagogique¶

Outil informatique pour calculer la variance de bloc.

📐 Variance théorique d’un bloc carré¶

La variance d’un bloc de volume $V$ est donnée par :

\sigma_V^2 = \frac{1}{V^2} \iint_V C(\| \mathbf{r} - \mathbf{r}' \|) \, d\mathbf{r} \, d\mathbf{r}'

(1)

où :

$V$ est le volume du bloc,
$\mathbf{r}$ et $\mathbf{r}'$ sont des positions dans ce volume,
$C(h)$ est la fonction de covariance, avec $h = \| \mathbf{r} - \mathbf{r}'$ | la distance entre les deux points.

⚙️ Description des paramètres¶

Dimension (dim) : dimension spatiale du bloc (1 pour ligne, 2 pour surface, 3 pour cube).
Plateau ( $c_1$ ) : variance structurée maximale du modèle de covariance.
Effet de pépite ( $c_0$ ) : variance indépendante du modèle, représentant la variance microscopique ou bruit blanc.
Portée en X, Y, Z ( $a_x$ , $a_y$ , $a_z$ ) : paramètres de portée (range) dans chaque direction définissant la portée spatiale de la corrélation.
Longueur du bloc en X, Y, Z ( $l_x$ , $l_y$ , $l_z$ ) : dimensions physiques du bloc.
Modèle de covariance : choix du modèle théorique (sphérique, exponentiel, gaussien).

🚨 Limitations¶

Le bloc est supposé orienté selon les directions d’anisotropie définies par les portées $a_x$ , $a_y$ , et $a_z$ .
Note : Si le bloc n’est pas orienté selon ces directions, il faut appliquer une rotation pour transformer le système en un cadre isotrope.

import numpy as np
import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import FloatText, Dropdown, Layout
# --- Covariances --- #
def spherical_covariance(h, range_, sill):
    cov = np.zeros_like(h)
    mask = h <= range_
    hr = h[mask] / range_
    cov[mask] = sill * (1 - 1.5 * hr + 0.5 * hr**3)
    return cov

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

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

def get_covariance_model(name):
    return {
        "Sphérique": spherical_covariance,
        "Exponentiel": exponential_covariance,
        "Gaussien": gaussian_covariance
    }[name]

# --- Distance anisotrope --- #
def anisotropic_distance(dx, dy, dz, ax, ay, az):
    return np.sqrt((dx/ax)**2 + (dy/ay)**2 + (dz/az)**2)

# --- Calcul de variance de bloc --- #
def block_variance(dim, sill, nugget, ax, ay, az, lx, ly, lz, model_name, n_points=20):
    x = np.linspace(0, lx, n_points)
    cov_func = get_covariance_model(model_name)

    if dim == 1:
        X = x[:, None]
        dX = X - X.T
        h = anisotropic_distance(dX, 0, 0, ax, ay, az)

    elif dim == 2:
        y = np.linspace(0, ly, n_points)
        X, Y = np.meshgrid(x, y)
        coords = np.stack([X.ravel(), Y.ravel()], axis=1)
        dx = coords[:, 0][:, None] - coords[:, 0][None, :]
        dy = coords[:, 1][:, None] - coords[:, 1][None, :]
        h = anisotropic_distance(dx, dy, 0, ax, ay, az)

    elif dim == 3:
        y = np.linspace(0, ly, n_points)
        z = np.linspace(0, lz, n_points)
        X, Y, Z = np.meshgrid(x, y, z)
        coords = np.stack([X.ravel(), Y.ravel(), Z.ravel()], axis=1)
        dx = coords[:, 0][:, None] - coords[:, 0][None, :]
        dy = coords[:, 1][:, None] - coords[:, 1][None, :]
        dz = coords[:, 2][:, None] - coords[:, 2][None, :]
        h = anisotropic_distance(dx, dy, dz, ax, ay, az)

    else:
        raise ValueError("Dimension doit être 1, 2 ou 3")

    C = cov_func(h, range_=1.0, sill=sill)
    np.fill_diagonal(C, C.diagonal() + nugget)
    return np.mean(C)

# --- Widgets à saisie manuelle --- #
label_width = '120px'  # plus large que la valeur par défaut (~120px)

dim_input = Dropdown(options=[1, 2, 3], value=2, description='Dimension', style={'description_width': label_width})
sill_input = FloatText(value=1.0, description='$c_1$', style={'description_width': label_width})
nugget_input = FloatText(value=0.0, description='Effet de pépite ($c_0$)', style={'description_width': label_width})
ax_input = FloatText(value=30, description='Portée X ($a_x$)', style={'description_width': label_width})
ay_input = FloatText(value=30, description='Portée Y ($a_y$)', style={'description_width': label_width})
az_input = FloatText(value=30, description='Portée Z ($a_z$)', style={'description_width': label_width})
lx_input = FloatText(value=10, description='Longueur X ($l_x$)', style={'description_width': label_width})
ly_input = FloatText(value=10, description='Longueur Y ($l_y$)', style={'description_width': label_width})
lz_input = FloatText(value=10, description='Longueur Z ($l_z$)', style={'description_width': label_width})
model_input = Dropdown(
    options=["Sphérique", "Exponentiel", "Gaussien"], 
    value="Sphérique", 
    description='Modèle', 
    style={'description_width': label_width}
)

# --- Bouton de calcul --- #
calc_button = widgets.Button(description="Calculer la variance", button_style='success')
output_area = widgets.Output()

# --- Action quand on clique sur le bouton --- #
def on_calculate_clicked(b):
    with output_area:
        clear_output()
        var = block_variance(
            dim_input.value,
            sill_input.value,
            nugget_input.value,
            ax_input.value,
            ay_input.value,
            az_input.value,
            lx_input.value,
            ly_input.value,
            lz_input.value,
            model_input.value
        )
        print(f"▣ Variance de bloc : {var:.6f}")

calc_button.on_click(on_calculate_clicked)

# --- Affichage final --- #
ui = widgets.VBox([
    dim_input, sill_input, nugget_input,
    ax_input, ay_input, az_input,
    lx_input, ly_input, lz_input,
    model_input, calc_button, output_area
])
display(ui)

Note de cours et ateliers interactifs

Atelier 2 - Discrétisation de la variance de bloc

Note de cours et ateliers interactifs

Atelier 4 - Influence du modèle