Atelier 6 - Calculateur de variogramme théorique (multi-composantes anisotropes)

🎯 But pédagogique¶

Calcul direct du variogramme pour tout type de structure.

Fonctionnalités :¶

Entrée : jusqu’à 3 composantes :
- Type de modèle : sphérique, exponentiel, gaussien, pépite
- Variance, portée, anisotropie ( $a_x$ , $a_y$ , angle) de chaque modèle
- Distance ( $h_x$ et $h_y$ )

💡 Note : Pour le modèle de pépite, entrée n’importe qu’elle valeur pour les paramètres du modèle.
💡 Note : Les distances peuvent être négatives.

import numpy as np
import ipywidgets as widgets
from IPython.display import display, clear_output

def cov_component(h, range_ag, range_ap, azimuth_deg, model, sill):
    theta = np.deg2rad(-(90 - azimuth_deg))  # rotation antihoraire
    
    hx, hy = h[0], h[1]
    h_rot_x = np.cos(theta) * hx - np.sin(theta) * hy
    h_rot_y = np.sin(theta) * hx + np.cos(theta) * hy
    
    h_scaled = h_rot_x / range_ag
    v_scaled = h_rot_y / range_ap
    
    dist = np.sqrt(h_scaled**2 + v_scaled**2)
    
    if model == 'Sphérique':
        if dist < 1:
            cov = sill * (1 - 1.5*dist + 0.5*dist**3)
        else:
            cov = 0
    elif model == 'Exponentiel':
        cov = sill * np.exp(-3 * dist)
    elif model == 'Gaussien':
        cov = sill * np.exp(-3 * dist**2)
    elif model == 'Pépite':
        cov = sill if dist == 0 else 0
    else:
        raise ValueError(f"Modèle inconnu : {model}")
    return cov

def covariance_multi(h, comps):
    cov_total = 0
    for c in comps:
        if c['active']:
            cov_total += cov_component(h, c['range_ag'], c['range_ap'], c['azimuth_deg'], c['model'], c['sill'])
    return cov_total

def variogram(h, comps):
    sill_total = sum(c['sill'] for c in comps if c['active'])
    return sill_total - covariance_multi(h, comps)

def create_component_widget(index):
    idx = index + 1
    active_cb = widgets.Checkbox(value=False, description=f"Comp. {idx}")
    model_dd = widgets.Dropdown(
        options=['Sphérique', 'Exponentiel', 'Gaussien', 'Pépite'],
        value='Sphérique',
        description=f'Type {idx}'
    )
    sill_txt = widgets.Text(value='1.0', description=f'Sill {idx}')
    range_ap_txt = widgets.Text(value='1.0', description=f'Portée ap {idx}')
    range_ag_txt = widgets.Text(value='1.0', description=f'Portée ag {idx}')
    azimuth_txt = widgets.Text(value='40.0', description=f'Angle {idx} (°)')
    
    def toggle_inputs(change):
        enabled = change['new']
        model_dd.disabled = not enabled
        sill_txt.disabled = not enabled
        range_ap_txt.disabled = not enabled
        range_ag_txt.disabled = not enabled
        azimuth_txt.disabled = not enabled
    
    active_cb.observe(toggle_inputs, names='value')
    toggle_inputs({'new': active_cb.value})
    
    comp_box = widgets.VBox([active_cb, model_dd, sill_txt, range_ap_txt, range_ag_txt, azimuth_txt])
    return comp_box, active_cb, model_dd, sill_txt, range_ap_txt, range_ag_txt, azimuth_txt

comp_widgets = []
for i in range(3):
    comp_widgets.append(create_component_widget(i))

hx_txt = widgets.Text(value='5.0', description='$h_x$:')
hy_txt = widgets.Text(value='0.0', description='$h_y$:')

btn_calc = widgets.Button(description="Calculer variogramme γ(h) et covariance")
output = widgets.Output()

def on_calc_clicked(b):
    with output:
        clear_output()
        warnings = []
        try:
            comps = []
            for (comp_box, active_cb, model_dd, sill_txt, rap_txt, rag_txt, az_txt) in comp_widgets:
                active = active_cb.value
                if active:
                    model = model_dd.value
                    sill = float(sill_txt.value)
                    range_ap = float(rap_txt.value)
                    range_ag = float(rag_txt.value)
                    azimuth_deg = float(az_txt.value)
                    if range_ag < range_ap:
                        warnings.append(f"⚠️ Composante {comp_box.children[0].description} : Portée ag ({range_ag}) < Portée ap ({range_ap})")
                else:
                    model = None
                    sill = 0
                    range_ap = 1
                    range_ag = 1
                    azimuth_deg = 0
                comps.append({
                    'active': active,
                    'model': model,
                    'sill': sill,
                    'range_ap': range_ap,
                    'range_ag': range_ag,
                    'azimuth_deg': azimuth_deg
                })
            
            for w in warnings:
                print(w)
            if warnings:
                print()
                
            hx = float(hx_txt.value)
            hy = float(hy_txt.value)
            h = np.array([hx, hy])
            gamma = variogram(h, comps)
            cov = covariance_multi(h, comps)
            print(f"Variogramme γ(h) pour h=({hx}, {hy}) = {gamma:.4f}")
            print(f"Covariance C(h) pour h=({hx}, {hy}) = {cov:.4f}")
        except Exception as e:
            print(f"Erreur : {e}")

btn_calc.on_click(on_calc_clicked)

ui = widgets.VBox([
    widgets.HTML("<h3>Calculateur de variogramme théorique (multi-composantes anisotropes)</h3>"),
    widgets.HBox([hx_txt, hy_txt]),
    widgets.HBox([w[0] for w in comp_widgets]),
    btn_calc,
    output
])

display(ui)

Note de cours et ateliers interactifs

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

Note de cours et ateliers interactifs

Atelier 7 - Cas cibles : erreurs et effets géologiques