Refactor compressed_reference_idx, fix relative index bug

src/dolphin/phase_link/_compress.py

@@ -9,6 +9,7 @@
 def compress(
     slc_stack: ArrayLike,
     pl_cpx_phase: ArrayLike,
+    first_real_slc_idx: int = 0,
     slc_mean: ArrayLike | None = None,
     reference_idx: int | None = None,
 ):
@@ -21,6 +22,10 @@ def compress(
     pl_cpx_phase : ArrayLike
         The estimated complex phase from phase linking.
         shape  = (nslc, rows // strides.y, cols // strides.x)
+    first_real_slc_idx : int
+        Index within `slc_stack` where the "real" SLCs start.
+        Indexes before this (which are assumed to be already Compressed SLCs) are
+        excluded from the dot product during compression
     slc_mean : ArrayLike, optional
         The mean SLC magnitude, shape (rows, cols), to use as output pixel magnitudes.
         If None, the mean is computed from the input SLC stack.
@@ -39,18 +44,20 @@ def compress(
         pl_referenced = pl_cpx_phase * pl_cpx_phase[reference_idx][None, :, :].conj()
     else:
         pl_referenced = pl_cpx_phase
+
+    # Slice away the compressed layers *after* doing the reference.
+    pl_referenced = pl_referenced[first_real_slc_idx:, :, :]
+    slcs = slc_stack[first_real_slc_idx:]
     # If the output is downsampled, we need to make `pl_cpx_phase` the same shape
     # as the output
-    pl_estimate_upsampled = upsample_nearest(pl_referenced, slc_stack.shape[1:])
+    pl_estimate_upsampled = upsample_nearest(pl_referenced, slcs.shape[1:])
     # For each pixel, project the SLCs onto the (normalized) estimated phase
     # by performing a pixel-wise complex dot product
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="invalid value encountered")
-        phase = np.angle(
-            np.nansum(slc_stack * np.conjugate(pl_estimate_upsampled), axis=0)
-        )
+        phase = np.angle(np.nansum(slcs * np.conjugate(pl_estimate_upsampled), axis=0))
     if slc_mean is None:
-        slc_mean = np.mean(np.abs(slc_stack), axis=0)
+        slc_mean = np.mean(np.abs(slcs), axis=0)
     # If the phase is invalid, set the mean to NaN
     slc_mean[phase == 0] = np.nan
     return slc_mean * np.exp(1j * phase)

src/dolphin/stack.py

@@ -5,6 +5,7 @@
 import logging
 import warnings
 from datetime import date, datetime
+from enum import Enum
 from os import fspath
 from pathlib import Path
 from typing import Optional, Sequence
@@ -20,6 +21,14 @@
 logger = logging.getLogger(__name__)
 
 
+class CompressedSlcPlan(str, Enum):
+    """Plan for creating Compressed SLCs during phase linking."""
+
+    ALWAYS_FIRST = "always_first"
+    FIRST_PER_MINISTACK = "first_per_ministack"
+    LAST_PER_MINISTACK = "last_per_ministack"
+
+
 class BaseStack(BaseModel):
     """Base class for mini- and full stack classes."""
 
@@ -342,6 +351,11 @@ class MiniStackInfo(BaseStack):
         ),
     )
 
+    def __rich_repr__(self):
+        yield from super().__rich_repr__()
+        yield "output_reference_idx", self.output_reference_idx
+        yield "compressed_reference_idx", self.compressed_reference_idx
+
     @property
     def output_reference_date(self):
         """Date of the reference phase of the stack."""
@@ -391,6 +405,7 @@ class MiniStackPlanner(BaseStack):
         0,
         description="Index of the SLC to use as reference during phase linking",
     )
+    compressed_slc_plan: CompressedSlcPlan = CompressedSlcPlan.ALWAYS_FIRST
 
     def plan(
         self, ministack_size: int, compressed_idx: int | None = None
@@ -448,12 +463,20 @@ def plan(
             )
             cur_output_folder = self.output_folder / new_date_str
 
-            if compressed_idx is None:
-                # TODO: To change to Ansari-style ministack references, change this
-                # so that a new `output_reference_idx` gets made
-                compressed_reference_idx = self.output_reference_idx
-            else:
+            if compressed_idx is not None:
                 compressed_reference_idx = compressed_idx
+            elif self.compressed_slc_plan == CompressedSlcPlan.ALWAYS_FIRST:
+                # Simplest operational version: CompSLCs have same base phase,
+                # but different "residual" added on
+                compressed_reference_idx = 0
+            elif self.compressed_slc_plan == CompressedSlcPlan.FIRST_PER_MINISTACK:
+                # Like Ansari, 2017 paper: each ministack is "self contained"
+                compressed_reference_idx = num_ccslc
+                # Ansari, 2017 also had output_reference_idx = num_ccslcs, and
+                # use the "Datum Adjustment" step to get outputs relative to day 0
+            elif self.compressed_slc_plan == CompressedSlcPlan.LAST_PER_MINISTACK:
+                # Alternative that allows sequential interferograms across ministacks
+                compressed_reference_idx = -1
 
             cur_ministack = MiniStackInfo(
                 file_list=combined_files,
@@ -469,3 +492,9 @@ def plan(
             compressed_slc_infos.append(cur_comp_slc)
 
         return output_ministacks
+
+    def __rich_repr__(self):
+        yield from super().__rich_repr__()
+        yield "max_num_compressed", self.max_num_compressed
+        yield "output_reference_idx", self.output_reference_idx
+        yield "compressed_slc_plan", self.compressed_slc_plan

src/dolphin/workflows/config/_common.py

@@ -19,6 +19,7 @@
 from dolphin import __version__ as _dolphin_version
 from dolphin._types import Bbox
 from dolphin.io import DEFAULT_HDF5_OPTIONS, DEFAULT_TIFF_OPTIONS
+from dolphin.stack import CompressedSlcPlan
 
 from ._enums import ShpMethod
 from ._yaml_model import YamlModel
@@ -129,6 +130,7 @@ class PhaseLinkingOptions(BaseModel, extra="forbid"):
             "`n` interferograms. `baseline_line` must be positive."
         ),
     )
+    compressed_slc_plan: CompressedSlcPlan = CompressedSlcPlan.ALWAYS_FIRST
 
 
 class InterferogramNetwork(BaseModel, extra="forbid"):

src/dolphin/workflows/sequential.py

@@ -19,7 +19,7 @@
 from dolphin._types import Filename
 from dolphin.io import VRTStack
 from dolphin.similarity import create_similarities
-from dolphin.stack import MiniStackPlanner
+from dolphin.stack import CompressedSlcPlan, MiniStackPlanner
 
 from .config import ShpMethod
 from .single import run_wrapped_phase_single
@@ -45,6 +45,7 @@ def run_wrapped_phase_sequential(
     shp_nslc: Optional[int] = None,
     use_evd: bool = False,
     beta: float = 0.00,
+    compressed_slc_plan: CompressedSlcPlan = CompressedSlcPlan.ALWAYS_FIRST,
     max_num_compressed: int = 100,
     output_reference_idx: int = 0,
     new_compressed_reference_idx: int | None = None,
@@ -69,6 +70,7 @@ def run_wrapped_phase_sequential(
         output_folder=output_folder,
         max_num_compressed=max_num_compressed,
         output_reference_idx=output_reference_idx,
+        compressed_slc_plan=compressed_slc_plan,
     )
     ministacks = ministack_planner.plan(
         ministack_size, compressed_idx=new_compressed_reference_idx

src/dolphin/workflows/single.py

@@ -237,18 +237,14 @@ def run_wrapped_phase_single(
         # Get the mean to set as pixel magnitudes
         abs_stack = np.abs(cur_data[first_real_slc_idx:, in_trim_rows, in_trim_cols])
         cur_data_mean, cur_amp_dispersion, _ = calc_ps_block(abs_stack)
-        # NOTE: the `ministack.compressed_reference_idx` was set relative to *all*
-        # input images, real and compressed.
-        # Since we pass in just the first real idx, we have to subtract that off
-        ref_index_relative = ministack.compressed_reference_idx - first_real_slc_idx
         cur_comp_slc = compress(
             # Get the inner portion of the full-res SLC data
-            cur_data[first_real_slc_idx:, in_trim_rows, in_trim_cols],
-            pl_output.cpx_phase[first_real_slc_idx:, out_trim_rows, out_trim_cols],
+            cur_data[:, in_trim_rows, in_trim_cols],
+            pl_output.cpx_phase[:, out_trim_rows, out_trim_cols],
+            first_real_slc_idx=first_real_slc_idx,
             slc_mean=cur_data_mean,
-            reference_idx=ref_index_relative,
+            reference_idx=ministack.compressed_reference_idx,
         )
-        # TODO: truncate
 
         # ### Save results ###
 

tests/test_show_versions.py

@@ -15,7 +15,6 @@ def test_get_deps_info():
     assert "osgeo.gdal" in deps_info
     assert "numpy" in deps_info
     assert "pydantic" in deps_info
-    assert "ruamel_yaml" in deps_info
     assert "h5py" in deps_info