x2sys_solve: warn and offer ridge regularization for unstable drift fits

[Esteban82]

Description of proposed changes

x2sys_solve -Ed (or -Et) can return per-track offsets hundreds of metres in size on a survey only a few hundred metres deep. Applying that correction table doesn't clean the data, it destroys it — and nothing warns you, since the crossover residuals stay small and the fit still looks fine.

The cause is conditioning, not a coding error: these models fit two parameters per track, and when a track's crossings are few or clustered together, the offset is barely constrained, so least squares is free to return a huge value that still satisfies the crossovers.

This PR adds:

1. A warning when a track's solved offset exceeds 20x the pre-correction crossover-error standard deviation.
2. An opt-in -E<mode>+r[<K>] modifier for -Ed/-Et that ridge-regularizes the solve, pulling poorly-constrained tracks back towards a plain offset while leaving well-constrained ones alone. Default K = 10; smaller K damps harder.

Implementation notes

Standard Tikhonov/ridge fix: before the gmt_gaussjordan solve, add lambda to the diagonal of N for each real unknown (not the Lagrange rows enforcing the per-cluster offset-sum-to-zero gauge):

lambda = mean(diag(N)) / K^2

Scaling to the mean diagonal keeps it invariant to data/weight units. K is a trust multiplier on |offset| relative to the pre-correction COE st.dev., default 10; smaller K damps harder. Well-constrained tracks (N[i][i] >> lambda) are essentially unaffected.

Done with Sonnet 5 and reviewed with Opus 5

Tested with

I personally tested it with my own data, and it works fine. Below is the script to test it with synthetic data.

Full script:

# x2sys_solve -Ed returns per-track offsets far larger than the crossover errors
# they are meant to remove, when the crossings cannot separate offset from drift.
# Fully synthetic, no external data. Deterministic (srand 42).
#
# Geometry: 10 N-S lines (100 km long) crossed by 3 E-W tie lines clustered
# between y=58 and y=62 km => 3 crossings per N-S line, all within 4 km of each
# other, so "offset at d=0" is a 58 km extrapolation from a 4 km baseline.
# Each track carries a known constant error of a few metres, plus +-3 m noise.
#
# Produces two figures: BUG_x2sys_drift.png (solved offset per track) and
# BUG_x2sys_drift_map.png (same bug/fix, gridded into an actual map).

rm -rf tracks SYNTH
mkdir -p tracks
export X2SYS_HOME=$PWD
region=0/110000/0/100000

awk 'BEGIN{OFS="\t"; srand(42); split("58000 60000 62000",Y," ")
	for (k=1;k<=10;k++) {x=k*10000; e=(k%2?2.0:-3.0); print "NS"sprintf("%02d",k), e > "truth.txt"
		f=sprintf("tracks/NS%02d.txt",k)
		for (y=0;y<=100000;y+=1000) printf "%.3f\t%.3f\t%.3f\t%.1f\n", x,y,-(200+3*x/1000+0.5*y/1000)+e+3*(rand()-0.5)*2, y > f}
	for (m=1;m<=3;m++) {y=Y[m]; e=(m==2?4.0:-1.5); print "EW"sprintf("%02d",m), e > "truth.txt"
		f=sprintf("tracks/EW%02d.txt",m)
		for (x=0;x<=110000;x+=1000) printf "%.3f\t%.3f\t%.3f\t%.1f\n", x,y,-(200+3*x/1000+0.5*y/1000)+e+3*(rand()-0.5)*2, x > f}}'

(cd tracks && ls *.txt) > tracks.lis

cat > xyzd.fmt <<'EOF'
# x, y, z, dist (Cartesian)
# ASCII
#name  type  NaN  NaN-proxy  scale  offset  oformat
x      a     N    0          1      0       -
y      a     N    0          1      0       -
z      a     N    0          1      0       -
dist   a     N    0          1      0       -
EOF

gmt x2sys_init SYNTH -Etxt -F -R$region -I50000 -Dxyzd
echo tracks >> SYNTH/SYNTH_paths.txt
gmt x2sys_cross -TSYNTH =tracks.lis -Qe -Ia > outs.txt
gmt x2sys_list -TSYNTH outs.txt -Cz -Fndc > coe_ndc.txt
gmt x2sys_list -TSYNTH outs.txt -Cz -Fnc  > coe_nc.txt

gmt x2sys_solve -TSYNTH coe_ndc.txt -Cz -Ed    > sol_Ed.txt	# unstable, no fix
gmt x2sys_solve -TSYNTH coe_ndc.txt -Cz -Ed+r  > sol_Edr.txt	# regularized, the fix
gmt x2sys_solve -TSYNTH coe_nc.txt  -Cz -Ec    > sol_Ec.txt	# reference

# --- Figure 1: solved offset per track vs truth ---

# solved offsets -> "index offset", one file per model, tracks in a common order.
# The gauge is "offsets sum to zero", so compare against truth minus its mean.
for f in Ed Edr Ec; do
	sort -k1,1 sol_$f.txt | awk -F'\t' '{gsub(/ /,"",$1); print NR, $3}' > p_$f.txt
done
sort -k1,1 truth.txt | awk '{n++; v[n]=$2; s+=$2} END{for(i=1;i<=n;i++) print i, v[i]-s/n}' > p_truth.txt

gmt begin BUG_x2sys_drift png
	gmt subplot begin 2x1 -Fs16c/7c -M0.2c/1.1c -Scol -T"x2sys_solve: solved offset per track (13 synthetic tracks, 3 clustered crossings each)"
	gmt subplot set 0
		gmt basemap -R0.5/13.5/-600/600 -JX16c/7c -Baf -BWSne+t"full range" -By+l"solved offset (m)"
		gmt plot p_Ed.txt   -Sc0.22c -Gred    -W0.25p -l"-Ed (unstable)"
		gmt plot p_Edr.txt  -Sc0.16c -Gdodgerblue    -l"-Ed ridge-regularized"
		gmt plot p_Ec.txt   -Sc0.10c -Gdarkgreen     -l"-Ec (reference)"
		gmt plot p_truth.txt -Sx0.25c -W1p,black     -l"true injected error"
		gmt legend -DjTR+o0.3c/0.3c+w4.8c -F+p0.5p+gwhite
	gmt subplot set 1
		gmt basemap -R0.5/13.5/-15/15 -JX16c/7c -Baf -BWSne+t"zoom to +/-15 m" -Bx+l"track #" -By+l"solved offset (m)"
		gmt plot p_Ed.txt   -Sc0.22c -Gred -W0.25p
		gmt plot p_Edr.txt  -Sc0.16c -Gdodgerblue
		gmt plot p_Ec.txt   -Sc0.10c -Gdarkgreen
		gmt plot p_truth.txt -Sx0.25c -W1p,black
	gmt subplot end
gmt end

# --- Figure 2: same bug/fix, gridded into an actual map ---

# Apply each correction table and grid the result -- same recipe (x2sys_datalist -L,
# blockmedian, surface -T0.35) as the real processing pipeline that first hit this bug.
for f in Ed Edr Ec; do
	gmt x2sys_datalist -TSYNTH -Lsol_$f.txt =tracks.lis -Fx,y,z > xyz_$f.txt
	gmt blockmedian xyz_$f.txt -R$region -I2000 > bm_$f.txt
	gmt surface bm_$f.txt -R$region -I2000 -Ggrid_$f.nc -T0.35
done

gmt begin BUG_x2sys_drift_map png
	gmt makecpt -Cbathy -T-800/-100
	gmt subplot begin 1x3 -Fs8c/7.3c -M0.3c/0.5c -R$region -JX8c/7.3c
	gmt subplot set 0
		gmt grdgradient grid_Ec.nc -Ne0.75 -A155 -Ggrid_Ec_int.nc
		gmt grdimage grid_Ec.nc -Igrid_Ec_int.nc -C -Baf -BWSne+t"-Ec (reference)"
	gmt subplot set 1
		gmt grdgradient grid_Ed.nc -Ne0.75 -A155 -Ggrid_Ed_int.nc
		gmt grdimage grid_Ed.nc -Igrid_Ed_int.nc -C -Baf -BWSne+t"-Ed (sin fix / no fix)"
	gmt subplot set 2
		gmt grdgradient grid_Edr.nc -Ne0.75 -A155 -Ggrid_Edr_int.nc
		gmt grdimage grid_Edr.nc -Igrid_Edr_int.nc -C -Baf -BWSne+t"-Ed+r (con fix / fixed)"
	gmt subplot end
	gmt colorbar -C -DJBC+o0/0.8c+w16c/0.4c -Baf -B+l"z (m)"
gmt end

# Clean up everything but the script and the two figures
rm -rf tracks SYNTH tracks.lis xyzd.fmt outs.txt coe_n*.txt sol_E*.txt p_E*.txt p_truth.txt truth.txt xyz_E*.txt bm_E*.txt grid_E*.nc gmt.history

Nine of the thirteen tracks are flagged by the new warning, e.g.:

x2sys_solve [WARNING]: Track NS01: solved offset 518.4 is 111x the pre-correction
COE st.dev. (4.677), backed by only 3 crossing(s) -- likely unstable/unreliable,
not a real systematic offset. Try appending +r to this -E option, or fall back to -Ec.