diff --git a/docs/gpu.md b/docs/gpu.md index 7c9b524e5..eafcd5235 100644 --- a/docs/gpu.md +++ b/docs/gpu.md @@ -1,6 +1,11 @@ -# Configure GPU acceleration for the network inversion # +# Configure GPU acceleration # -The `invert_network` step (in `ifgram_inversion.py`) ships an opt-in GPU solver that batches the per-pixel weighted least-squares inversion as normal-equations + Cholesky on a CUDA device via PyTorch. This is a partial GPU implementation: only `invert_network` is offloaded to the GPU; every other step in `smallbaselineApp.py` continues to run on the CPU. The solver is opt-in — the default `mintpy.networkInversion.solver = auto` resolves to `cpu`, so existing setups are unaffected. +Two `smallbaselineApp.py` steps ship an opt-in GPU solver that batches their per-pixel inversion as normal-equations + Cholesky on a CUDA device via PyTorch: + +- **`invert_network`** (`ifgram_inversion.py`) — weighted least-squares network inversion. Toggle: `mintpy.networkInversion.solver = torch` / `--solver torch`. +- **`correct_topography`** (`dem_error.py`) — pixelwise DEM-error fit. Toggle: `mintpy.topographicResidual.solver = torch` / `--solver torch`. Only the pixelwise-geometry branch (`pixelwiseGeometry = yes`, the production default) is GPU-dispatched; the mean-geometry branch is already pixel-batched on CPU and stays there. + +Every other step in `smallbaselineApp.py` continues to run on the CPU. Both solvers are opt-in — `mintpy.*.solver = auto` resolves to `cpu`, so existing setups are unaffected. The `torch` solver is orthogonal to Dask parallel processing (see [dask.md](./dask.md)): the former replaces the per-pixel CPU loop with a single batched Cholesky on one CUDA device, the latter distributes that same per-pixel loop across multiple worker processes. The two paths are not currently combined; pick one. @@ -8,7 +13,7 @@ The `torch` solver is orthogonal to Dask parallel processing (see [dask.md](./da See [installation.md](./installation.md) section 2.4 for installing the `[gpu]` extras with the matching CUDA wheel index. Selecting `solver = torch` on a host without a visible CUDA device is a hard error (no silent CPU fallback). -## 2. Enable ## +## 2. Enable on `invert_network` ## #### 2.1 via command line #### @@ -49,7 +54,41 @@ ifgram_inversion.py inputs/ifgramStack.h5 -w no --solver torch The two outputs should agree to float32 round-off (RMS on the order of 1e-5). -## 3. Behavior notes ## +## 3. Enable on `correct_topography` ## + +The same `--solver torch` opt-in is exposed on `dem_error.py`, controlled by `mintpy.topographicResidual.solver`. Only the pixelwise-geometry branch (`pixelwiseGeometry = yes`) is GPU-dispatched. + +#### 3.1 via command line #### + +```bash +dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver torch +dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver torch --gpu-chunk-size 200000 +``` + +#### 3.2 via template file #### + +```cfg +mintpy.topographicResidual.solver = torch #[cpu / torch], auto for cpu +mintpy.topographicResidual.gpuChunkSize = auto #[int >= 0], auto for 0 (auto-size from free VRAM) +``` + +then: + +```bash +smallbaselineApp.py smallbaselineApp.cfg +``` + +#### 3.3 Testing using example data #### + +```bash +cd FernandinaSenDT128/mintpy +dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver cpu -o ts_demErr_cpu.h5 +dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5 --solver torch -o ts_demErr_gpu.h5 +``` + +CPU and GPU outputs should agree to float32 round-off (rms / |cpu|.max on the order of 1e-6 for `delta_z`, 1e-8 for `ts_cor`). + +## 4. Behavior notes ## + **VRAM auto-sizing.** `gpuChunkSize = 0` (auto) probes free VRAM at runtime and chooses a per-chunk pixel count with a fixed headroom factor. Set an explicit integer to override (e.g. for reproducible chunking across hosts with different VRAM). @@ -59,10 +98,12 @@ The two outputs should agree to float32 round-off (RMS on the order of 1e-5). + **No silent CPU fallback.** Selecting `solver = torch` on a host without a visible CUDA device raises immediately rather than silently falling back to CPU; this keeps performance regressions visible. -## 4. Performance ## +## 5. Performance ## Benchmarks on RTX 5080 (Blackwell sm_120, CUDA 12.8, PyTorch 2.11) are tracked in the sibling [`mintpy-benchmark`](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark) repository. +### `invert_network` + Tutorial-scale on FernandinaSenDT128 (270k pixels, 288 ifgs): + [report_torch.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/report_torch.md) — `cpu` vs `torch` end-to-end @@ -73,3 +114,14 @@ Tutorial-scale on FernandinaSenDT128 (270k pixels, 288 ifgs): Large-scene on GalapagosSenDT128 (3.4M pixels, 475 kept ifgs; ~12.6× pixels and 1.65× ifgs vs Fernandina): + [report_large_scene.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/report_large_scene.md) — `solver=torch` reaches **36.4× step wall** / **44.4× internal** speedup on `invert_network` (cpu 6189 s → torch 170 s on RTX 5080 / SSD), confirming the speedup grows at scale; output equivalence preserved at float32 round-off (abs RMS max ~16 µm) + +### `correct_topography` + +Five-scene survey (Fernandina + Galapagos + three Zenodo Tier 1 scenes covering ISCE2 / GMTSAR / ARIA / ROI_PAC ingest pipelines, two wavelengths, and D ∈ [24, 333]): + ++ [report_bench_survey.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_bench_survey.md) — combined speedup-vs-K curve, numeric gate validation, axes coverage analysis ++ Per-scene: [report_fernandina.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_fernandina.md), [report_galapagos.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_galapagos.md), [report_sanfranbay.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_sanfranbay.md), [report_sanfran_aria.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_sanfran_aria.md), [report_kuju.md](https://github.com/s-sasaki-earthsea-wizard/mintpy-benchmark/blob/main/reports/dem_error/report_kuju.md) + +Headline: on local SSD, tutorial-scale scenes (~250–325k pixels) typically see **2–2.5×** speedup; production-scale (K = 3.4M Galapagos) reaches **6.15×**. CPU/GPU numeric agreement at `rms / |cpu|.max < 1e-5` on all 5 scenes. Note that the GPU path is **more I/O sensitive** than CPU because its compute portion is sub-second — networked storage (CIFS / NAS) can erase the speedup at small scenes (Fernandina drops from 2.56× on SSD to 0.97× on NAS). + +The closed-form speedup-vs-(K,D,P) model fitted against the survey lives in the wiki: [GPU Speedup Scaling Model](https://github.com/s-sasaki-earthsea-wizard/MintPy/wiki/GPU-Speedup-Scaling-Model). Short version: K is the GPU-parallelisable batch dimension; D and P are per-pixel internal dimensions that scale CPU and GPU walls similarly. Expect ≳ 4× speedup once K ≳ 1M. diff --git a/src/mintpy/cli/dem_error.py b/src/mintpy/cli/dem_error.py index 637df236e..e8fbe7bfa 100755 --- a/src/mintpy/cli/dem_error.py +++ b/src/mintpy/cli/dem_error.py @@ -73,6 +73,21 @@ def create_parser(subparsers=None): '1) output time-series file already exists, readable '+ 'and newer than input interferograms file\n' + '2) all configuration parameters are the same.') + + # solver (pixelwise-geometry branch only; mean-geometry stays CPU) + solver = parser.add_argument_group('solver', 'per-pixel solver for the pixelwise-geometry branch') + solver.add_argument('--solver', dest='solver', default='cpu', + choices={'cpu', 'torch'}, + help='per-pixel solver: cpu (scipy.linalg.lstsq, default) ' + 'or torch (CUDA-batched normal-equation + Cholesky via ' + 'PyTorch). torch requires the [gpu] extras and a visible ' + 'CUDA device; absence is a hard error. Only the ' + 'pixelwise-geometry branch is affected; the ' + 'mean-geometry branch is already pixel-batched on CPU.') + solver.add_argument('--gpu-chunk-size', dest='gpuChunkSize', type=int, default=0, + help='pixels per GPU chunk for --solver=torch ' + '(0=auto-size from free VRAM; default).') + # computing parser = arg_utils.add_memory_argument(parser) parser = arg_utils.add_parallel_argument(parser) @@ -145,6 +160,10 @@ def read_template2inps(template_file, inps): iDict[key] = int(value) elif key in ['excludeDate','stepDate']: iDict[key] = ptime.yyyymmdd(value.split(',')) + elif key in ['solver']: + iDict[key] = value + elif key in ['gpuChunkSize']: + iDict[key] = int(value) # computing configurations dask_key_prefix = 'mintpy.compute.' diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg index 02b4b65da..cb6615828 100644 --- a/src/mintpy/defaults/smallbaselineApp.cfg +++ b/src/mintpy/defaults/smallbaselineApp.cfg @@ -283,6 +283,16 @@ mintpy.topographicResidual.phaseVelocity = auto #[yes / no], auto for no - mintpy.topographicResidual.stepDate = auto #[20080529,20190704T1733 / no], auto for no, date of step jump mintpy.topographicResidual.excludeDate = auto #[20070321 / txtFile / no], auto for exclude_date.txt mintpy.topographicResidual.pixelwiseGeometry = auto #[yes / no], auto for yes, use pixel-wise geometry info +## Per-pixel solver for the pixelwise-geometry branch (GPU path is opt-in): +## a. cpu - scipy.linalg.lstsq, per-pixel (default, original behavior) +## b. torch - batched normal-equation + Cholesky on CUDA via PyTorch. +## Requires the [gpu] extras (CUDA-enabled torch build) and a +## visible CUDA device; absence is a hard error (no silent CPU +## fallback). Only the pixelwise-geometry branch is affected; +## the mean-geometry branch is already pixel-batched on CPU. +## See docs/installation.md and docs/gpu.md. +mintpy.topographicResidual.solver = auto #[cpu / torch], auto for cpu +mintpy.topographicResidual.gpuChunkSize = auto #[int >= 0], auto for 0 (auto-size from free VRAM) ########## 11.1 residual_RMS (root mean squares for noise evaluation) diff --git a/src/mintpy/defaults/smallbaselineApp_auto.cfg b/src/mintpy/defaults/smallbaselineApp_auto.cfg index 654b9a226..796e13afa 100644 --- a/src/mintpy/defaults/smallbaselineApp_auto.cfg +++ b/src/mintpy/defaults/smallbaselineApp_auto.cfg @@ -122,6 +122,8 @@ mintpy.topographicResidual.phaseVelocity = no mintpy.topographicResidual.stepDate = no mintpy.topographicResidual.excludeDate = exclude_date.txt mintpy.topographicResidual.pixelwiseGeometry = yes +mintpy.topographicResidual.solver = cpu +mintpy.topographicResidual.gpuChunkSize = 0 ########## residual_RMS diff --git a/src/mintpy/dem_error.py b/src/mintpy/dem_error.py index d415f53c3..bcfc9a7c1 100644 --- a/src/mintpy/dem_error.py +++ b/src/mintpy/dem_error.py @@ -281,7 +281,8 @@ def estimate_dem_error(ts0, G0, tbase, date_flag=None, phase_velocity=False, def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, - date_flag=None, phase_velocity=False): + date_flag=None, phase_velocity=False, + solver='cpu', chunk_size=None): """ Correct one path of a time-series for DEM error. @@ -291,6 +292,13 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, box - tuple of 4 int in (x0, y0, x1, y1) for the area of interest date_flag - 1D np.ndarray in bool in size of (num_date), dates used for the estimation phase_velocity - bool, minimize the resdiual phase or phase velocity + solver - str, 'cpu' (default, scipy.linalg.lstsq per pixel) or + 'torch' (batched Cholesky on CUDA via mintpy.gpu.dem_error). + Only affects the pixelwise-geometry branch; the + mean-geometry branch is already pixel-batched on CPU. + chunk_size - int or None. Pixels per GPU chunk for solver='torch'. + None / <=0 selects auto-sizing from free VRAM. + Ignored for solver='cpu'. Returns: delta_z - 2D np.ndarray in size of (num_row, num_col) ts_cor - 3D np.ndarray in size of (num_date, num_row, num_col) ts_res - 3D np.ndarray in size of (num_date, num_row, num_col) @@ -386,6 +394,40 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, ts_cor[:, mask] = ts_cor_i ts_res[:, mask] = ts_res_i + elif solver != 'cpu': + print(f'estimating DEM error pixel-wisely via solver={solver!r} ...') + from mintpy.gpu.dem_error import ( + estimate_dem_error_pixelwise_batch, + is_solver_available, + ) + if not is_solver_available(solver): + raise RuntimeError( + f"solver={solver!r} requested but not available " + "(missing torch / CUDA)" + ) + + ts_valid = ts_data[:, mask] + if pbase.shape[1] == 1: + pbase_valid = pbase + else: + pbase_valid = pbase[:, mask] + + delta_z_valid, ts_cor_valid, ts_res_valid = estimate_dem_error_pixelwise_batch( + ts_valid=ts_valid, + pbase_valid=pbase_valid, + range_dist_valid=range_dist[mask], + sin_inc_angle_valid=sin_inc_angle[mask], + G_defo=G_defo, + tbase=tbase, + date_flag=date_flag, + phase_velocity=phase_velocity, + solver=solver, + chunk_size=chunk_size, + ) + delta_z[mask] = delta_z_valid + ts_cor[:, mask] = ts_cor_valid + ts_res[:, mask] = ts_res_valid + else: print('estimating DEM error pixel-wisely ...') prog_bar = ptime.progressBar(maxValue=num_pixel2inv) @@ -499,12 +541,17 @@ def correct_dem_error(inps): ) # 3.2 prepare the input arguments for *_patch() + # chunk_size <=0 means "auto from free VRAM" inside the GPU dispatch; + # the patch function passes it through to mintpy.gpu.dem_error. + chunk_size = getattr(inps, 'gpuChunkSize', 0) or 0 data_kwargs = { 'G_defo' : G_defo, 'ts_file' : inps.ts_file, 'geom_file' : inps.geom_file, 'date_flag' : date_flag, 'phase_velocity' : inps.phaseVelocity, + 'solver' : getattr(inps, 'solver', 'cpu'), + 'chunk_size' : chunk_size if chunk_size > 0 else None, } # 3.3 invert / write block-by-block diff --git a/src/mintpy/gpu/_common.py b/src/mintpy/gpu/_common.py new file mode 100644 index 000000000..bd891d619 --- /dev/null +++ b/src/mintpy/gpu/_common.py @@ -0,0 +1,129 @@ +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# Shared helpers for the mintpy.gpu batched solvers # +############################################################ +# Recommend import: +# from mintpy.gpu import _common as gpu_common + + +"""Shared device/solver helpers for the mintpy.gpu batched solvers. + +Both ``mintpy.gpu.ifgram_inversion`` and ``mintpy.gpu.dem_error`` reduce +a per-pixel weighted (or unweighted) least-squares problem to batched +normal equations + cuSolver-batched Cholesky. This module factors out +the pieces that do not depend on the specific solver: + +* PyTorch / CUDA availability probing (``HAS_TORCH``, + ``is_solver_available``, ``get_torch_device``). +* VRAM-aware chunk sizing (``auto_chunk_size``). +* The normal-equations + rank-deficient fallback core + (``solve_normal_equations_batched``). + +Callers are still responsible for assembling the per-pixel design +matrix and right-hand-side (e.g. applying sqrt-weights, broadcasting a +shared geometric column, ...) before passing them in. +""" + + +try: + import torch + HAS_TORCH = True +except ImportError: + HAS_TORCH = False + + +SUPPORTED_SOLVERS = ('cpu', 'torch') + +# default chunk size when caller does not provide one and VRAM probing fails +DEFAULT_CHUNK_SIZE = 20000 + +# safety factor applied to free VRAM when auto-sizing the chunk +VRAM_SAFETY = 0.4 + + +def is_solver_available(solver): + """Return True if the named solver is importable and usable.""" + if solver == 'cpu': + return True + if solver == 'torch': + return HAS_TORCH and torch.cuda.is_available() + return False + + +def get_torch_device(solver): + """Return ``torch.device('cuda')`` for a GPU solver, else raise.""" + if not HAS_TORCH: + raise ImportError( + f"solver='{solver}' requires PyTorch. " + "Install with `pip install -e \".[gpu]\" " + "--extra-index-url https://download.pytorch.org/whl/cu128 " + "--index-strategy unsafe-best-match`." + ) + if not torch.cuda.is_available(): + raise RuntimeError( + f"solver='{solver}' requires CUDA, but torch.cuda.is_available() is False." + ) + return torch.device('cuda') + + +def auto_chunk_size(num_rows, num_cols, dtype_bytes=4): + """Pick chunk size from free VRAM. + + The dominant per-pixel allocation in both batched solvers is the + design matrix tensor of shape ``(n, num_rows, num_cols)``; ~3x of + that empirically covers temporaries (normal matrix, Cholesky + factor, residual). Returns ``DEFAULT_CHUNK_SIZE`` when CUDA is + unavailable. + """ + if not (HAS_TORCH and torch.cuda.is_available()): + return DEFAULT_CHUNK_SIZE + free_b, _ = torch.cuda.mem_get_info() + per_pixel_bytes = 3 * num_rows * num_cols * dtype_bytes + n = int(VRAM_SAFETY * free_b / max(per_pixel_bytes, 1)) + return max(1, n) + + +def solve_normal_equations_batched(G_batch, y_batch, print_msg=True): + """Per-pixel least-squares via normal equations + batched Cholesky. + + For each pixel k with system ``G_k @ x_k = y_k`` (shapes + ``(num_row, num_col)`` and ``(num_row,)``), solve the normal + equations + ``(G_k^T @ G_k) x_k = G_k^T @ y_k`` + via cuSolver-batched Cholesky. Compared to a per-pixel QR path this + collapses ~num_col Householder iterations into a single batched + factorization. + + Rank-deficient pixels are detected via ``cholesky_ex`` info codes; + their factor is replaced with identity and right-hand-side with + zeros so the downstream ``cholesky_solve`` produces an all-zero + solution for those pixels and never propagates NaN/Inf. + + Args: + G_batch: (n, num_row, num_col) per-pixel design matrix. + Callers solving a weighted system must pre-apply + sqrt-weights to ``G_batch`` and ``y_batch``. + y_batch: (n, num_row) per-pixel observation vector. + print_msg: bool. Warn on rank-deficient pixels. + + Returns: + (n, num_col) per-pixel solution. + """ + G_T = G_batch.transpose(-1, -2) + N = G_T @ G_batch + r = G_T @ y_batch.unsqueeze(-1) + + L, info = torch.linalg.cholesky_ex(N) + fail_mask = info != 0 + if fail_mask.any(): + n_fail = int(fail_mask.sum().item()) + if print_msg: + print(f'WARNING: {n_fail} rank-deficient pixel(s) in chunk; ' + 'setting solution to zero') + eye = torch.eye(N.shape[-1], device=L.device, dtype=L.dtype) + L = torch.where(fail_mask.view(-1, 1, 1), eye, L) + r = torch.where(fail_mask.view(-1, 1, 1), torch.zeros_like(r), r) + + X = torch.cholesky_solve(r, L) + return X.squeeze(-1) diff --git a/src/mintpy/gpu/dem_error.py b/src/mintpy/gpu/dem_error.py new file mode 100644 index 000000000..6d446b976 --- /dev/null +++ b/src/mintpy/gpu/dem_error.py @@ -0,0 +1,238 @@ +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# GPU-accelerated DEM error correction # +############################################################ +# Recommend import: +# from mintpy.gpu import dem_error as gpu_demerr + + +"""GPU-batched solver for the pixelwise DEM-error correction. + +This module provides ``estimate_dem_error_pixelwise_batch``, an opt-in +replacement for the per-pixel CPU loop in ``correct_dem_error_patch`` +(``mintpy.dem_error``). All pixels in a chunk are solved on CUDA via +batched normal equations + cuSolver-batched Cholesky; rank-deficient +pixels are detected via ``cholesky_ex`` info codes and zeroed so NaN/Inf +never propagate. + +Compared to the SBAS network inversion solver in +``mintpy.gpu.ifgram_inversion``, the pixelwise DEM-error design matrix +differs per pixel only in its first column (the geometric scaling +``pbase / (range_dist * sin_inc_angle)``). The Phase-1 implementation +materializes the full ``(K, D, P)`` ``G_batch`` for clarity and code +symmetry with the SBAS solver; a Phase-2 optimisation that exploits the +shared structure is left for a follow-up if benchmarks show the naive +form is leaving performance on the table. + +Currently ``phase_velocity=True`` is not implemented and will raise. The +production default is ``phase_velocity=False`` (see +``smallbaselineApp.cfg``: ``mintpy.topographicResidual.phaseVelocity = +no``). +""" + + +import numpy as np + +from mintpy.gpu._common import ( + HAS_TORCH, + SUPPORTED_SOLVERS, + auto_chunk_size, + get_torch_device, + is_solver_available, + solve_normal_equations_batched, +) + +if HAS_TORCH: + import torch + + +def estimate_dem_error_pixelwise_batch( + ts_valid, + pbase_valid, + range_dist_valid, + sin_inc_angle_valid, + G_defo, + tbase=None, + date_flag=None, + phase_velocity=False, + chunk_size=None, + solver='torch', + print_msg=True, +): + """Batched GPU solver for the pixelwise DEM-error correction. + + Mirrors the CPU per-pixel loop in ``correct_dem_error_patch`` over + the ``mask`` of invertable pixels: builds, per pixel ``k``, + ``G0_k = [pbase_k / (R_k * sinθ_k) | G_defo]`` (D, P) + where ``P = 1 + G_defo.shape[1]``, then solves + ``lstsq(G0_k[date_flag], ts_k[date_flag])`` and assembles + ``ts_cor``/``ts_res`` from the full (un-filtered) ``G0_k``. + + Args: + ts_valid: (num_date, num_pixel) float. Time-series for + the pixels to invert (already mask-filtered + by caller). NaN-free. + pbase_valid: (num_date, 1) for scene-mean baseline, or + (num_date, num_pixel) for per-pixel baseline + (column-aligned with ``ts_valid``). + range_dist_valid: (num_pixel,) slant range distance in metres. + sin_inc_angle_valid:(num_pixel,) sin of incidence angle. + G_defo: (num_date, num_param_defo) shared deformation + design matrix (polynomial / step / etc). + tbase: (num_date,) or (num_date, 1) temporal + baseline in years. Only consumed when + ``phase_velocity=True`` (currently + unsupported). Kept for API parity with the + CPU path. + date_flag: (num_date,) bool. Dates used for the fit. + ``None`` means all dates. + phase_velocity: bool. Only ``False`` is currently + implemented; ``True`` raises + ``NotImplementedError``. + chunk_size: int or None. Pixels per GPU chunk. ``None`` + => auto from free VRAM. + solver: str. Only ``'torch'`` is implemented. + print_msg: bool. + + Returns: + delta_z: (num_pixel,) float32. Estimated DEM residual per + pixel (first parameter of X). + ts_cor: (num_date, num_pixel) float32. Topography-corrected + time-series, ``ts - G0[:, 0] * delta_z``. + ts_res: (num_date, num_pixel) float32. Fit residual, + ``ts - G0 @ X``. + """ + if solver != 'torch': + raise ValueError( + f"unsupported solver={solver!r}; choose from {SUPPORTED_SOLVERS}" + ) + if phase_velocity: + raise NotImplementedError( + "phase_velocity=True is not yet supported on the torch solver; " + "fall back to solver='cpu'" + ) + device = get_torch_device(solver) + + ts_valid = np.asarray(ts_valid, dtype=np.float32) + if ts_valid.ndim != 2: + raise ValueError(f'ts_valid must be 2D (num_date, num_pixel); got {ts_valid.shape}') + num_date, num_pixel = ts_valid.shape + + G_defo_np = np.asarray(G_defo, dtype=np.float32) + if G_defo_np.shape[0] != num_date: + raise ValueError( + f'G_defo first dim {G_defo_np.shape[0]} must match num_date={num_date}' + ) + num_param_defo = G_defo_np.shape[1] + num_param = 1 + num_param_defo + + pbase_np = np.asarray(pbase_valid, dtype=np.float32) + if pbase_np.shape == (num_date,): + pbase_np = pbase_np.reshape(-1, 1) + if pbase_np.ndim != 2 or pbase_np.shape[0] != num_date: + raise ValueError( + f'pbase_valid must be (num_date, 1) or (num_date, num_pixel); ' + f'got {pbase_np.shape}' + ) + pbase_per_pixel = pbase_np.shape[1] != 1 + if pbase_per_pixel and pbase_np.shape[1] != num_pixel: + raise ValueError( + f'pbase_valid column count {pbase_np.shape[1]} must be 1 or ' + f'num_pixel={num_pixel}' + ) + + range_dist_np = np.asarray(range_dist_valid, dtype=np.float32).reshape(-1) + sin_inc_angle_np = np.asarray(sin_inc_angle_valid, dtype=np.float32).reshape(-1) + if range_dist_np.shape != (num_pixel,) or sin_inc_angle_np.shape != (num_pixel,): + raise ValueError( + f'range_dist_valid / sin_inc_angle_valid must be 1D of length num_pixel=' + f'{num_pixel}; got {range_dist_np.shape} / {sin_inc_angle_np.shape}' + ) + + if date_flag is None: + date_flag_np = np.ones(num_date, dtype=np.bool_) + else: + date_flag_np = np.asarray(date_flag, dtype=np.bool_).reshape(-1) + if date_flag_np.shape != (num_date,): + raise ValueError( + f'date_flag must be of length num_date={num_date}; ' + f'got {date_flag_np.shape}' + ) + + delta_z = np.zeros(num_pixel, dtype=np.float32) + ts_cor = np.zeros((num_date, num_pixel), dtype=np.float32) + ts_res = np.zeros((num_date, num_pixel), dtype=np.float32) + + if num_pixel == 0: + return delta_z, ts_cor, ts_res + + if chunk_size is None or chunk_size <= 0: + chunk_size = auto_chunk_size(num_date, num_param) + if print_msg: + free_gib = torch.cuda.mem_get_info()[0] / 2**30 + print(f'GPU auto chunk_size = {chunk_size} pixels ' + f'(free VRAM {free_gib:.1f} GiB)') + else: + chunk_size = int(chunk_size) + + num_chunk = (num_pixel + chunk_size - 1) // chunk_size + if print_msg: + print(f'estimating DEM error via {solver} batched LSQ ' + f'in {num_chunk} chunk(s) of up to {chunk_size} pixels ...') + + # move shared design tensors to GPU once + G_defo_dev = torch.as_tensor(G_defo_np, device=device) # (D, P_defo) + date_flag_dev = torch.as_tensor(date_flag_np, dtype=torch.bool, device=device) + pbase_shared_dev = None + if not pbase_per_pixel: + pbase_shared_dev = torch.as_tensor(pbase_np.reshape(-1), device=device) # (D,) + + for ci in range(num_chunk): + c0 = ci * chunk_size + c1 = min(c0 + chunk_size, num_pixel) + n = c1 - c0 + + # per-pixel scaling c_k = 1 / (R_k * sin θ_k) + c_chunk = 1.0 / (range_dist_np[c0:c1] * sin_inc_angle_np[c0:c1]) + c_dev = torch.as_tensor(c_chunk, device=device) # (n,) + + # G_geom: (n, D, 1) + if pbase_per_pixel: + pbase_chunk = pbase_np[:, c0:c1] # (D, n) + pbase_chunk_dev = torch.as_tensor(pbase_chunk, device=device) + G_geom = (pbase_chunk_dev * c_dev.unsqueeze(0)).t().unsqueeze(-1) + else: + G_geom = (pbase_shared_dev.unsqueeze(0) * c_dev.unsqueeze(-1)).unsqueeze(-1) + + # broadcast G_defo across pixels via expand (view, no copy) + G_defo_b = G_defo_dev.unsqueeze(0).expand(n, -1, -1) # (n, D, P_defo) + G0_batch = torch.cat([G_geom, G_defo_b], dim=-1).contiguous() # (n, D, P) + + ts_chunk_dev = torch.as_tensor(ts_valid[:, c0:c1], device=device).t() # (n, D) + + # apply date_flag for the fit; keep full-D for ts_cor / ts_res assembly + G_fit = G0_batch[:, date_flag_dev, :] # (n, D', P) + y_fit = ts_chunk_dev[:, date_flag_dev] # (n, D') + + X = solve_normal_equations_batched(G_fit, y_fit, print_msg=print_msg) # (n, P) + + delta_z_chunk = X[:, 0] # (n,) + G0_col0 = G0_batch[:, :, 0] # (n, D) + ts_cor_chunk = ts_chunk_dev - G0_col0 * delta_z_chunk.unsqueeze(-1) + ts_pred = (G0_batch @ X.unsqueeze(-1)).squeeze(-1) # (n, D) + ts_res_chunk = ts_chunk_dev - ts_pred + + delta_z[c0:c1] = delta_z_chunk.detach().cpu().numpy().astype(np.float32) + ts_cor[:, c0:c1] = ts_cor_chunk.t().detach().cpu().numpy().astype(np.float32) + ts_res[:, c0:c1] = ts_res_chunk.t().detach().cpu().numpy().astype(np.float32) + + del G0_batch, G_fit, y_fit, G_geom, G_defo_b, X, ts_chunk_dev + del ts_cor_chunk, ts_res_chunk, ts_pred, G0_col0 + + if print_msg: + chunk_step = max(1, num_chunk // 5) + if (ci + 1) % chunk_step == 0 or ci == num_chunk - 1: + print(f'chunk {ci + 1} / {num_chunk}') + + return delta_z, ts_cor, ts_res diff --git a/src/mintpy/gpu/ifgram_inversion.py b/src/mintpy/gpu/ifgram_inversion.py index 9342ce0d2..a179ae884 100644 --- a/src/mintpy/gpu/ifgram_inversion.py +++ b/src/mintpy/gpu/ifgram_inversion.py @@ -26,79 +26,26 @@ import numpy as np -try: +from mintpy.gpu._common import ( + HAS_TORCH, + SUPPORTED_SOLVERS, + auto_chunk_size, + get_torch_device, + is_solver_available, + solve_normal_equations_batched, +) + +if HAS_TORCH: import torch - HAS_TORCH = True -except ImportError: - HAS_TORCH = False - - -SUPPORTED_SOLVERS = ('cpu', 'torch') - -# default chunk size when caller does not provide one and VRAM probing fails -DEFAULT_CHUNK_SIZE = 20000 - -# safety factor applied to free VRAM when auto-sizing the chunk -VRAM_SAFETY = 0.4 - - -def is_solver_available(solver): - """Return True if the named WLS solver is importable and usable.""" - if solver == 'cpu': - return True - if solver == 'torch': - return HAS_TORCH and torch.cuda.is_available() - return False - - -def _get_torch_device(solver): - if not HAS_TORCH: - raise ImportError( - f"solver='{solver}' requires PyTorch. " - "Install with `pip install -e \".[gpu]\" " - "--extra-index-url https://download.pytorch.org/whl/cu128 " - "--index-strategy unsafe-best-match`." - ) - if not torch.cuda.is_available(): - raise RuntimeError( - f"solver='{solver}' requires CUDA, but torch.cuda.is_available() is False." - ) - return torch.device('cuda') - - -def _auto_chunk_size(num_pair, num_unknown, dtype_bytes=4): - """Pick chunk size from free VRAM. - - The dominant per-pixel allocations during one chunk are: - Gw (n, num_pair, num_unknown) ~ num_pair * num_unknown * dtype_bytes - N, L (n, num_unknown, num_unknown) ~ 1/num_pair of the above - residual / temp ~ num_pair * dtype_bytes - Empirically, ~3x of Gw covers the rest. Use VRAM_SAFETY as the - headroom factor. - """ - if not (HAS_TORCH and torch.cuda.is_available()): - return DEFAULT_CHUNK_SIZE - free_b, _ = torch.cuda.mem_get_info() - per_pixel_bytes = 3 * num_pair * num_unknown * dtype_bytes - n = int(VRAM_SAFETY * free_b / max(per_pixel_bytes, 1)) - return max(1, n) def _solve_cholesky(G_dev, w_dev, y_dev): """Per-pixel weighted least-squares via normal equations + batched Cholesky. - For each pixel k with weighted system Gw_k @ x_k = yw_k where - Gw_k = diag(w_k) @ G and yw_k = w_k * y_k, solve the normal equations - (Gw_k^T @ Gw_k) x_k = Gw_k^T @ yw_k - via cuSolver-batched Cholesky. Compared to a gels (QR) path this - collapses ~num_unknown Householder iterations per pixel into a single - batched factorization, reducing per-chunk kernel launches from - ~num_pixel * num_unknown to ~5. - - Rank-deficient pixels are detected via cholesky_ex info codes; their - factor is replaced with identity and right-hand-side with zeros so the - downstream cholesky_solve produces an all-zero solution for those - pixels and never propagates NaN/Inf. + Assembles the weighted system + ``Gw_k = diag(w_k) @ G``, ``yw_k = w_k * y_k`` + in one broadcast and dispatches to the shared + ``solve_normal_equations_batched`` core. Args: G_dev: (num_pair, num_unknown) design matrix on GPU. @@ -109,23 +56,8 @@ def _solve_cholesky(G_dev, w_dev, y_dev): (n, num_unknown) per-pixel solution. """ Gw = w_dev.t().unsqueeze(-1) * G_dev.unsqueeze(0) - yw = (w_dev * y_dev).t().unsqueeze(-1) - Gw_T = Gw.transpose(-1, -2) - N = Gw_T @ Gw - r = Gw_T @ yw - - L, info = torch.linalg.cholesky_ex(N) - fail_mask = info != 0 - if fail_mask.any(): - n_fail = int(fail_mask.sum().item()) - print(f'WARNING: {n_fail} rank-deficient pixel(s) in chunk; ' - 'setting solution to zero') - eye = torch.eye(N.shape[-1], device=L.device, dtype=L.dtype) - L = torch.where(fail_mask.view(-1, 1, 1), eye, L) - r = torch.where(fail_mask.view(-1, 1, 1), torch.zeros_like(r), r) - - X = torch.cholesky_solve(r, L) - return X.squeeze(-1) + yw = (w_dev * y_dev).t() + return solve_normal_equations_batched(Gw, yw) def estimate_timeseries_batch( @@ -182,7 +114,7 @@ def estimate_timeseries_batch( raise ValueError( f"unsupported solver={solver!r}; choose from {SUPPORTED_SOLVERS}" ) - device = _get_torch_device(solver) + device = get_torch_device(solver) G = B if min_norm_velocity else A num_pair, num_unknown = G.shape @@ -201,7 +133,7 @@ def estimate_timeseries_batch( # decide chunk size if chunk_size is None or chunk_size <= 0: - chunk_size = _auto_chunk_size(num_pair, num_unknown) + chunk_size = auto_chunk_size(num_pair, num_unknown) if print_msg: free_gib = torch.cuda.mem_get_info()[0] / 2**30 print(f'GPU auto chunk_size = {chunk_size} pixels ' diff --git a/tests/test_gpu_dem_error.py b/tests/test_gpu_dem_error.py new file mode 100644 index 000000000..89b77ec2e --- /dev/null +++ b/tests/test_gpu_dem_error.py @@ -0,0 +1,284 @@ +"""Numerical comparison tests for src/mintpy/gpu/dem_error.py. + +Compare the batched-Cholesky GPU solver against the per-pixel CPU +reference (``scipy.linalg.lstsq`` via +``mintpy.dem_error.estimate_dem_error``) on small synthetic fixtures. +This file is intentionally lightweight: real-data validation on +FernandinaSenDT128 is a follow-up task. We assert finiteness and shape +only; the RMS / max-abs differences between CPU and GPU outputs are +printed for inspection (use ``pytest -s`` to see them). + +Tests are skipped automatically when PyTorch / CUDA are unavailable. +""" + +import numpy as np +import pytest + +torch = pytest.importorskip("torch") + +from mintpy.dem_error import estimate_dem_error +from mintpy.gpu.dem_error import estimate_dem_error_pixelwise_batch + +requires_cuda = pytest.mark.skipif( + not torch.cuda.is_available(), + reason="CUDA-capable GPU required for mintpy.gpu.dem_error tests", +) + + +# --------------------------------------------------------------------------- +# synthetic fixture builders +# --------------------------------------------------------------------------- + +def _build_design_matrices(num_date, num_param_defo, seed): + """Return (G_defo, tbase, pbase_scene) at the scale of a real SBAS run. + + G_defo columns: [vel * t, 0.5 * acc * t^2, ...] truncated at + num_param_defo. tbase is cumulative random year fractions starting + at 0. pbase_scene is a (D, 1) perpendicular-baseline column in + metres in roughly Sentinel-1 range. + """ + rng = np.random.default_rng(seed) + dt = rng.uniform(0.05, 0.2, size=num_date - 1).astype(np.float32) + tbase = np.concatenate([[0.0], np.cumsum(dt)]).astype(np.float32) # (D,) + cols = [] + for k in range(num_param_defo): + # column = t^(k+1) / (k+1)! — polynomial deformation basis + from math import factorial + cols.append((tbase ** (k + 1) / factorial(k + 1)).astype(np.float32)) + G_defo = np.stack(cols, axis=1) # (D, P_defo) + pbase_scene = rng.uniform(-200.0, 200.0, size=(num_date, 1)).astype(np.float32) + return G_defo, tbase.reshape(-1, 1), pbase_scene + + +def _build_geometry(num_pixel, seed): + """Return per-pixel slant range and sin(inc) in roughly Sentinel-1 + swath range (R ~ 700 km, inc ~ 30-45 deg).""" + rng = np.random.default_rng(seed) + range_dist = rng.uniform(7.0e5, 7.5e5, size=num_pixel).astype(np.float32) + inc_deg = rng.uniform(30.0, 45.0, size=num_pixel).astype(np.float32) + sin_inc = np.sin(np.deg2rad(inc_deg)).astype(np.float32) + return range_dist, sin_inc + + +def _synthesize_ts(G_defo, pbase_valid, range_dist, sin_inc_angle, + num_pixel, noise_std, seed): + """Forward-model ts_valid = G0_k @ X_true_k + Gaussian noise.""" + rng = np.random.default_rng(seed) + num_date, num_param_defo = G_defo.shape + num_param = 1 + num_param_defo + + # X_true: column 0 is delta_z (metres of DEM error), rest are + # polynomial deformation coefficients. Scale them so the geometric + # term and deformation term are comparable in radians. + X_true = np.zeros((num_param, num_pixel), dtype=np.float32) + X_true[0, :] = rng.normal(0.0, 5.0, size=num_pixel) # ~5 m DEM error + for p in range(1, num_param): + X_true[p, :] = rng.normal(0.0, 0.01, size=num_pixel) # small defo amps + + ts = np.zeros((num_date, num_pixel), dtype=np.float32) + pbase_per_pixel = pbase_valid.shape[1] != 1 + for k in range(num_pixel): + pbase_k = pbase_valid if not pbase_per_pixel else pbase_valid[:, k:k + 1] + G_geom = pbase_k / (range_dist[k] * sin_inc_angle[k]) # (D, 1) + G0 = np.hstack((G_geom, G_defo)) # (D, P) + ts[:, k] = (G0 @ X_true[:, k]).astype(np.float32) + ts += rng.normal(0.0, noise_std, size=ts.shape).astype(np.float32) + return ts, X_true + + +def _cpu_reference(ts_valid, pbase_valid, range_dist, sin_inc_angle, + G_defo, tbase, date_flag=None): + """Per-pixel CPU reference mirroring correct_dem_error_patch's + pixelwise branch (phase_velocity=False). + """ + num_date, num_pixel = ts_valid.shape + delta_z = np.zeros(num_pixel, dtype=np.float32) + ts_cor = np.zeros((num_date, num_pixel), dtype=np.float32) + ts_res = np.zeros((num_date, num_pixel), dtype=np.float32) + pbase_per_pixel = pbase_valid.shape[1] != 1 + for k in range(num_pixel): + pbase_k = pbase_valid if not pbase_per_pixel else pbase_valid[:, k:k + 1] + G_geom = pbase_k / (range_dist[k] * sin_inc_angle[k]) + G0 = np.hstack((G_geom, G_defo)).astype(np.float32) + dz, tc, tr = estimate_dem_error( + ts0=ts_valid[:, k], G0=G0, tbase=tbase, + date_flag=date_flag, phase_velocity=False, + ) + delta_z[k] = np.asarray(dz).reshape(-1).item() + ts_cor[:, k] = np.asarray(tc).flatten() + ts_res[:, k] = np.asarray(tr).flatten() + return delta_z, ts_cor, ts_res + + +# --------------------------------------------------------------------------- +# assertion / comparison helpers +# --------------------------------------------------------------------------- + +def _assert_finite_and_shaped(delta_z, ts_cor, ts_res, num_date, num_pixel): + assert delta_z.shape == (num_pixel,), \ + f'delta_z shape {delta_z.shape} != ({num_pixel},)' + assert ts_cor.shape == (num_date, num_pixel), \ + f'ts_cor shape {ts_cor.shape} != ({num_date}, {num_pixel})' + assert ts_res.shape == (num_date, num_pixel), \ + f'ts_res shape {ts_res.shape} != ({num_date}, {num_pixel})' + assert np.all(np.isfinite(delta_z)), 'delta_z has non-finite entries' + assert np.all(np.isfinite(ts_cor)), 'ts_cor has non-finite entries' + assert np.all(np.isfinite(ts_res)), 'ts_res has non-finite entries' + + +def _print_diff(name, cpu, gpu): + """Print rms / max-abs diff between CPU and GPU output for ``name``. + + The point is qualitative inspection, not gating; numerical + equivalence is not asserted in this session. + """ + diff = (cpu.astype(np.float64) - gpu.astype(np.float64)) + scale = float(np.abs(cpu).max()) if cpu.size else 0.0 + scale_safe = max(scale, 1e-12) + rms = float(np.sqrt(np.mean(diff ** 2))) + mx = float(np.abs(diff).max()) if diff.size else 0.0 + print( + f' {name:8s}: rms={rms:.3e} max|diff|={mx:.3e} ' + f'|cpu|max={scale:.3e} rms/|cpu|max={rms / scale_safe:.3e}' + ) + + +# --------------------------------------------------------------------------- +# fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture(scope='module') +def synthetic_inputs(): + """Single synthetic fixture shared across tests. + + Sized small (num_pixel=1000) to keep the test fast; FernandinaSenDT128- + scale validation is a follow-up task. + """ + num_date = 98 + num_pixel = 1000 + num_param_defo = 2 # vel + acc polynomial basis + noise_std = 1e-3 + G_defo, tbase, pbase_scene = _build_design_matrices( + num_date=num_date, num_param_defo=num_param_defo, seed=0, + ) + range_dist, sin_inc = _build_geometry(num_pixel=num_pixel, seed=1) + return dict( + num_date=num_date, + num_pixel=num_pixel, + G_defo=G_defo, + tbase=tbase, + pbase_scene=pbase_scene, + range_dist=range_dist, + sin_inc=sin_inc, + noise_std=noise_std, + ) + + +# --------------------------------------------------------------------------- +# tests +# --------------------------------------------------------------------------- + +@requires_cuda +def test_pixelwise_scene_mean_pbase(synthetic_inputs, capsys): + """Pixelwise GPU path with scene-mean baseline (pbase shape (D, 1)).""" + s = synthetic_inputs + ts, _ = _synthesize_ts( + G_defo=s['G_defo'], pbase_valid=s['pbase_scene'], + range_dist=s['range_dist'], sin_inc_angle=s['sin_inc'], + num_pixel=s['num_pixel'], noise_std=s['noise_std'], seed=2, + ) + cpu = _cpu_reference( + ts_valid=ts, pbase_valid=s['pbase_scene'], + range_dist=s['range_dist'], sin_inc_angle=s['sin_inc'], + G_defo=s['G_defo'], tbase=s['tbase'], + ) + gpu = estimate_dem_error_pixelwise_batch( + ts_valid=ts, pbase_valid=s['pbase_scene'], + range_dist_valid=s['range_dist'], sin_inc_angle_valid=s['sin_inc'], + G_defo=s['G_defo'], tbase=s['tbase'], + phase_velocity=False, solver='torch', print_msg=False, + ) + _assert_finite_and_shaped(*gpu, num_date=s['num_date'], num_pixel=s['num_pixel']) + + with capsys.disabled(): + print(f"\n[scene-mean pbase (D, 1), num_date={s['num_date']}, " + f"num_pixel={s['num_pixel']}]") + _print_diff('delta_z', cpu[0], gpu[0]) + _print_diff('ts_cor', cpu[1], gpu[1]) + _print_diff('ts_res', cpu[2], gpu[2]) + + +@requires_cuda +def test_pixelwise_per_pixel_pbase(synthetic_inputs, capsys): + """Pixelwise GPU path with per-pixel baseline (pbase shape (D, K)).""" + s = synthetic_inputs + rng = np.random.default_rng(3) + # Perturb the scene-mean pbase per-pixel to simulate + # geometryRadar.h5's bperp dataset. + pbase_per_pixel = ( + s['pbase_scene'] + + rng.uniform(-20.0, 20.0, size=(s['num_date'], s['num_pixel'])).astype(np.float32) + ) + ts, _ = _synthesize_ts( + G_defo=s['G_defo'], pbase_valid=pbase_per_pixel, + range_dist=s['range_dist'], sin_inc_angle=s['sin_inc'], + num_pixel=s['num_pixel'], noise_std=s['noise_std'], seed=4, + ) + cpu = _cpu_reference( + ts_valid=ts, pbase_valid=pbase_per_pixel, + range_dist=s['range_dist'], sin_inc_angle=s['sin_inc'], + G_defo=s['G_defo'], tbase=s['tbase'], + ) + gpu = estimate_dem_error_pixelwise_batch( + ts_valid=ts, pbase_valid=pbase_per_pixel, + range_dist_valid=s['range_dist'], sin_inc_angle_valid=s['sin_inc'], + G_defo=s['G_defo'], tbase=s['tbase'], + phase_velocity=False, solver='torch', print_msg=False, + ) + _assert_finite_and_shaped(*gpu, num_date=s['num_date'], num_pixel=s['num_pixel']) + + with capsys.disabled(): + print(f"\n[per-pixel pbase (D, K), num_date={s['num_date']}, " + f"num_pixel={s['num_pixel']}]") + _print_diff('delta_z', cpu[0], gpu[0]) + _print_diff('ts_cor', cpu[1], gpu[1]) + _print_diff('ts_res', cpu[2], gpu[2]) + + +@requires_cuda +def test_phase_velocity_true_raises(synthetic_inputs): + """phase_velocity=True is not yet implemented and must raise.""" + s = synthetic_inputs + ts, _ = _synthesize_ts( + G_defo=s['G_defo'], pbase_valid=s['pbase_scene'], + range_dist=s['range_dist'], sin_inc_angle=s['sin_inc'], + num_pixel=8, noise_std=0.0, seed=5, + ) + with pytest.raises(NotImplementedError, match='phase_velocity=True'): + estimate_dem_error_pixelwise_batch( + ts_valid=ts, + pbase_valid=s['pbase_scene'], + range_dist_valid=s['range_dist'][:8], + sin_inc_angle_valid=s['sin_inc'][:8], + G_defo=s['G_defo'], tbase=s['tbase'], + phase_velocity=True, solver='torch', print_msg=False, + ) + + +@requires_cuda +def test_unsupported_solver_raises(synthetic_inputs): + s = synthetic_inputs + ts, _ = _synthesize_ts( + G_defo=s['G_defo'], pbase_valid=s['pbase_scene'], + range_dist=s['range_dist'], sin_inc_angle=s['sin_inc'], + num_pixel=8, noise_std=0.0, seed=6, + ) + with pytest.raises(ValueError, match='unsupported solver'): + estimate_dem_error_pixelwise_batch( + ts_valid=ts, + pbase_valid=s['pbase_scene'], + range_dist_valid=s['range_dist'][:8], + sin_inc_angle_valid=s['sin_inc'][:8], + G_defo=s['G_defo'], tbase=s['tbase'], + phase_velocity=False, solver='cupy', print_msg=False, + )