NA Inversion Example (updated physics)
This notebook demonstrates kinematic inversion using the Neighbourhood Algorithm (NA) with the updated physical implementation:
\[M_0 = \mu(z) \times A \times slip\]
Workflow:
Load configuration and forward model
Generate controlled synthetic observed waveforms from a 7-parameter ellipse model
Run NA search
Analyze best model and compare \(M_0\)/\(M_w\)
Visualize convergence and export results
Parameters inverted:
a1, a2: Ellipse semi-axes (km)
theta: Rotation angle (x π)
np, tp: Center position
dmax: Maximum slip (m)
vr: Rupture velocity (km/s)
[ ]:
from pathlib import Path
import sys
import numpy as np
import matplotlib.pyplot as plt
def find_project_root(start: Path) -> Path:
for p in [start, *start.parents]:
if (p / 'kdellipspy').exists():
return p
raise FileNotFoundError('No se encontro PROJECT_ROOT con carpeta kdellipspy.')
PROJECT_ROOT = find_project_root(Path.cwd().resolve())
KIN_ROOT = PROJECT_ROOT / 'Kinematic_inversion'
INPUT_CTL = KIN_ROOT / 'input.ctl'
if not INPUT_CTL.exists():
raise FileNotFoundError(f'No se encontro input.ctl en {INPUT_CTL}')
if str(PROJECT_ROOT) not in sys.path:
sys.path.insert(0, str(PROJECT_ROOT))
from kdellipspy import ConfigParser
from kdellipspy import AxitraForwardModel
from kdellipspy import NAInversionModel, NAConfig, MisfitCalculator
from kdellipspy import load_and_filter_observed_data, bandpass_filter_waveforms
# Compatibilidad con celdas existentes
root = KIN_ROOT
input_ctl = INPUT_CTL
✓ All imports successful
Step 1: Load configuration
[2]:
input_ctl = root / 'input.ctl'
cfg = ConfigParser(str(input_ctl))
print(f"✓ Loaded configuration from {input_ctl}")
print(f"\nInversion parameters (from input.ctl):")
for i, param in enumerate(cfg.inversion_params.parameters, 1):
status = "INVERT" if param.flag else "FIXED"
print(f" {i}. {param.name:30s} [{param.min_val:8.3f}, {param.max_val:8.3f}] {status}")
print(f"\nInversion process parameters:")
print(f" Algorithm: {'NA' if cfg.inversion_process.algorithm_type == 0 else 'MC'}")
print(f" Iterations: {cfg.inversion_process.num_iterations}")
print(f" Initial samples: {cfg.inversion_process.ss1}")
print(f" Iteration samples: {cfg.inversion_process.ss_other}")
print(f" Resample cells: {cfg.inversion_process.cells_resample}")
✓ Loaded configuration from [PROJECT_ROOT]/Kinematic_inversion/input.ctl
Inversion parameters (from input.ctl):
1. Length of axis 1 (km) [ 5.000, 10.000] INVERT
2. Length of axis 2 (km) [ 5.000, 10.000] INVERT
3. Rotation angle (x pi) [ 0.000, 2.000] INVERT
4. Position of the center np [ 0.000, 1.000] INVERT
5. Position of the center tp (x 2pi) [ 0.000, 1.000] INVERT
6. Maximum slip (Dmax) (m) [ 1.000, 3.000] INVERT
7. Rupture velocity (Vr) (km/s) [ 0.500, 3.500] INVERT
Inversion process parameters:
Algorithm: NA
Iterations: 10
Initial samples: 100
Iteration samples: 30
Resample cells: 7
Step 2: Load REAL DATA and check it.
This apply for root / ‘DATA’ / ‘real_disp_x’
[3]:
from src.signal_utils import load_and_filter_observed_data
observed, time = load_and_filter_observed_data(
input_ctl_path=root / 'input.ctl',
data_dir=root / 'DATA'
)
print(f"✓ Observed data loaded")
print(f" Shape: {observed.shape} (n_stations={observed.shape[0]}, 3_components, npts={observed.shape[2]})")
print(f" Time range: {time[0]:.3f} - {time[-1]:.3f} s")
print(f" Sampling rate: {1/(time[1]-time[0]):.1f} Hz")
print(f" Frequency band: {float(cfg.ellipse.freq1):.3f} - {float(cfg.ellipse.freq2):.3f} Hz")
✓ Observed data loaded
Shape: (10, 3, 512) (n_stations=10, 3_components, npts=512)
Time range: 0.000 - 127.750 s
Sampling rate: 4.0 Hz
Frequency band: 0.020 - 0.100 Hz
Plot the data to visualize
[4]:
# Visualize observed data
fig, axes = plt.subplots(3, 3, figsize=(14, 8))
n_stations_show = min(3, observed.shape[0])
for i in range(n_stations_show):
axes[i, 0].plot(time, observed[i, 0], 'g-', linewidth=1.5)
axes[i, 0].set_title(f'Station {i+1} - X Component')
axes[i, 0].grid(True, alpha=0.3)
axes[i, 1].plot(time, observed[i, 1], 'g-', linewidth=1.5)
axes[i, 1].set_title(f'Station {i+1} - Y Component')
axes[i, 1].grid(True, alpha=0.3)
axes[i, 2].plot(time, observed[i, 2], 'g-', linewidth=1.5)
axes[i, 2].set_title(f'Station {i+1} - Z Component')
axes[i, 2].grid(True, alpha=0.3)
plt.suptitle('Observed Waveforms (Filtered)', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
print(f"✓ Observed data visualized ({n_stations_show} stations)")
✓ Observed data visualized (3 stations)
[6]:
# Initialize forward model and use midpoint model parameters
fwd = AxitraForwardModel(str(root / 'input.ctl'))
# Use midpoint of parameter ranges from input.ctl
midpoint_model = np.array([
0.5 * (float(p.min_val) + float(p.max_val))
for p in cfg.inversion_params.parameters
], dtype=float)
print(f"✓ Forward model initialized")
print(f"\n✓ Using midpoint model parameters:")
param_names = ['a1', 'a2', 'theta', 'np', 'tp', 'dmax', 'vr']
for name, val in zip(param_names, midpoint_model):
print(f" {name:6s} = {val:.4f}")
# Build geometry
geometry = fwd.build_geometry_with_ellipse_slip(midpoint_model)
m0, mw = fwd.estimate_total_moment_and_mw(midpoint_model, geometry)
print(f"\n✓ Geometry built:")
print(f" Total moment M0 = {m0:.3e} N.m")
print(f" Moment magnitude Mw = {mw:.2f}")
✓ Forward model initialized
axitra dir: [PROJECT_ROOT]/Kinematic_inversion/axitra
✓ Using midpoint model parameters:
a1 = 7.5000
a2 = 7.5000
theta = 1.0000
np = 0.5000
tp = 0.5000
dmax = 2.0000
vr = 2.0000
✓ Geometry built:
Total moment M0 = 1.838e+19 N.m
Moment magnitude Mw = 6.78
Step 3: Build forward model and generate synthetic waveforms
This section demonstrates the complete workflow for computing synthetic waveforms and filtering them to match the observed data’s frequency band.
[ ]:
# Create output directory
from datetime import datetime
output_dir = root / 'output' / f'na_inversion_{datetime.now().strftime("%Y%m%d_%H%M%S")}'
output_dir.mkdir(parents=True, exist_ok=True)
print(f"✓ Output directory created: {output_dir}")
# Export results using the built-in method
result.export_results(output_dir / 'inversion_results.json')
print(f"✓ Results exported to: {output_dir / 'inversion_results.json'}")
# Save additional analysis
import json
summary = {
'timestamp': datetime.now().isoformat(),
'total_models': len(result.all_models),
'iterations': int(na_config.n_iterations),
'best_misfit': float(result.best_model.misfit),
'best_iteration': int(result.best_model.iteration),
'best_model': result.best_model.model.tolist(),
'best_model_m0': float(m0_best),
'best_model_mw': float(mw_best),
'initial_model': midpoint_model.tolist(),
'initial_misfit': float(misfit_test),
'search_time_seconds': float(t_elapsed),
'misfit_improvement_percent': float(improvement),
'config': {
'freq1': float(freq1),
'freq2': float(freq2),
't0': float(cfg.ellipse.t0),
'time_window_s': 20.0,
}
}
summary_path = output_dir / 'inversion_summary.json'
with open(summary_path, 'w') as f:
json.dump(summary, f, indent=2)
print(f"✓ Summary saved to: {summary_path}")
# Save waveform comparison
waveform_data = {
'time': time.tolist(),
'observed': observed.tolist(),
'synthetic_best': synthetic_best_filtered.tolist(),
'residual': (observed - synthetic_best_filtered).tolist(),
}
waveform_path = output_dir / 'waveform_comparison.json'
with open(waveform_path, 'w') as f:
json.dump(waveform_data, f)
print(f"✓ Waveform data saved to: {waveform_path}")
# Print summary
print("\n" + "=" * 70)
print("INVERSION SUMMARY")
print("=" * 70)
print(f"Total models evaluated: {len(result.all_models)}")
print(f"Search iterations: {na_config.n_iterations}")
print(f"Search time: {t_elapsed:.1f} seconds")
print(f"\nInitial model misfit: {misfit_test:.6f}")
print(f"Best model misfit: {result.best_model.misfit:.6f}")
print(f"Improvement: {improvement:.1f}%")
print(f"\nBest model:")
for name, val in zip(param_names, result.best_model.model):
print(f" {name:6s} = {val:.4f}")
print(f"\nBest model moment:")
print(f" M0 = {m0_best:.3e} N.m")
print(f" Mw = {mw_best:.2f}")
print(f"\nResults saved to: {output_dir}")
print("=" * 70)
Step 7: Export Results
[ ]:
# Generate and visualize synthetic waveforms from best model
print("Generating synthetic waveforms from best model...")
# Build geometry with best model
geometry_best = fwd.build_geometry_with_ellipse_slip(result.best_model.model)
# Compute Green functions and convolve
axitra_best = fwd.build_axitra(geometry_best, latlon=False, freesurface=True)
ap_best = fwd.green(axitra_best, quiet=True)
result_best = fwd.conv(ap_best, geometry_best, source_type=1, t0=float(cfg.ellipse.t0), quiet=True)
if isinstance(result_best, tuple):
_, sx_best, sy_best, sz_best = result_best
synthetic_best = np.array([sx_best, sy_best, sz_best], dtype=float)
synthetic_best = np.transpose(synthetic_best, (1, 0, 2)) # (nsta, 3, npts)
else:
synthetic_best = result_best
# Apply filtering
synthetic_best_filtered = bandpass_filter_waveforms(
synthetic_best, time,
freq1=freq1, freq2=freq2,
corners=4, zerophase=True
)
print("✓ Synthetic waveforms from best model generated")
# Calculate misfit for best model
misfit_best = misfit_calc.l2_misfit(synthetic_best_filtered)
print(f" Misfit (best model): {misfit_best:.6f}")
# Visualize comparison
station_idx = 0
fig, axes = plt.subplots(2, 3, figsize=(14, 7))
components = ['X', 'Y', 'Z']
# Top row: Time series
for icomp in range(3):
ax = axes[0, icomp]
ax.plot(time, observed[station_idx, icomp], 'g-', linewidth=2, label='Observed', alpha=0.8)
ax.plot(time, synthetic_best_filtered[station_idx, icomp], 'r--', linewidth=1.5, label='Best synthetic', alpha=0.8)
ax.set_title(f'Station {station_idx+1} - {components[icomp]} Component')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.legend()
ax.grid(True, alpha=0.3)
# Bottom row: Residuals
for icomp in range(3):
ax = axes[1, icomp]
residual = observed[station_idx, icomp] - synthetic_best_filtered[station_idx, icomp]
ax.plot(time, residual, 'orange', linewidth=1.5)
ax.axhline(0, color='k', linestyle=':', alpha=0.5)
ax.set_title(f'Residual - {components[icomp]}')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.grid(True, alpha=0.3)
plt.suptitle(f'Best Model Fit (Misfit: {misfit_best:.6f})', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
print("✓ Best model visualization complete")
[ ]:
# Parameter distribution: show range explored
fig, axes = plt.subplots(2, 4, figsize=(16, 8))
axes = axes.flatten()
for iparam in range(len(param_names)):
param_name = param_names[iparam]
param_vals = all_models[:, iparam]
param_min = float(cfg.inversion_params.parameters[iparam].min_val)
param_max = float(cfg.inversion_params.parameters[iparam].max_val)
ax = axes[iparam]
ax.hist(param_vals, bins=20, alpha=0.7, color='steelblue', edgecolor='black')
ax.axvline(result.best_model.model[iparam], color='red', linestyle='--', linewidth=2, label='Best')
ax.axvline(midpoint_model[iparam], color='orange', linestyle=':', linewidth=2, label='Initial')
ax.set_xlim([param_min, param_max])
ax.set_xlabel(param_name)
ax.set_ylabel('Frequency')
ax.set_title(f'{param_name} distribution (n={len(param_vals)})')
ax.legend()
ax.grid(True, alpha=0.3, axis='y')
# Hide extra subplot
axes[7].axis('off')
plt.suptitle('Parameter Space Exploration', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
print("✓ Parameter distribution plot complete")
[ ]:
# Extract all models and misfits for analysis
all_models = np.array([m.model for m in result.all_models])
all_misfits = np.array([m.misfit for m in result.all_models])
all_iterations = np.array([m.iteration for m in result.all_models])
print(f"✓ Extracted {len(all_models)} models from search")
print(f" Misfit range: [{all_misfits.min():.6f}, {all_misfits.max():.6f}]")
print(f" Iterations: {all_iterations.min()} to {all_iterations.max()}")
# Plot convergence
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
# Misfit convergence (all models)
axes[0, 0].scatter(range(len(all_misfits)), all_misfits, alpha=0.6, s=30)
axes[0, 0].axhline(result.best_model.misfit, color='r', linestyle='--', linewidth=2, label='Best')
axes[0, 0].set_xlabel('Model Index')
axes[0, 0].set_ylabel('Misfit (L2)')
axes[0, 0].set_title('Convergence: All Models')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# Misfit vs iteration
axes[0, 1].scatter(all_iterations, all_misfits, alpha=0.6, s=30, c=all_iterations, cmap='viridis')
axes[0, 1].set_xlabel('Iteration')
axes[0, 1].set_ylabel('Misfit (L2)')
axes[0, 1].set_title('Misfit by Iteration')
axes[0, 1].grid(True, alpha=0.3)
# Best misfit evolution by iteration
best_per_iteration = []
for it in range(int(all_iterations.max()) + 1):
mask = all_iterations == it
if mask.any():
best_per_iteration.append(np.min(all_misfits[mask]))
axes[1, 0].plot(range(len(best_per_iteration)), best_per_iteration, 'o-', linewidth=2, markersize=8)
axes[1, 0].set_xlabel('Iteration')
axes[1, 0].set_ylabel('Best Misfit in Iteration')
axes[1, 0].set_title('Best Misfit Evolution')
axes[1, 0].grid(True, alpha=0.3)
# Parameter space: best models
n_best = min(20, len(result.all_models))
best_indices = np.argsort(all_misfits)[:n_best]
best_models_subset = all_models[best_indices]
best_misfits_subset = all_misfits[best_indices]
axes[1, 1].scatter(best_models_subset[:, 0], best_models_subset[:, 1],
c=best_misfits_subset, cmap='RdYlGn_r', s=50, alpha=0.8)
axes[1, 1].scatter(result.best_model.model[0], result.best_model.model[1],
color='red', s=200, marker='*', label='Best model', edgecolor='black', linewidth=2)
axes[1, 1].set_xlabel('a1 (km)')
axes[1, 1].set_ylabel('a2 (km)')
axes[1, 1].set_title('Parameter Space (a1 vs a2)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle('NA Inversion Convergence Analysis', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
print("✓ Convergence analysis plot complete")
Step 6: Analyze Results and Visualize
Examine the inversion results and visualize convergence and parameter distribution.
[ ]:
# Run NA search
print("=" * 70)
print("STARTING NEIGHBOURHOOD ALGORITHM SEARCH")
print("=" * 70)
import time as time_module
t0 = time_module.time()
result = na_inversion.run_na_search(na_config=na_config)
t_elapsed = time_module.time() - t0
print("=" * 70)
print("NA SEARCH COMPLETE")
print("=" * 70)
print(f"\nSearch completed in {t_elapsed:.1f} seconds")
print(f"\n✓ Best model found:")
print(f" Misfit: {result.best_model.misfit:.6f}")
print(f" Iteration: {result.best_model.iteration}")
# Show best model parameters
print(f"\n✓ Best model parameters:")
best_params = result.best_model.model
for name, val in zip(param_names, best_params):
print(f" {name:6s} = {val:.4f}")
# Estimate moment magnitude for best model
m0_best, mw_best = fwd.estimate_total_moment_and_mw(best_params)
print(f"\n✓ Best model moment:")
print(f" M0 = {m0_best:.3e} N.m")
print(f" Mw = {mw_best:.2f}")
# Compare with initial midpoint model
print(f"\n✓ Comparison: Initial vs Best")
print(f" Initial misfit: {misfit_test:.6f}")
print(f" Best misfit: {result.best_model.misfit:.6f}")
improvement = (misfit_test - result.best_model.misfit) / misfit_test * 100
print(f" Improvement: {improvement:.1f}%")
======================================================================
STARTING NEIGHBOURHOOD ALGORITHM SEARCH
======================================================================
[NA] Adjusting n_samples_iteration from 30 to 35 to satisfy neighpy constraint (multiple of n_cells_resample=7).
[NA] Starting search: ni=100, ns=35, n=10, nr=7 -> expected evaluations=450
NAI - Initial Random Search
=========================
[NA] iter=000 eval=00001 misfit=2.769579e+02 best=2.769579e+02
=========================
[NA] iter=000 eval=00002 misfit=7.618012e+02 best=2.769579e+02
=========================
[NA] iter=000 eval=00003 misfit=2.285451e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00004 misfit=4.513069e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00005 misfit=6.275357e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00006 misfit=4.133923e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00007 misfit=4.572672e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00008 misfit=6.211029e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00009 misfit=6.818219e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00010 misfit=5.649620e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00011 misfit=3.688018e+02 best=2.285451e+02
=========================
[NA] iter=000 eval=00012 misfit=1.902730e+02 best=1.902730e+02
=========================
[NA] iter=000 eval=00013 misfit=1.746109e+02 best=1.746109e+02
=========================
[NA] iter=000 eval=00014 misfit=5.637256e+02 best=1.746109e+02
=========================
[NA] iter=000 eval=00015 misfit=1.008369e+03 best=1.746109e+02
=========================
[NA] iter=000 eval=00016 misfit=8.488220e+01 best=8.488220e+01
=========================
[NA] iter=000 eval=00017 misfit=2.297137e+02 best=8.488220e+01
=========================
[NA] iter=000 eval=00018 misfit=2.990200e+02 best=8.488220e+01
=========================
[NA] iter=000 eval=00019 misfit=3.661088e+02 best=8.488220e+01
=========================
[NA] iter=000 eval=00020 misfit=1.361535e+02 best=8.488220e+01
=========================
[NA] iter=000 eval=00021 misfit=8.510534e+02 best=8.488220e+01
=========================
[NA] iter=000 eval=00022 misfit=1.214458e+02 best=8.488220e+01
=========================
[NA] iter=000 eval=00023 misfit=6.576112e+01 best=6.576112e+01
=========================
[NA] iter=000 eval=00024 misfit=4.581872e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00025 misfit=7.365747e+01 best=6.576112e+01
=========================
[NA] iter=000 eval=00026 misfit=1.512020e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00027 misfit=4.657471e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00028 misfit=9.151568e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00029 misfit=1.917657e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00030 misfit=4.426384e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00031 misfit=1.443503e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00032 misfit=4.627750e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00033 misfit=7.229074e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00034 misfit=4.095387e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00035 misfit=2.421097e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00036 misfit=3.756428e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00037 misfit=1.435802e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00038 misfit=1.147086e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00039 misfit=6.602224e+01 best=6.576112e+01
=========================
[NA] iter=000 eval=00040 misfit=8.644089e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00041 misfit=4.220766e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00042 misfit=8.174270e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00043 misfit=4.995017e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00044 misfit=9.882750e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00045 misfit=6.887415e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00046 misfit=4.356548e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00047 misfit=4.201063e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00048 misfit=3.499546e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00049 misfit=9.236685e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00050 misfit=4.875200e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00051 misfit=7.304803e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00052 misfit=2.215242e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00053 misfit=5.268643e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00054 misfit=1.043245e+03 best=6.576112e+01
=========================
[NA] iter=000 eval=00055 misfit=1.671066e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00056 misfit=3.826328e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00057 misfit=6.115258e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00058 misfit=8.224414e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00059 misfit=1.786400e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00060 misfit=9.329408e+01 best=6.576112e+01
=========================
[NA] iter=000 eval=00061 misfit=9.415829e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00062 misfit=7.946245e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00063 misfit=3.471845e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00064 misfit=1.377649e+03 best=6.576112e+01
=========================
[NA] iter=000 eval=00065 misfit=2.057245e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00066 misfit=9.334342e+01 best=6.576112e+01
=========================
[NA] iter=000 eval=00067 misfit=3.017572e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00068 misfit=3.656392e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00069 misfit=2.566395e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00070 misfit=1.590060e+03 best=6.576112e+01
=========================
[NA] iter=000 eval=00071 misfit=4.428784e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00072 misfit=2.531554e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00073 misfit=5.598516e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00074 misfit=6.005383e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00075 misfit=7.093158e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00076 misfit=6.846966e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00077 misfit=7.099600e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00078 misfit=1.378297e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00079 misfit=3.805219e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00080 misfit=1.052804e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00081 misfit=1.283322e+03 best=6.576112e+01
=========================
[NA] iter=000 eval=00082 misfit=2.716042e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00083 misfit=1.764444e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00084 misfit=2.277101e+03 best=6.576112e+01
=========================
[NA] iter=000 eval=00085 misfit=3.485750e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00086 misfit=5.811730e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00087 misfit=1.424227e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00088 misfit=2.894929e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00089 misfit=3.624547e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00090 misfit=2.402042e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00091 misfit=6.369413e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00092 misfit=1.045446e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00093 misfit=3.678291e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00094 misfit=4.154266e+02 best=6.576112e+01
=========================
[NA] iter=000 eval=00095 misfit=6.264214e+02 best=6.576112e+01
[32]:
# Configure NA search
na_config = NAConfig(
n_samples_initial=int(cfg.inversion_process.ss1), # Initial samples
n_samples_iteration=int(cfg.inversion_process.ss_other), # Samples per iteration
n_iterations=int(cfg.inversion_process.num_iterations), # Number of iterations
n_cells_resample=int(cfg.inversion_process.cells_resample),
n_jobs=-1,
random_seed=42,
keep_axitra_files=False,
)
print("✓ NA Configuration:")
print(f" Initial samples: {na_config.n_samples_initial}")
print(f" Samples per iteration: {na_config.n_samples_iteration}")
print(f" Iterations: {na_config.n_iterations}")
print(f" Cells to resample: {na_config.n_cells_resample}")
print(f" Total models: {na_config.n_samples_initial + na_config.n_samples_iteration * na_config.n_iterations}")
# Initialize NA inversion model
na_inversion = NAInversionModel(
input_ctl_path=str(root / 'input.ctl'),
,
observed_waveforms=observed,
time_array=time,
azi_times_array=azi_times,
)
print("✓ NA inversion model initialized")
✓ NA Configuration:
Initial samples: 100
Samples per iteration: 30
Iterations: 10
Cells to resample: 7
Total models: 400
✓ NA inversion model initialized
Step 5: Run Neighbourhood Algorithm (NA) Inversion
Configure and execute the NA search to find the best model parameters.
[29]:
# Build azi_times array and initialize MisfitCalculator
from src.signal_utils import build_azi_times_array
print("Building azi_times array...")
azi_times = build_azi_times_array(input_ctl_path=root / 'input.ctl')
print(f"✓ azi_times array built: shape {azi_times.shape}")
# Initialize misfit calculator
misfit_calc = MisfitCalculator(
observed_waveforms=observed,
time_array=time,
azi_times_array=azi_times,
time_window_s=20.0,
)
print("✓ MisfitCalculator initialized")
# Calculate misfit for the synthetic waveforms
misfit_test = misfit_calc.l2_misfit(synthetic_filtered)
print(f"\n✓ Misfit calculation (synthetic vs observed):")
print(f" L2 misfit: {misfit_test:.6f}")
# Show diagnostics
print("\n--- Misfit Diagnostics ---")
diag_text = misfit_calc.diagnostics_summary(synthetic_filtered, max_stations=3)
print(diag_text)
Building azi_times array...
✓ azi_times array built: shape (10, 3)
✓ MisfitCalculator initialized
✓ Misfit calculation (synthetic vs observed):
L2 misfit: 1.026894
--- Misfit Diagnostics ---
[MISFIT DIAG] rms_global(obs)=1.971e-04 rms_global(syn)=4.901e-06 syn/obs=2.487e-02
[MISFIT DIAG] sta=01 P(R): obs=1.152e-04 syn=3.150e-05 | P(Z): obs=1.359e-04 syn=2.343e-05 | S(T): obs=5.550e-04 syn=4.990e-06
[MISFIT DIAG] sta=02 P(R): obs=2.337e-04 syn=2.002e-05 | P(Z): obs=8.872e-05 syn=9.913e-06 | S(T): obs=1.207e-04 syn=1.922e-05
[MISFIT DIAG] sta=03 P(R): obs=1.309e-04 syn=9.204e-06 | P(Z): obs=1.160e-04 syn=5.515e-06 | S(T): obs=1.258e-04 syn=8.568e-06
[MISFIT DIAG] window_energy(obs)=4.150e-05 window_energy(syn)=2.329e-07 syn/obs=5.611e-03
Step 4: Compute Misfit and Build Misfit Calculator
Now calculate the misfit between observed and filtered synthetic waveforms.