"""Cross-check the NumPy SBAS inversion in insar_subsidence.py with MintPy.

insar_subsidence.py writes the referenced LiCSAR stack to data/ifgramStack.h5 in MintPy's layout.
Here MintPy inverts the same network (ifgram_inversion), fits the velocity (timeseries2velocity)
and we compare the result pixel by pixel with our NumPy version.
Run with:  python mintpy_crosscheck.py   (after insar_subsidence.py)
"""
import contextlib
import os
from pathlib import Path

import h5py
import mintpy
import numpy as np
from mintpy.cli import ifgram_inversion, timeseries2velocity

HERE = Path(__file__).parent
DATA = HERE / "data"
os.chdir(DATA)
print(f"MintPy {mintpy.__version__}")

# Unweighted least squares: the same estimator as np.linalg.lstsq. MintPy's progress log goes to mintpy.log.
with open("mintpy.log", "w") as log, contextlib.redirect_stdout(log):
    ifgram_inversion.main(["ifgramStack.h5", "--weight-func", "no", "--update"])
    timeseries2velocity.main(["timeseries.h5", "-o", "velocity.h5", "--update"])

with h5py.File("timeseries.h5") as f:
    ts = f["timeseries"][:]                    # metres, positive = towards the satellite
    dates = [d.decode() for d in f["date"][:]]
with h5py.File("velocity.h5") as f:
    v_mintpy = f["velocity"][:] * 1000         # m/yr -> mm/yr
mask = np.load("mask.npy")
v_numpy = np.load("v_los_numpy.npy")
diff = v_mintpy[mask] - v_numpy[mask]
print(f"\nMintPy time series: {ts.shape[0]} dates, {dates[0]} to {dates[-1]}")
print(f"LOS velocity, MintPy: 1st percentile {np.percentile(v_mintpy[mask], 1):.0f} mm/yr, "
      f"median {np.median(v_mintpy[mask]):.0f} mm/yr")
print(f"MintPy minus NumPy over {mask.sum():,} pixels: mean {diff.mean():+.2f} mm/yr, "
      f"max |difference| {np.abs(diff).max():.2f} mm/yr, correlation {np.corrcoef(v_mintpy[mask], v_numpy[mask])[0, 1]:.4f}")