Repository: isce-framework/dolphin (language: python)

Issue:
Partially empty blocks produce SHP count of 1 instead of 0
Attempting to stitch SHP rasters per burst is failing
![image](https://github.com/isce-framework/dolphin/assets/8291800/d54a1d85-3643-4207-8145-67e00d7e0681)
Hovering over the dark pixels near the valid data shows `1` instead of `0`. This is likely from some coherence matrix or SHP counting which outputs 1, even when there's no data, for blocks which are partially covered by data
Partially empty blocks produce SHP count of 1 instead of 0
Attempting to stitch SHP rasters per burst is failing
![image](https://github.com/isce-framework/dolphin/assets/8291800/d54a1d85-3643-4207-8145-67e00d7e0681)
Hovering over the dark pixels near the valid data shows `1` instead of `0`. This is likely from some coherence matrix or SHP counting which outputs 1, even when there's no data, for blocks which are partially covered by data


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6c2c6abf..06229ad4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,24 @@
-# [Unreleased](https://github.com/isce-framework/dolphin/compare/v0.17.0...main)
+# [Unreleased](https://github.com/isce-framework/dolphin/compare/v0.18.0...main)
+
+**Removed**
+- the KL-divergence SHP estimator has been removed. GLRT is recommended instead.
+
+**Fixed**
+- `reproject_bounds` uses the `rasterio` version, which densifies points along edges for a more accurate bounding box
+- The output SHP rasters now output 0 if there was no valid input data
+
+**Changed**
+- Unweighted time series inversion will make one batch call, providing a large speedup over the `vmap` version for weighted least squares
+
+# [v0.18.0](https://github.com/isce-framework/dolphin/compare/v0.17.0...v0.18.0)
+
+**Added**
+- `dolphin timeseries` command line tool for inverting unwrapped interferogram network and estimating velocity
+
+**Fixed**
+- Parse the file names correctly to find compressed SLCs and read dates based on production file naming convention
+
+
 # [v0.17.0](https://github.com/isce-framework/dolphin/compare/v0.16.3...v0.17.0) - 2024-04-10)
 **Added**
 - Added Goldstein filtering for unwrapping
diff --git a/src/dolphin/io/_writers.py b/src/dolphin/io/_writers.py
index 71c2d3fe..10eb05c1 100644
--- a/src/dolphin/io/_writers.py
+++ b/src/dolphin/io/_writers.py
@@ -41,10 +41,10 @@ class BackgroundBlockWriter(BackgroundWriter):
     """Class to write data to multiple files in the background using `gdal` bindings."""
 
     def __init__(self, *, max_queue: int = 0, debug: bool = False, **kwargs):
-        if debug is False:
-            super().__init__(nq=max_queue, name="Writer", **kwargs)
-        else:
-            # Don't start a background thread. Just synchronously write data
+        super().__init__(nq=max_queue, name="Writer")
+        if debug:
+            #  background thread. Just synchronously write data
+            self.notify_finished()
             self.queue_write = self.write  # type: ignore[assignment]
 
     def write(
@@ -326,11 +326,12 @@ class BackgroundRasterWriter(BackgroundWriter, DatasetWriter):
     def __init__(
         self, filename: Filename, *, max_queue: int = 0, debug: bool = False, **kwargs
     ):
-        if debug is False:
-            super().__init__(nq=max_queue, name="Writer")
-        else:
-            # Don't start a background thread. Just synchronously write data
+        super().__init__(nq=max_queue, name="Writer")
+        if debug:
+            #  background thread. Just synchronously write data
+            self.notify_finished()
             self.queue_write = self.write  # type: ignore[assignment]
+
         if Path(filename).exists():
             self._raster = RasterWriter(filename)
         else:
@@ -397,10 +398,10 @@ def __init__(
     ):
         from dolphin.io import write_arr
 
-        if debug is False:
-            super().__init__(nq=max_queue, name="GdalStackWriter")
-        else:
-            # Don't start a background thread. Just synchronously write data
+        super().__init__(nq=max_queue, name="GdalStackWriter")
+        if debug:
+            # Stop background thread. Just synchronously write data
+            self.notify_finished()
             self.queue_write = self.write  # type: ignore[assignment]
 
         for fn in file_list:
diff --git a/src/dolphin/phase_link/_core.py b/src/dolphin/phase_link/_core.py
index 27500f1e..bdcf61f7 100644
--- a/src/dolphin/phase_link/_core.py
+++ b/src/dolphin/phase_link/_core.py
@@ -45,6 +45,9 @@ class PhaseLinkOutput(NamedTuple):
     A goodness of fit parameter from 0 to 1 at each pixel.
     """
 
+    shp_counts: np.ndarray
+    """Number of neighbor pixels used in adaptive multilooking."""
+
     eigenvalues: np.ndarray
     """The smallest (largest) eigenvalue resulting from EMI (EVD)."""
 
@@ -210,12 +213,13 @@ def run_phase_linking(
     temp_coh[mask_looked] = np.nan
 
     return PhaseLinkOutput(
-        cpx_phase,
-        temp_coh,
+        cpx_phase=cpx_phase,
+        temp_coh=temp_coh,
+        shp_counts=np.asarray(cpl_out.shp_counts),
         # Convert the rest to numpy for writing
-        np.array(cpl_out.eigenvalues),
-        np.array(cpl_out.estimator),
-        cpl_out.avg_coh,
+        eigenvalues=np.asarray(cpl_out.eigenvalues),
+        estimator=np.asarray(cpl_out.estimator),
+        avg_coh=cpl_out.avg_coh,
     )
 
 
@@ -300,6 +304,15 @@ def run_cpl(
     # Get the temporal coherence
     temp_coh = metrics.estimate_temp_coh(cpx_phase, C_arrays)
 
+    # Reshape the (rows, cols, nslcs) output to be same as input stack
+    cpx_phase_reshaped = jnp.moveaxis(cpx_phase, -1, 0)
+
+    # Get the SHP counts for each pixel (if not using Rect window)
+    if neighbor_arrays is None:
+        shp_counts = jnp.zeros(temp_coh.shape, dtype=np.int16)
+    else:
+        shp_counts = jnp.sum(neighbor_arrays, axis=(-2, -1))
+
     if calc_average_coh:
         # If requested, average the Cov matrix at each row for reference selection
         avg_coh_per_date = jnp.abs(C_arrays).mean(axis=3)
@@ -307,10 +320,13 @@ def run_cpl(
     else:
         avg_coh = None
 
-    # Reshape the (rows, cols, nslcs) output to be same as input stack
-    cpx_phase_reshaped = jnp.moveaxis(cpx_phase, -1, 0)
     return PhaseLinkOutput(
-        cpx_phase_reshaped, temp_coh, eigenvalues, estimator, avg_coh
+        cpx_phase=cpx_phase_reshaped,
+        temp_coh=temp_coh,
+        shp_counts=shp_counts,
+        eigenvalues=eigenvalues,
+        estimator=estimator,
+        avg_coh=avg_coh,
     )
 
 
diff --git a/src/dolphin/shp/__init__.py b/src/dolphin/shp/__init__.py
index 600e3c2e..86cd4107 100644
--- a/src/dolphin/shp/__init__.py
+++ b/src/dolphin/shp/__init__.py
@@ -8,7 +8,7 @@
 from dolphin._log import get_log
 from dolphin.workflows import ShpMethod
 
-from . import _glrt, _kld, _ks
+from . import _glrt, _ks
 
 logger = get_log(__name__)
 
@@ -66,7 +66,7 @@ def estimate_neighbors(
     Raises
     ------
     ValueError
-        - nslc is not provided for GLRT/KLD methods
+        - nslc is not provided for GLRT method
         - amp_stack is not provided for the KS method.
         - `method` not a valid `ShpMethod`
 
@@ -76,11 +76,6 @@ def estimate_neighbors(
     if strides is None:
         strides = {"x": 1, "y": 1}
     logger.debug(f"NUMBA THREADS: {numba.get_num_threads()}")
-    if method.lower() in (ShpMethod.KLD, ShpMethod.GLRT):
-        if mean is None:
-            mean = np.mean(amp_stack, axis=0)
-        if var is None:
-            var = np.var(amp_stack, axis=0)
 
     if method == ShpMethod.RECT:
         # No estimation needed
@@ -90,6 +85,10 @@ def estimate_neighbors(
         if nslc is None:
             msg = "`nslc` must be provided for GLRT method"
             raise ValueError(msg)
+        if mean is None:
+            mean = np.mean(amp_stack, axis=0)
+        if var is None:
+            var = np.var(amp_stack, axis=0)
         neighbor_arrays = _glrt.estimate_neighbors(
             mean=mean,
             var=var,
@@ -99,20 +98,6 @@ def estimate_neighbors(
             alpha=alpha,
             prune_disconnected=prune_disconnected,
         )
-    elif method.lower() == ShpMethod.KLD:
-        logger.debug("Estimating SHP neighbors using KLD")
-        if nslc is None:
-            msg = "`nslc` must be provided for GLRT method"
-            raise ValueError(msg)
-        neighbor_arrays = _kld.estimate_neighbors(
-            mean=mean,
-            var=var,
-            halfwin_rowcol=halfwin_rowcol,
-            strides=strides,
-            nslc=nslc,
-            alpha=alpha,
-            prune_disconnected=prune_disconnected,
-        )
     elif method.lower() == ShpMethod.KS:
         if amp_stack is None:
             msg = "amp_stack must be provided for KS method"
diff --git a/src/dolphin/shp/_common.py b/src/dolphin/shp/_common.py
index e5e374fe..4951865b 100644
--- a/src/dolphin/shp/_common.py
+++ b/src/dolphin/shp/_common.py
@@ -1,18 +1,9 @@
 from __future__ import annotations
 
-import csv
-from functools import lru_cache
-from pathlib import Path
-from typing import Callable
-
 import numba
 import numpy as np
 from numpy.typing import ArrayLike
 
-from dolphin.utils import _get_slices
-
-_get_slices = numba.njit(_get_slices)
-
 
 @numba.njit(nogil=True)
 def remove_unconnected(data: ArrayLike, inplace: bool = True) -> np.ndarray:
@@ -86,90 +77,3 @@ def remove_unconnected(data: ArrayLike, inplace: bool = True) -> np.ndarray:
     # Finally only keep ones from the original that were visited
     data &= visited
     return data
-
-
-# Factory function to make a JIT-ed function calling `compute_test_stat`
-# https://numba.readthedocs.io/en/stable/user/faq.html#can-i-pass-a-function-as-an-argument-to-a-jitted-function
-
-
-def _make_loop_function(
-    compute_test_stat: Callable,
-):
-    """Create a JIT-ed function computing a test statistic for each pixel.
-
-    `compute_test_stat` should have the signature:
-
-    compute_test_stat(rayleigh_scale1: float, rayleigh_scale2: float) -> float
-    """
-
-    @numba.njit(nogil=True, parallel=True)
-    def _loop_over_pixels(
-        mean: ArrayLike,
-        var: ArrayLike,
-        halfwin_rowcol: tuple[int, int],
-        strides_rowcol: tuple[int, int],
-        threshold: float,
-        prune_disconnected: bool,
-        is_shp: np.ndarray,
-    ) -> np.ndarray:
-        """Loop common to SHP tests using only mean and variance."""
-        half_row, half_col = halfwin_rowcol
-        row_strides, col_strides = strides_rowcol
-        # location to start counting from in the larger input
-        r0, c0 = row_strides // 2, col_strides // 2
-        in_rows, in_cols = mean.shape
-        out_rows, out_cols = is_shp.shape[:2]
-
-        # Convert mean/var to the Rayleigh scale parameter
-        scale_squared = (var + mean**2) / 2
-
-        for out_r in numba.prange(out_rows):
-            for out_c in range(out_cols):
-                in_r = r0 + out_r * row_strides
-                in_c = c0 + out_c * col_strides
-
-                scale_1 = scale_squared[in_r, in_c]
-                # Clamp the window to the image bounds
-                (r_start, r_end), (c_start, c_end) = _get_slices(
-                    half_row, half_col, in_r, in_c, in_rows, in_cols
-                )
-
-                for in_r2 in range(r_start, r_end):
-                    for in_c2 in range(c_start, c_end):
-                        # window offsets for dims 3,4 of `is_shp`
-                        r_off = in_r2 - r_start
-                        c_off = in_c2 - c_start
-
-                        # itself is always a neighbor
-                        if in_r2 == in_r and in_c2 == in_c:
-                            is_shp[out_r, out_c, r_off, c_off] = True
-                            continue
-                        scale_2 = scale_squared[in_r2, in_c2]
-
-                        T = compute_test_stat(scale_1, scale_2)
-
-                        is_shp[out_r, out_c, r_off, c_off] = threshold > T
-                if prune_disconnected:
-                    # For this pixel, prune the groups not connected to the center
-                    remove_unconnected(is_shp[out_r, out_c], inplace=True)
-
-        return is_shp
-
-    return _loop_over_pixels
-
-
-@lru_cache
-def _read_cutoff_csv(test_name):
-    # Replace with the actual filename
-    filename = Path(__file__).parent / f"{test_name}_cutoffs.csv"
-
-    result = {}
-    with open(filename) as file:
-        reader = csv.DictReader(file)
-        for row in reader:
-            n = int(row["N"])
-            alpha = float(row["alpha"])
-            cutoff = float(row["cutoff"])
-            result[(n, alpha)] = cutoff
-
-    return result
diff --git a/src/dolphin/shp/_glrt.py b/src/dolphin/shp/_glrt.py
index 366f935d..13dc6d6f 100644
--- a/src/dolphin/shp/_glrt.py
+++ b/src/dolphin/shp/_glrt.py
@@ -1,6 +1,9 @@
 from __future__ import annotations
 
+import csv
+from functools import lru_cache
 from math import log
+from pathlib import Path
 from typing import Optional
 
 import numba
@@ -8,19 +11,11 @@
 from numpy.typing import ArrayLike
 
 from dolphin._types import Strides
-from dolphin.utils import compute_out_shape
+from dolphin.utils import _get_slices, compute_out_shape
 
-from ._common import _make_loop_function, _read_cutoff_csv
+from ._common import remove_unconnected
 
-
-@numba.njit(nogil=True)
-def _compute_glrt_test_stat(scale_1, scale_2):
-    """Compute the GLRT test statistic."""
-    scale_pooled = (scale_1 + scale_2) / 2
-    return 2 * log(scale_pooled) - log(scale_1) - log(scale_2)
-
-
-_loop_over_pixels = _make_loop_function(_compute_glrt_test_stat)
+_get_slices = numba.njit(_get_slices)
 
 
 def estimate_neighbors(
@@ -123,9 +118,89 @@ def get_cutoff(alpha: float, N: int) -> float:
         Cutoff value for the GLRT test statistic.
 
     """
-    n_alpha_to_cutoff = _read_cutoff_csv("glrt")
+    n_alpha_to_cutoff = _read_cutoff_csv()
     try:
         return n_alpha_to_cutoff[(N, alpha)]
     except KeyError as e:
         msg = f"Not implemented for {N = }, {alpha = }"
         raise NotImplementedError(msg) from e
+
+
+@numba.njit(nogil=True)
+def _compute_glrt_test_stat(scale_1, scale_2):
+    """Compute the GLRT test statistic."""
+    scale_pooled = (scale_1 + scale_2) / 2
+    return 2 * log(scale_pooled) - log(scale_1) - log(scale_2)
+
+
+@numba.njit(nogil=True, parallel=True)
+def _loop_over_pixels(
+    mean: ArrayLike,
+    var: ArrayLike,
+    halfwin_rowcol: tuple[int, int],
+    strides_rowcol: tuple[int, int],
+    threshold: float,
+    prune_disconnected: bool,
+    is_shp: np.ndarray,
+) -> np.ndarray:
+    """Loop common to SHP tests using only mean and variance."""
+    half_row, half_col = halfwin_rowcol
+    row_strides, col_strides = strides_rowcol
+    # location to start counting from in the larger input
+    r0, c0 = row_strides // 2, col_strides // 2
+    in_rows, in_cols = mean.shape
+    out_rows, out_cols = is_shp.shape[:2]
+
+    # Convert mean/var to the Rayleigh scale parameter
+    scale_squared = (var + mean**2) / 2
+
+    for out_r in numba.prange(out_rows):
+        for out_c in range(out_cols):
+            in_r = r0 + out_r * row_strides
+            in_c = c0 + out_c * col_strides
+
+            scale_1 = scale_squared[in_r, in_c]
+            # Clamp the window to the image bounds
+            (r_start, r_end), (c_start, c_end) = _get_slices(
+                half_row, half_col, in_r, in_c, in_rows, in_cols
+            )
+            if mean[in_r, in_c] == 0:
+                # Skip nodata pixels
+                continue
+
+            for in_r2 in range(r_start, r_end):
+                for in_c2 in range(c_start, c_end):
+                    # window offsets for dims 3,4 of `is_shp`
+                    r_off = in_r2 - r_start
+                    c_off = in_c2 - c_start
+
+                    # Don't count itself as a neighbor
+                    if in_r2 == in_r and in_c2 == in_c:
+                        is_shp[out_r, out_c, r_off, c_off] = False
+                        continue
+                    scale_2 = scale_squared[in_r2, in_c2]
+
+                    T = _compute_glrt_test_stat(scale_1, scale_2)
+
+                    is_shp[out_r, out_c, r_off, c_off] = threshold > T
+            if prune_disconnected:
+                # For this pixel, prune the groups not connected to the center
+                remove_unconnected(is_shp[out_r, out_c], inplace=True)
+
+    return is_shp
+
+
+@lru_cache
+def _read_cutoff_csv():
+    filename = Path(__file__).parent / "glrt_cutoffs.csv"
+
+    result = {}
+    with open(filename) as file:
+        reader = csv.DictReader(file)
+        for row in reader:
+            n = int(row["N"])
+            alpha = float(row["alpha"])
+            cutoff = float(row["cutoff"])
+            result[(n, alpha)] = cutoff
+
+    return result
diff --git a/src/dolphin/shp/_kld.py b/src/dolphin/shp/_kld.py
deleted file mode 100644
index 9355ff14..00000000
--- a/src/dolphin/shp/_kld.py
+++ /dev/null
@@ -1,131 +0,0 @@
-from __future__ import annotations
-
-from typing import Optional
-
-import numba
-import numpy as np
-from numpy.typing import ArrayLike
-
-from dolphin._types import Strides
-from dolphin.utils import compute_out_shape
-
-from ._common import _make_loop_function, _read_cutoff_csv
-
-
-@numba.njit(nogil=True)
-def _kld(sigma2_p, sigma2_q):
-    """Compute KL divergence of two Rayleigh PDFs given their scale parameters."""
-    return np.log(sigma2_q / sigma2_p) + (sigma2_p / sigma2_q) - 1
-
-
-@numba.njit(nogil=True)
-def _compute_test_stat_js(sigma2_p, sigma2_q):
-    """Compute the Jensen-Shannon divergence."""
-    return (_kld(sigma2_p, sigma2_q) + _kld(sigma2_q, sigma2_p)) / 2
-
-
-_loop_over_pixels = _make_loop_function(_compute_test_stat_js)
-
-
-def estimate_neighbors(
-    mean: ArrayLike,
-    var: ArrayLike,
-    halfwin_rowcol: tuple[int, int],
-    nslc: int,
-    strides: Optional[dict] = None,
-    alpha: float = 0.05,
-    prune_disconnected: bool = False,
-):
-    """Estimate the number of neighbors using the Jensen-Shannon.
-
-    Parameters
-    ----------
-    mean : ArrayLike, 2D
-        Mean amplitude of each pixel.
-    var: ArrayLike, 2D
-        Variance of each pixel's amplitude.
-    halfwin_rowcol : tuple[int, int]
-        Half the size of the block in (row, col) dimensions
-    nslc : int
-        Number of images in the stack used to compute `mean` and `var`.
-        Used to compute the degrees of freedom for the t- and F-tests to
-        determine the critical values.
-    strides: dict, optional
-        The (x, y) strides (in pixels) to use for the sliding window.
-        By default {"x": 1, "y": 1}
-    alpha : float, default=0.05
-        Significance level at which to reject the null hypothesis.
-        Rejecting means declaring a neighbor is not a SHP.
-    prune_disconnected : bool, default=False
-        If True, keeps only SHPs that are 8-connected to the current pixel.
-        Otherwise, any pixel within the window may be considered an SHP, even
-        if it is not directly connected.
-
-
-    Notes
-    -----
-    When `strides` is not (1, 1), the output first two dimensions
-    are smaller than `mean` and `var` by a factor of `strides`. This
-    will match the downstream shape of the strided phase linking results.
-
-    Returns
-    -------
-    is_shp : np.ndarray, 4D
-        Boolean array marking which neighbors are SHPs for each pixel in the block.
-        Shape is (out_rows, out_cols, window_rows, window_cols), where
-            `out_rows` and `out_cols` are computed by
-            `[dolphin.io.compute_out_shape][]`
-            `window_rows = 2 * halfwin_rowcol[0] + 1`
-            `window_cols = 2 * halfwin_rowcol[1] + 1`
-
-    """
-    if strides is None:
-        strides = {"x": 1, "y": 1}
-    half_row, half_col = halfwin_rowcol
-    rows, cols = mean.shape
-
-    threshold = get_cutoff(nslc, alpha)
-
-    strides_rowcol = (strides["y"], strides["x"])
-    out_rows, out_cols = compute_out_shape((rows, cols), Strides(*strides_rowcol))
-    is_shp = np.zeros(
-        (out_rows, out_cols, 2 * half_row + 1, 2 * half_col + 1), dtype=np.bool_
-    )
-    return _loop_over_pixels(
-        mean,
-        var,
-        halfwin_rowcol,
-        strides_rowcol,
-        threshold,
-        prune_disconnected,
-        is_shp,
-    )
-
-
-def get_cutoff(N: int, alpha: float) -> float:
-    """Get the pre-computed test statistic cutoff.
-
-    Parameters
-    ----------
-    N : int
-        number of SLCs used for mean/variance.
-    alpha : float
-        Significance level
-
-    Returns
-    -------
-    float
-        Threshold, above which to reject the null hypothesis.
-
-    Raises
-    ------
-    ValueError
-        If a (N, alpha) combination is passed which hasn't been precomputed.
-
-    """
-    n_alpha_to_cutoff = _read_cutoff_csv("kld")
-    try:
-        return n_alpha_to_cutoff[(N, alpha)]
-    except KeyError as e:
-        msg = f"Not implemented for {N = }, {alpha = }"
-        raise NotImplementedError(msg) from e
diff --git a/src/dolphin/shp/_ks.py b/src/dolphin/shp/_ks.py
index df8316f7..dbab3cf0 100644
--- a/src/dolphin/shp/_ks.py
+++ b/src/dolphin/shp/_ks.py
@@ -156,7 +156,7 @@ def _set_neighbors(amp_block, halfwin_rowcol, ecdf_dist_cutoff, neighbors):
     for i in range(rows):
         for j in range(cols):
             if i == r_c and j == c_c:
-                neighbors[i, j] = True
+                neighbors[i, j] = False
                 continue
             x2 = amp_block[:, i, j]
 
diff --git a/src/dolphin/shp/kld_cutoffs.csv b/src/dolphin/shp/kld_cutoffs.csv
deleted file mode 100644
index 26b54a64..00000000
--- a/src/dolphin/shp/kld_cutoffs.csv
+++ /dev/null
@@ -1,1201 +0,0 @@
-N,alpha,cutoff
-1,0.001,1102.4824
-1,0.005,191.8067
-1,0.01,97.8216
-1,0.05,18.5767
-2,0.001,34.8548
-2,0.005,15.2274
-2,0.01,10.3598
-2,0.05,3.8784
-3,0.001,12.2842
-3,0.005,6.3788
-3,0.01,4.7782
-3,0.05,2.022
-4,0.001,6.1861
-4,0.005,3.6642
-4,0.01,2.8543
-4,0.05,1.3493
-5,0.001,4.249
-5,0.005,2.6062
-5,0.01,2.017
-5,0.05,1.0
-6,0.001,3.1531
-6,0.005,1.9924
-6,0.01,1.5858
-6,0.05,0.7924
-7,0.001,2.5549
-7,0.005,1.6359
-7,0.01,1.284
-7,0.05,0.6596
-8,0.001,2.0569
-8,0.005,1.349
-8,0.01,1.0874
-8,0.05,0.5619
-9,0.001,1.6515
-9,0.005,1.151
-9,0.01,0.9333
-9,0.05,0.4857
-10,0.001,1.4608
-10,0.005,0.9865
-10,0.01,0.818
-10,0.05,0.4392
-11,0.001,1.2871
-11,0.005,0.8768
-11,0.01,0.7103
-11,0.05,0.3933
-12,0.001,1.1548
-12,0.005,0.7929
-12,0.01,0.6482
-12,0.05,0.3582
-13,0.001,1.0266
-13,0.005,0.7242
-13,0.01,0.5978
-13,0.05,0.3277
-14,0.001,0.9862
-14,0.005,0.6706
-14,0.01,0.5462
-14,0.05,0.3007
-15,0.001,0.8918
-15,0.005,0.6114
-15,0.01,0.4991
-15,0.05,0.2792
-16,0.001,0.8367
-16,0.005,0.5767
-16,0.01,0.4692
-16,0.05,0.2608
-17,0.001,0.7586
-17,0.005,0.5339
-17,0.01,0.438
-17,0.05,0.2423
-18,0.001,0.7061
-18,0.005,0.5032
-18,0.01,0.4135
-18,0.05,0.229
-19,0.001,0.6757
-19,0.005,0.4664
-19,0.01,0.3859
-19,0.05,0.2176
-20,0.001,0.609
-20,0.005,0.4485
-20,0.01,0.3676
-20,0.05,0.2083
-21,0.001,0.5973
-21,0.005,0.4173
-21,0.01,0.3456
-21,0.05,0.1926
-22,0.001,0.5707
-22,0.005,0.4028
-22,0.01,0.333
-22,0.05,0.1868
-23,0.001,0.5216
-23,0.005,0.3819
-23,0.01,0.3154
-23,0.05,0.1769
-24,0.001,0.512
-24,0.005,0.3686
-24,0.01,0.3027
-24,0.05,0.168
-25,0.001,0.5072
-25,0.005,0.3463
-25,0.01,0.2893
-25,0.05,0.1625
-26,0.001,0.4699
-26,0.005,0.3315
-26,0.01,0.2731
-26,0.05,0.1553
-27,0.001,0.4493
-27,0.005,0.3198
-27,0.01,0.2652
-27,0.05,0.1494
-28,0.001,0.4402
-28,0.005,0.3074
-28,0.01,0.255
-28,0.05,0.1434
-29,0.001,0.42
-29,0.005,0.2983
-29,0.01,0.248
-29,0.05,0.1382
-30,0.001,0.397
-30,0.005,0.2843
-30,0.01,0.2346
-30,0.05,0.1346
-31,0.001,0.3801
-31,0.005,0.2764
-31,0.01,0.2289
-31,0.05,0.1279
-32,0.001,0.3793
-32,0.005,0.2637
-32,0.01,0.2184
-32,0.05,0.1244
-33,0.001,0.3626
-33,0.005,0.2582
-33,0.01,0.2137
-33,0.05,0.1198
-34,0.001,0.3484
-34,0.005,0.2482
-34,0.01,0.2044
-34,0.05,0.1158
-35,0.001,0.3462
-35,0.005,0.2417
-35,0.01,0.2004
-35,0.05,0.1146
-36,0.001,0.3259
-36,0.005,0.2326
-36,0.01,0.1934
-36,0.05,0.1099
-37,0.001,0.3096
-37,0.005,0.2252
-37,0.01,0.1866
-37,0.05,0.1061
-38,0.001,0.3052
-38,0.005,0.2204
-38,0.01,0.1839
-38,0.05,0.1037
-39,0.001,0.3117
-39,0.005,0.2131
-39,0.01,0.176
-39,0.05,0.1001
-40,0.001,0.296
-40,0.005,0.2064
-40,0.01,0.173
-40,0.05,0.0992
-41,0.001,0.2775
-41,0.005,0.204
-41,0.01,0.1698
-41,0.05,0.0965
-42,0.001,0.2785
-42,0.005,0.1996
-42,0.01,0.1655
-42,0.05,0.0937
-43,0.001,0.2729
-43,0.005,0.1915
-43,0.01,0.1598
-43,0.05,0.0912
-44,0.001,0.2647
-44,0.005,0.1893
-44,0.01,0.1572
-44,0.05,0.0898
-45,0.001,0.257
-45,0.005,0.1793
-45,0.01,0.1525
-45,0.05,0.0874
-46,0.001,0.2512
-46,0.005,0.1771
-46,0.01,0.1492
-46,0.05,0.0846
-47,0.001,0.2451
-47,0.005,0.174
-47,0.01,0.1477
-47,0.05,0.0839
-48,0.001,0.243
-48,0.005,0.1718
-48,0.01,0.1441
-48,0.05,0.0819
-49,0.001,0.2344
-49,0.005,0.168
-49,0.01,0.1401
-49,0.05,0.0806
-50,0.001,0.2284
-50,0.005,0.1613
-50,0.01,0.1371
-50,0.05,0.0789
-51,0.001,0.2167
-51,0.005,0.1588
-51,0.01,0.1336
-51,0.05,0.0764
-52,0.001,0.2185
-52,0.005,0.1582
-52,0.01,0.1332
-52,0.05,0.0757
-53,0.001,0.2122
-53,0.005,0.1521
-53,0.01,0.1293
-53,0.05,0.0739
-54,0.001,0.2057
-54,0.005,0.1493
-54,0.01,0.1262
-54,0.05,0.072
-55,0.001,0.2043
-55,0.005,0.1456
-55,0.01,0.1224
-55,0.05,0.0716
-56,0.001,0.2008
-56,0.005,0.1446
-56,0.01,0.1234
-56,0.05,0.0709
-57,0.001,0.198
-57,0.005,0.1453
-57,0.01,0.1226
-57,0.05,0.0688
-58,0.001,0.1976
-58,0.005,0.1433
-58,0.01,0.1173
-58,0.05,0.0673
-59,0.001,0.1877
-59,0.005,0.1382
-59,0.01,0.1169
-59,0.05,0.0662
-60,0.001,0.1836
-60,0.005,0.1334
-60,0.01,0.1124
-60,0.05,0.0649
-61,0.001,0.1887
-61,0.005,0.1325
-61,0.01,0.1083
-61,0.05,0.0629
-62,0.001,0.1803
-62,0.005,0.1264
-62,0.01,0.1084
-62,0.05,0.0624
-63,0.001,0.1798
-63,0.005,0.1303
-63,0.01,0.1073
-63,0.05,0.0611
-64,0.001,0.1744
-64,0.005,0.1286
-64,0.01,0.107
-64,0.05,0.0604
-65,0.001,0.1646
-65,0.005,0.1241
-65,0.01,0.1038
-65,0.05,0.0596
-66,0.001,0.1741
-66,0.005,0.1273
-66,0.01,0.1045
-66,0.05,0.0602
-67,0.001,0.166
-67,0.005,0.1189
-67,0.01,0.1002
-67,0.05,0.0582
-68,0.001,0.169
-68,0.005,0.1204
-68,0.01,0.1004
-68,0.05,0.0583
-69,0.001,0.1703
-69,0.005,0.1153
-69,0.01,0.0967
-69,0.05,0.0562
-70,0.001,0.1609
-70,0.005,0.1141
-70,0.01,0.0966
-70,0.05,0.0551
-71,0.001,0.1548
-71,0.005,0.113
-71,0.01,0.094
-71,0.05,0.0551
-72,0.001,0.1536
-72,0.005,0.1116
-72,0.01,0.0928
-72,0.05,0.0533
-73,0.001,0.1521
-73,0.005,0.1108
-73,0.01,0.0933
-73,0.05,0.0534
-74,0.001,0.1495
-74,0.005,0.1102
-74,0.01,0.0929
-74,0.05,0.0527
-75,0.001,0.1475
-75,0.005,0.1094
-75,0.01,0.0913
-75,0.05,0.0522
-76,0.001,0.1475
-76,0.005,0.1054
-76,0.01,0.0882
-76,0.05,0.0504
-77,0.001,0.1512
-77,0.005,0.1044
-77,0.01,0.0879
-77,0.05,0.0506
-78,0.001,0.1368
-78,0.005,0.102
-78,0.01,0.087
-78,0.05,0.0501
-79,0.001,0.1394
-79,0.005,0.1004
-79,0.01,0.0839
-79,0.05,0.0492
-80,0.001,0.1321
-80,0.005,0.101
-80,0.01,0.0843
-80,0.05,0.0479
-81,0.001,0.1337
-81,0.005,0.1003
-81,0.01,0.0851
-81,0.05,0.048
-82,0.001,0.1349
-82,0.005,0.0982
-82,0.01,0.0818
-82,0.05,0.0475
-83,0.001,0.1345
-83,0.005,0.0985
-83,0.01,0.0812
-83,0.05,0.0466
-84,0.001,0.1239
-84,0.005,0.0953
-84,0.01,0.0801
-84,0.05,0.0459
-85,0.001,0.1299
-85,0.005,0.0949
-85,0.01,0.0795
-85,0.05,0.046
-86,0.001,0.1282
-86,0.005,0.0925
-86,0.01,0.0772
-86,0.05,0.0449
-87,0.001,0.1246
-87,0.005,0.0908
-87,0.01,0.0757
-87,0.05,0.0448
-88,0.001,0.1216
-88,0.005,0.0912
-88,0.01,0.077
-88,0.05,0.0446
-89,0.001,0.121
-89,0.005,0.0904
-89,0.01,0.0764
-89,0.05,0.0437
-90,0.001,0.1255
-90,0.005,0.0879
-90,0.01,0.0744
-90,0.05,0.0431
-91,0.001,0.1282
-91,0.005,0.0895
-91,0.01,0.0741
-91,0.05,0.0426
-92,0.001,0.1154
-92,0.005,0.0867
-92,0.01,0.0728
-92,0.05,0.0422
-93,0.001,0.1208
-93,0.005,0.0864
-93,0.01,0.0725
-93,0.05,0.0425
-94,0.001,0.1194
-94,0.005,0.0868
-94,0.01,0.073
-94,0.05,0.0415
-95,0.001,0.1169
-95,0.005,0.0853
-95,0.01,0.0705
-95,0.05,0.0408
-96,0.001,0.1121
-96,0.005,0.0837
-96,0.01,0.0692
-96,0.05,0.04
-97,0.001,0.1134
-97,0.005,0.0832
-97,0.01,0.0687
-97,0.05,0.0397
-98,0.001,0.1129
-98,0.005,0.0827
-98,0.01,0.0692
-98,0.05,0.0401
-99,0.001,0.1132
-99,0.005,0.0788
-99,0.01,0.0674
-99,0.05,0.0387
-100,0.001,0.1083
-100,0.005,0.0787
-100,0.01,0.0658
-100,0.05,0.0383
-101,0.001,0.1077
-101,0.005,0.0797
-101,0.01,0.0658
-101,0.05,0.0381
-102,0.001,0.1083
-102,0.005,0.0771
-102,0.01,0.0653
-102,0.05,0.038
-103,0.001,0.1024
-103,0.005,0.0773
-103,0.01,0.0648
-103,0.05,0.0378
-104,0.001,0.1053
-104,0.005,0.0758
-104,0.01,0.0644
-104,0.05,0.0371
-105,0.001,0.1022
-105,0.005,0.0758
-105,0.01,0.064
-105,0.05,0.0367
-106,0.001,0.1045
-106,0.005,0.0759
-106,0.01,0.064
-106,0.05,0.0371
-107,0.001,0.1015
-107,0.005,0.0748
-107,0.01,0.0626
-107,0.05,0.0363
-108,0.001,0.1057
-108,0.005,0.0748
-108,0.01,0.0632
-108,0.05,0.0359
-109,0.001,0.1053
-109,0.005,0.0759
-109,0.01,0.0629
-109,0.05,0.0357
-110,0.001,0.1012
-110,0.005,0.0727
-110,0.01,0.061
-110,0.05,0.0346
-111,0.001,0.0968
-111,0.005,0.0722
-111,0.01,0.0608
-111,0.05,0.0348
-112,0.001,0.1
-112,0.005,0.0712
-112,0.01,0.0603
-112,0.05,0.0348
-113,0.001,0.0961
-113,0.005,0.0701
-113,0.01,0.059
-113,0.05,0.0339
-114,0.001,0.0944
-114,0.005,0.07
-114,0.01,0.0592
-114,0.05,0.034
-115,0.001,0.0953
-115,0.005,0.0696
-115,0.01,0.0582
-115,0.05,0.0336
-116,0.001,0.0937
-116,0.005,0.0684
-116,0.01,0.0581
-116,0.05,0.0336
-117,0.001,0.0943
-117,0.005,0.0676
-117,0.01,0.0577
-117,0.05,0.0332
-118,0.001,0.0913
-118,0.005,0.0689
-118,0.01,0.0576
-118,0.05,0.033
-119,0.001,0.0895
-119,0.005,0.0654
-119,0.01,0.0562
-119,0.05,0.0323
-120,0.001,0.0917
-120,0.005,0.0658
-120,0.01,0.0556
-120,0.05,0.0323
-121,0.001,0.0933
-121,0.005,0.0663
-121,0.01,0.0562
-121,0.05,0.0321
-122,0.001,0.0882
-122,0.005,0.0653
-122,0.01,0.0544
-122,0.05,0.0316
-123,0.001,0.086
-123,0.005,0.0651
-123,0.01,0.0544
-123,0.05,0.0319
-124,0.001,0.0886
-124,0.005,0.0662
-124,0.01,0.0543
-124,0.05,0.0308
-125,0.001,0.0891
-125,0.005,0.064
-125,0.01,0.0539
-125,0.05,0.0309
-126,0.001,0.0852
-126,0.005,0.0618
-126,0.01,0.0528
-126,0.05,0.0307
-127,0.001,0.0819
-127,0.005,0.0621
-127,0.01,0.0522
-127,0.05,0.0306
-128,0.001,0.0828
-128,0.005,0.062
-128,0.01,0.0525
-128,0.05,0.0301
-129,0.001,0.0872
-129,0.005,0.0619
-129,0.01,0.0515
-129,0.05,0.03
-130,0.001,0.0837
-130,0.005,0.0612
-130,0.01,0.0505
-130,0.05,0.0293
-131,0.001,0.0816
-131,0.005,0.0607
-131,0.01,0.0502
-131,0.05,0.0294
-132,0.001,0.0869
-132,0.005,0.0606
-132,0.01,0.0504
-132,0.05,0.0291
-133,0.001,0.0816
-133,0.005,0.0605
-133,0.01,0.0509
-133,0.05,0.0291
-134,0.001,0.0816
-134,0.005,0.0593
-134,0.01,0.0499
-134,0.05,0.0291
-135,0.001,0.0823
-135,0.005,0.0588
-135,0.01,0.0493
-135,0.05,0.0288
-136,0.001,0.0795
-136,0.005,0.0571
-136,0.01,0.048
-136,0.05,0.0284
-137,0.001,0.0772
-137,0.005,0.0575
-137,0.01,0.0486
-137,0.05,0.0282
-138,0.001,0.0805
-138,0.005,0.0574
-138,0.01,0.0487
-138,0.05,0.028
-139,0.001,0.0763
-139,0.005,0.0563
-139,0.01,0.0481
-139,0.05,0.028
-140,0.001,0.0804
-140,0.005,0.0578
-140,0.01,0.0483
-140,0.05,0.0276
-141,0.001,0.0767
-141,0.005,0.0554
-141,0.01,0.0472
-141,0.05,0.0269
-142,0.001,0.0808
-142,0.005,0.0563
-142,0.01,0.0476
-142,0.05,0.027
-143,0.001,0.076
-143,0.005,0.0557
-143,0.01,0.0465
-143,0.05,0.0272
-144,0.001,0.0749
-144,0.005,0.0552
-144,0.01,0.0467
-144,0.05,0.0267
-145,0.001,0.074
-145,0.005,0.0531
-145,0.01,0.0454
-145,0.05,0.0264
-146,0.001,0.0746
-146,0.005,0.054
-146,0.01,0.0451
-146,0.05,0.0263
-147,0.001,0.0753
-147,0.005,0.053
-147,0.01,0.0455
-147,0.05,0.0261
-148,0.001,0.0757
-148,0.005,0.0546
-148,0.01,0.0455
-148,0.05,0.0259
-149,0.001,0.0729
-149,0.005,0.0534
-149,0.01,0.0441
-149,0.05,0.0259
-150,0.001,0.0706
-150,0.005,0.0527
-150,0.01,0.0447
-150,0.05,0.0258
-151,0.001,0.0752
-151,0.005,0.053
-151,0.01,0.0449
-151,0.05,0.0257
-152,0.001,0.0725
-152,0.005,0.0528
-152,0.01,0.0442
-152,0.05,0.0255
-153,0.001,0.0716
-153,0.005,0.0524
-153,0.01,0.0437
-153,0.05,0.0253
-154,0.001,0.0721
-154,0.005,0.0519
-154,0.01,0.0434
-154,0.05,0.0251
-155,0.001,0.0723
-155,0.005,0.0519
-155,0.01,0.0433
-155,0.05,0.025
-156,0.001,0.0705
-156,0.005,0.0515
-156,0.01,0.043
-156,0.05,0.0249
-157,0.001,0.068
-157,0.005,0.0507
-157,0.01,0.0416
-157,0.05,0.0244
-158,0.001,0.0687
-158,0.005,0.0508
-158,0.01,0.0423
-158,0.05,0.0242
-159,0.001,0.0699
-159,0.005,0.0494
-159,0.01,0.0417
-159,0.05,0.0239
-160,0.001,0.0702
-160,0.005,0.0494
-160,0.01,0.041
-160,0.05,0.024
-161,0.001,0.0668
-161,0.005,0.0492
-161,0.01,0.0412
-161,0.05,0.0238
-162,0.001,0.0693
-162,0.005,0.0492
-162,0.01,0.041
-162,0.05,0.0236
-163,0.001,0.0646
-163,0.005,0.0475
-163,0.01,0.041
-163,0.05,0.0236
-164,0.001,0.0672
-164,0.005,0.0499
-164,0.01,0.0415
-164,0.05,0.0235
-165,0.001,0.0673
-165,0.005,0.0481
-165,0.01,0.0404
-165,0.05,0.0232
-166,0.001,0.0667
-166,0.005,0.0487
-166,0.01,0.0404
-166,0.05,0.023
-167,0.001,0.0648
-167,0.005,0.0476
-167,0.01,0.0406
-167,0.05,0.023
-168,0.001,0.0645
-168,0.005,0.0481
-168,0.01,0.0395
-168,0.05,0.0232
-169,0.001,0.0637
-169,0.005,0.0468
-169,0.01,0.0389
-169,0.05,0.0229
-170,0.001,0.066
-170,0.005,0.0473
-170,0.01,0.0401
-170,0.05,0.0228
-171,0.001,0.0635
-171,0.005,0.0464
-171,0.01,0.039
-171,0.05,0.0227
-172,0.001,0.063
-172,0.005,0.0466
-172,0.01,0.0385
-172,0.05,0.0223
-173,0.001,0.0665
-173,0.005,0.0471
-173,0.01,0.0383
-173,0.05,0.0222
-174,0.001,0.0627
-174,0.005,0.0454
-174,0.01,0.0384
-174,0.05,0.0222
-175,0.001,0.063
-175,0.005,0.0449
-175,0.01,0.0375
-175,0.05,0.022
-176,0.001,0.0611
-176,0.005,0.0456
-176,0.01,0.0381
-176,0.05,0.0218
-177,0.001,0.0616
-177,0.005,0.045
-177,0.01,0.0378
-177,0.05,0.022
-178,0.001,0.0656
-178,0.005,0.0455
-178,0.01,0.0381
-178,0.05,0.0215
-179,0.001,0.0602
-179,0.005,0.0431
-179,0.01,0.0364
-179,0.05,0.0215
-180,0.001,0.0616
-180,0.005,0.0441
-180,0.01,0.0371
-180,0.05,0.0214
-181,0.001,0.06
-181,0.005,0.0439
-181,0.01,0.0369
-181,0.05,0.0212
-182,0.001,0.0615
-182,0.005,0.0443
-182,0.01,0.0373
-182,0.05,0.0215
-183,0.001,0.0614
-183,0.005,0.0431
-183,0.01,0.0368
-183,0.05,0.0209
-184,0.001,0.0582
-184,0.005,0.0438
-184,0.01,0.037
-184,0.05,0.021
-185,0.001,0.0581
-185,0.005,0.0433
-185,0.01,0.0366
-185,0.05,0.0206
-186,0.001,0.0581
-186,0.005,0.0425
-186,0.01,0.0356
-186,0.05,0.0205
-187,0.001,0.0587
-187,0.005,0.0422
-187,0.01,0.0356
-187,0.05,0.0205
-188,0.001,0.0585
-188,0.005,0.0421
-188,0.01,0.0357
-188,0.05,0.0205
-189,0.001,0.0577
-189,0.005,0.0425
-189,0.01,0.0355
-189,0.05,0.0202
-190,0.001,0.0575
-190,0.005,0.0416
-190,0.01,0.0349
-190,0.05,0.0198
-191,0.001,0.0573
-191,0.005,0.042
-191,0.01,0.0347
-191,0.05,0.0201
-192,0.001,0.0563
-192,0.005,0.0416
-192,0.01,0.035
-192,0.05,0.0198
-193,0.001,0.057
-193,0.005,0.041
-193,0.01,0.0346
-193,0.05,0.0202
-194,0.001,0.0578
-194,0.005,0.042
-194,0.01,0.0349
-194,0.05,0.0199
-195,0.001,0.0582
-195,0.005,0.0411
-195,0.01,0.0341
-195,0.05,0.0198
-196,0.001,0.0567
-196,0.005,0.0399
-196,0.01,0.0338
-196,0.05,0.0197
-197,0.001,0.0547
-197,0.005,0.0399
-197,0.01,0.0334
-197,0.05,0.0195
-198,0.001,0.0555
-198,0.005,0.0398
-198,0.01,0.0338
-198,0.05,0.0195
-199,0.001,0.055
-199,0.005,0.0391
-199,0.01,0.0332
-199,0.05,0.0192
-200,0.001,0.0555
-200,0.005,0.0398
-200,0.01,0.0332
-200,0.05,0.0193
-201,0.001,0.0554
-201,0.005,0.0404
-201,0.01,0.0335
-201,0.05,0.0193
-202,0.001,0.0553
-202,0.005,0.0393
-202,0.01,0.0331
-202,0.05,0.019
-203,0.001,0.0548
-203,0.005,0.0405
-203,0.01,0.0335
-203,0.05,0.019
-204,0.001,0.0541
-204,0.005,0.0387
-204,0.01,0.033
-204,0.05,0.0189
-205,0.001,0.053
-205,0.005,0.0386
-205,0.01,0.0326
-205,0.05,0.0187
-206,0.001,0.0542
-206,0.005,0.0394
-206,0.01,0.0328
-206,0.05,0.0186
-207,0.001,0.0522
-207,0.005,0.0386
-207,0.01,0.032
-207,0.05,0.0185
-208,0.001,0.0528
-208,0.005,0.0382
-208,0.01,0.0322
-208,0.05,0.0185
-209,0.001,0.0521
-209,0.005,0.0384
-209,0.01,0.0321
-209,0.05,0.0184
-210,0.001,0.0533
-210,0.005,0.038
-210,0.01,0.0318
-210,0.05,0.018
-211,0.001,0.0529
-211,0.005,0.0381
-211,0.01,0.032
-211,0.05,0.0182
-212,0.001,0.0525
-212,0.005,0.0376
-212,0.01,0.0317
-212,0.05,0.0179
-213,0.001,0.051
-213,0.005,0.0368
-213,0.01,0.0313
-213,0.05,0.018
-214,0.001,0.0521
-214,0.005,0.0366
-214,0.01,0.031
-214,0.05,0.0179
-215,0.001,0.0519
-215,0.005,0.0367
-215,0.01,0.0308
-215,0.05,0.0177
-216,0.001,0.0509
-216,0.005,0.037
-216,0.01,0.031
-216,0.05,0.0177
-217,0.001,0.0505
-217,0.005,0.0355
-217,0.01,0.0302
-217,0.05,0.0175
-218,0.001,0.0521
-218,0.005,0.0372
-218,0.01,0.0309
-218,0.05,0.0176
-219,0.001,0.0492
-219,0.005,0.0361
-219,0.01,0.0307
-219,0.05,0.0174
-220,0.001,0.0518
-220,0.005,0.036
-220,0.01,0.0301
-220,0.05,0.0173
-221,0.001,0.0492
-221,0.005,0.0361
-221,0.01,0.0304
-221,0.05,0.0172
-222,0.001,0.0503
-222,0.005,0.0366
-222,0.01,0.0304
-222,0.05,0.0172
-223,0.001,0.0497
-223,0.005,0.0356
-223,0.01,0.0299
-223,0.05,0.0173
-224,0.001,0.0487
-224,0.005,0.0349
-224,0.01,0.0296
-224,0.05,0.0174
-225,0.001,0.0481
-225,0.005,0.0345
-225,0.01,0.0293
-225,0.05,0.0169
-226,0.001,0.0493
-226,0.005,0.0348
-226,0.01,0.0291
-226,0.05,0.0169
-227,0.001,0.0454
-227,0.005,0.0333
-227,0.01,0.0285
-227,0.05,0.0166
-228,0.001,0.0479
-228,0.005,0.0349
-228,0.01,0.0296
-228,0.05,0.0168
-229,0.001,0.0474
-229,0.005,0.0335
-229,0.01,0.0288
-229,0.05,0.0165
-230,0.001,0.0471
-230,0.005,0.0334
-230,0.01,0.0286
-230,0.05,0.0167
-231,0.001,0.0491
-231,0.005,0.034
-231,0.01,0.0288
-231,0.05,0.0166
-232,0.001,0.0473
-232,0.005,0.0349
-232,0.01,0.0287
-232,0.05,0.0166
-233,0.001,0.0465
-233,0.005,0.0337
-233,0.01,0.0287
-233,0.05,0.0165
-234,0.001,0.0469
-234,0.005,0.0335
-234,0.01,0.0281
-234,0.05,0.0163
-235,0.001,0.0454
-235,0.005,0.0332
-235,0.01,0.0283
-235,0.05,0.0164
-236,0.001,0.046
-236,0.005,0.0332
-236,0.01,0.028
-236,0.05,0.016
-237,0.001,0.0455
-237,0.005,0.033
-237,0.01,0.0278
-237,0.05,0.0164
-238,0.001,0.0468
-238,0.005,0.0331
-238,0.01,0.028
-238,0.05,0.016
-239,0.001,0.0462
-239,0.005,0.033
-239,0.01,0.0277
-239,0.05,0.0161
-240,0.001,0.0461
-240,0.005,0.033
-240,0.01,0.0278
-240,0.05,0.016
-241,0.001,0.0456
-241,0.005,0.0325
-241,0.01,0.0274
-241,0.05,0.0159
-242,0.001,0.0468
-242,0.005,0.0334
-242,0.01,0.0273
-242,0.05,0.0158
-243,0.001,0.0469
-243,0.005,0.0331
-243,0.01,0.0276
-243,0.05,0.016
-244,0.001,0.046
-244,0.005,0.0328
-244,0.01,0.0273
-244,0.05,0.0156
-245,0.001,0.0454
-245,0.005,0.0322
-245,0.01,0.027
-245,0.05,0.0157
-246,0.001,0.0449
-246,0.005,0.0324
-246,0.01,0.0272
-246,0.05,0.0156
-247,0.001,0.0448
-247,0.005,0.0325
-247,0.01,0.027
-247,0.05,0.0157
-248,0.001,0.0438
-248,0.005,0.0313
-248,0.01,0.0265
-248,0.05,0.0154
-249,0.001,0.0445
-249,0.005,0.031
-249,0.01,0.0263
-249,0.05,0.0156
-250,0.001,0.0444
-250,0.005,0.0318
-250,0.01,0.0268
-250,0.05,0.0156
-251,0.001,0.0426
-251,0.005,0.0319
-251,0.01,0.0262
-251,0.05,0.0151
-252,0.001,0.0434
-252,0.005,0.0315
-252,0.01,0.0262
-252,0.05,0.0152
-253,0.001,0.0429
-253,0.005,0.0309
-253,0.01,0.0262
-253,0.05,0.0151
-254,0.001,0.0429
-254,0.005,0.0306
-254,0.01,0.0256
-254,0.05,0.015
-255,0.001,0.0404
-255,0.005,0.0305
-255,0.01,0.0258
-255,0.05,0.0152
-256,0.001,0.0438
-256,0.005,0.031
-256,0.01,0.0257
-256,0.05,0.0151
-257,0.001,0.0389
-257,0.005,0.0304
-257,0.01,0.0257
-257,0.05,0.0149
-258,0.001,0.0431
-258,0.005,0.031
-258,0.01,0.0256
-258,0.05,0.0148
-259,0.001,0.0413
-259,0.005,0.0297
-259,0.01,0.0253
-259,0.05,0.0147
-260,0.001,0.0427
-260,0.005,0.0308
-260,0.01,0.0258
-260,0.05,0.0148
-261,0.001,0.0429
-261,0.005,0.0303
-261,0.01,0.0255
-261,0.05,0.0147
-262,0.001,0.0408
-262,0.005,0.0299
-262,0.01,0.0252
-262,0.05,0.0145
-263,0.001,0.0414
-263,0.005,0.03
-263,0.01,0.0257
-263,0.05,0.0146
-264,0.001,0.0405
-264,0.005,0.0297
-264,0.01,0.0251
-264,0.05,0.0145
-265,0.001,0.0401
-265,0.005,0.0296
-265,0.01,0.0248
-265,0.05,0.0145
-266,0.001,0.0409
-266,0.005,0.0295
-266,0.01,0.0249
-266,0.05,0.0143
-267,0.001,0.0408
-267,0.005,0.0301
-267,0.01,0.0254
-267,0.05,0.0142
-268,0.001,0.0416
-268,0.005,0.0294
-268,0.01,0.0247
-268,0.05,0.0143
-269,0.001,0.0402
-269,0.005,0.0291
-269,0.01,0.0248
-269,0.05,0.0144
-270,0.001,0.0384
-270,0.005,0.0289
-270,0.01,0.0242
-270,0.05,0.0141
-271,0.001,0.0406
-271,0.005,0.029
-271,0.01,0.0245
-271,0.05,0.0142
-272,0.001,0.0404
-272,0.005,0.0293
-272,0.01,0.0248
-272,0.05,0.0142
-273,0.001,0.0396
-273,0.005,0.0288
-273,0.01,0.0243
-273,0.05,0.014
-274,0.001,0.0416
-274,0.005,0.0293
-274,0.01,0.0248
-274,0.05,0.0141
-275,0.001,0.04
-275,0.005,0.0289
-275,0.01,0.0243
-275,0.05,0.014
-276,0.001,0.0404
-276,0.005,0.029
-276,0.01,0.0242
-276,0.05,0.0139
-277,0.001,0.0387
-277,0.005,0.0285
-277,0.01,0.024
-277,0.05,0.0139
-278,0.001,0.0398
-278,0.005,0.0283
-278,0.01,0.0236
-278,0.05,0.0137
-279,0.001,0.0397
-279,0.005,0.0276
-279,0.01,0.0234
-279,0.05,0.0137
-280,0.001,0.0393
-280,0.005,0.0289
-280,0.01,0.0238
-280,0.05,0.0137
-281,0.001,0.0385
-281,0.005,0.0283
-281,0.01,0.0237
-281,0.05,0.0136
-282,0.001,0.0382
-282,0.005,0.0288
-282,0.01,0.0238
-282,0.05,0.0135
-283,0.001,0.0402
-283,0.005,0.0281
-283,0.01,0.0233
-283,0.05,0.0136
-284,0.001,0.038
-284,0.005,0.0272
-284,0.01,0.0228
-284,0.05,0.0134
-285,0.001,0.0391
-285,0.005,0.0278
-285,0.01,0.0229
-285,0.05,0.0134
-286,0.001,0.0384
-286,0.005,0.0278
-286,0.01,0.0234
-286,0.05,0.0136
-287,0.001,0.0377
-287,0.005,0.0279
-287,0.01,0.0232
-287,0.05,0.0136
-288,0.001,0.0387
-288,0.005,0.0277
-288,0.01,0.0229
-288,0.05,0.0131
-289,0.001,0.0373
-289,0.005,0.0268
-289,0.01,0.0228
-289,0.05,0.0131
-290,0.001,0.0384
-290,0.005,0.0275
-290,0.01,0.0227
-290,0.05,0.0134
-291,0.001,0.0371
-291,0.005,0.0273
-291,0.01,0.0228
-291,0.05,0.0131
-292,0.001,0.0389
-292,0.005,0.0274
-292,0.01,0.0226
-292,0.05,0.013
-293,0.001,0.0378
-293,0.005,0.027
-293,0.01,0.0226
-293,0.05,0.013
-294,0.001,0.0383
-294,0.005,0.0275
-294,0.01,0.0233
-294,0.05,0.0131
-295,0.001,0.0367
-295,0.005,0.0265
-295,0.01,0.0225
-295,0.05,0.0129
-296,0.001,0.0363
-296,0.005,0.0264
-296,0.01,0.0222
-296,0.05,0.0129
-297,0.001,0.0374
-297,0.005,0.0267
-297,0.01,0.0224
-297,0.05,0.0131
-298,0.001,0.0363
-298,0.005,0.026
-298,0.01,0.022
-298,0.05,0.013
-299,0.001,0.0357
-299,0.005,0.0267
-299,0.01,0.0222
-299,0.05,0.0128
-300,0.001,0.0359
-300,0.005,0.0261
-300,0.01,0.0223
-300,0.05,0.0127
diff --git a/src/dolphin/workflows/config/_enums.py b/src/dolphin/workflows/config/_enums.py
index 178c946c..f9d857c6 100644
--- a/src/dolphin/workflows/config/_enums.py
+++ b/src/dolphin/workflows/config/_enums.py
@@ -12,7 +12,6 @@ class ShpMethod(str, Enum):
 
     GLRT = "glrt"
     KS = "ks"
-    KLD = "kld"
     RECT = "rect"
     # Alias for no SHP search
     NONE = "rect"
diff --git a/src/dolphin/workflows/displacement.py b/src/dolphin/workflows/displacement.py
index 2dcc0242..5af6a48d 100755
--- a/src/dolphin/workflows/displacement.py
+++ b/src/dolphin/workflows/displacement.py
@@ -35,6 +35,7 @@ class OutputPaths(NamedTuple):
     stitched_temp_coh_file: Path
     stitched_ps_file: Path
     stitched_amp_dispersion_file: Path
+    stitched_shp_count_file: Path
     unwrapped_paths: list[Path] | None
     conncomp_paths: list[Path] | None
     tropospheric_corrections: list[Path] | None
@@ -136,6 +137,7 @@ def run(
     temp_coh_file_list: list[Path] = []
     ps_file_list: list[Path] = []
     amp_dispersion_file_list: list[Path] = []
+    shp_count_file_list: list[Path] = []
     # The comp_slc tracking object is a dict, since we'll need to organize
     # multiple comp slcs by burst (they'll have the same filename)
     comp_slc_dict: dict[str, list[Path]] = {}
@@ -169,12 +171,15 @@ def run(
         for fut in fut_to_burst:
             burst = fut_to_burst[fut]
 
-            cur_ifg_list, comp_slcs, temp_coh, ps_file, amp_disp_file = fut.result()
+            cur_ifg_list, comp_slcs, temp_coh, ps_file, amp_disp_file, shp_count = (
+                fut.result()
+            )
             ifg_file_list.extend(cur_ifg_list)
             comp_slc_dict[burst] = comp_slcs
             temp_coh_file_list.append(temp_coh)
             ps_file_list.append(ps_file)
             amp_dispersion_file_list.append(amp_disp_file)
+            shp_count_file_list.append(shp_count)
 
     # ###################################
     # 2. Stitch burst-wise interferograms
@@ -190,11 +195,13 @@ def run(
         stitched_temp_coh_file,
         stitched_ps_file,
         stitched_amp_dispersion_file,
+        stitched_shp_count_file,
     ) = stitching_bursts.run(
         ifg_file_list=ifg_file_list,
         temp_coh_file_list=temp_coh_file_list,
         ps_file_list=ps_file_list,
         amp_dispersion_list=amp_dispersion_file_list,
+        shp_count_file_list=shp_count_file_list,
         stitched_ifg_dir=cfg.interferogram_network._directory,
         output_options=cfg.output_options,
         file_date_fmt=cfg.input_options.cslc_date_fmt,
@@ -214,6 +221,7 @@ def run(
             stitched_temp_coh_file=stitched_temp_coh_file,
             stitched_ps_file=stitched_ps_file,
             stitched_amp_dispersion_file=stitched_amp_dispersion_file,
+            stitched_shp_count_file=stitched_shp_count_file,
             unwrapped_paths=None,
             conncomp_paths=None,
             tropospheric_corrections=None,
@@ -328,6 +336,7 @@ def run(
         stitched_temp_coh_file=stitched_temp_coh_file,
         stitched_ps_file=stitched_ps_file,
         stitched_amp_dispersion_file=stitched_amp_dispersion_file,
+        stitched_shp_count_file=stitched_shp_count_file,
         unwrapped_paths=unwrapped_paths,
         # TODO: Let's keep the uwrapped_paths since all the outputs are
         # corresponding to those and if we have a network unwrapping, the
diff --git a/src/dolphin/workflows/sequential.py b/src/dolphin/workflows/sequential.py
index 92c24e6d..f0820dc4 100644
--- a/src/dolphin/workflows/sequential.py
+++ b/src/dolphin/workflows/sequential.py
@@ -44,7 +44,7 @@ def run_wrapped_phase_sequential(
     beta: float = 0.00,
     block_shape: tuple[int, int] = (512, 512),
     **tqdm_kwargs,
-) -> tuple[list[Path], list[Path], Path]:
+) -> tuple[list[Path], list[Path], Path, Path]:
     """Estimate wrapped phase using batches of ministacks."""
     if strides is None:
         strides = {"x": 1, "y": 1}
@@ -62,8 +62,9 @@ def run_wrapped_phase_sequential(
 
     # list where each item is [output_slc_files] from a ministack
     output_slc_files: list[list] = []
-    # Each item is the temp_coh file from a ministack
+    # Each item is the temp_coh/shp_count file from a ministack
     temp_coh_files: list[Path] = []
+    shp_count_files: list[Path] = []
 
     # function to check if a ministack has already been processed
     def already_processed(d: Path, search_ext: str = ".tif") -> bool:
@@ -113,34 +114,24 @@ def already_processed(d: Path, search_ext: str = ".tif") -> bool:
                 **tqdm_kwargs,
             )
 
-        cur_output_files, cur_comp_slc_file, temp_coh_file = _get_outputs_from_folder(
-            cur_output_folder
+        cur_output_files, cur_comp_slc_file, temp_coh_file, shp_count_file = (
+            _get_outputs_from_folder(cur_output_folder)
         )
         output_slc_files.append(cur_output_files)
         temp_coh_files.append(temp_coh_file)
+        shp_count_files.append(shp_count_file)
 
     ##############################################
 
     # Average the temporal coherence files in each ministack
     full_span = ministack_planner.real_slc_date_range_str
     output_temp_coh_file = output_folder / f"temporal_coherence_average_{full_span}.tif"
+    output_shp_count_file = output_folder / f"shp_counts_average_{full_span}.tif"
+
     # we can pass the list of files to gdal_calc, which interprets it
     # as a multi-band file
-    if len(temp_coh_files) > 1:
-        logger.info(f"Averaging temporal coherence files into: {output_temp_coh_file}")
-        gdal_calc.Calc(
-            NoDataValue=0,
-            format="GTiff",
-            outfile=fspath(output_temp_coh_file),
-            type="Float32",
-            quiet=True,
-            overwrite=True,
-            creation_options=io.DEFAULT_TIFF_OPTIONS,
-            A=temp_coh_files,
-            calc="numpy.nanmean(A, axis=0)",
-        )
-    else:
-        temp_coh_files[0].rename(output_temp_coh_file)
+    _average_rasters(temp_coh_files, output_temp_coh_file, "Float32")
+    _average_rasters(shp_count_files, output_shp_count_file, "Int16")
 
     # Combine the separate SLC output lists into a single list
     all_slc_files = list(chain.from_iterable(output_slc_files))
@@ -151,16 +142,40 @@ def already_processed(d: Path, search_ext: str = ".tif") -> bool:
         slc_fname.rename(output_folder / slc_fname.name)
         out_pl_slcs.append(output_folder / slc_fname.name)
 
-    comp_outputs = []
+    comp_slc_outputs = []
     for p in all_comp_slc_files:
         p.rename(output_folder / p.name)
-        comp_outputs.append(output_folder / p.name)
+        comp_slc_outputs.append(output_folder / p.name)
 
-    return out_pl_slcs, comp_outputs, output_temp_coh_file
+    return out_pl_slcs, comp_slc_outputs, output_temp_coh_file, output_shp_count_file
 
 
-def _get_outputs_from_folder(output_folder: Path):
+def _get_outputs_from_folder(
+    output_folder: Path,
+) -> tuple[list[Path], Path, Path, Path]:
     cur_output_files = sorted(output_folder.glob("2*.slc.tif"))
     cur_comp_slc_file = next(output_folder.glob("compressed_*"))
     temp_coh_file = next(output_folder.glob("temporal_coherence_*"))
-    return cur_output_files, cur_comp_slc_file, temp_coh_file
+    shp_count_file = next(output_folder.glob("shp_counts_*"))
+    # Currently ignoring to not stitch:
+    # eigenvalues, estimator, avg_coh
+    return cur_output_files, cur_comp_slc_file, temp_coh_file, shp_count_file
+
+
+def _average_rasters(file_list: list[Path], outfile: Path, output_type: str):
+    if len(file_list) == 1:
+        file_list[0].rename(outfile)
+        return
+
+    logger.info(f"Averaging {len(file_list)} files into {outfile}")
+    gdal_calc.Calc(
+        NoDataValue=0,
+        format="GTiff",
+        outfile=fspath(outfile),
+        type=output_type,
+        quiet=True,
+        overwrite=True,
+        creation_options=io.DEFAULT_TIFF_OPTIONS,
+        A=file_list,
+        calc="numpy.nanmean(A, axis=0)",
+    )
diff --git a/src/dolphin/workflows/single.py b/src/dolphin/workflows/single.py
index e43a7ba1..24b93b4a 100644
--- a/src/dolphin/workflows/single.py
+++ b/src/dolphin/workflows/single.py
@@ -119,10 +119,10 @@ def run_wrapped_phase_single(
         OutputFile(
             output_folder / f"temporal_coherence_{start_end}.tif", np.float32, strides
         ),
+        OutputFile(output_folder / f"shp_counts_{start_end}.tif", np.uint16, strides),
         OutputFile(output_folder / f"eigenvalues_{start_end}.tif", np.float32, strides),
         OutputFile(output_folder / f"estimator_{start_end}.tif", np.int8, strides),
         OutputFile(output_folder / f"avg_coh_{start_end}.tif", np.uint16, strides),
-        OutputFile(output_folder / f"shp_counts_{start_end}.tif", np.uint16, strides),
     ]
     for op in output_files:
         io.write_arr(
@@ -230,12 +230,6 @@ def run_wrapped_phase_single(
         ):
             writer.queue_write(img, f, out_rows.start, out_cols.start)
 
-        # Get the SHP counts for each pixel (if not using Rect window)
-        if neighbor_arrays is None:
-            shp_counts = np.zeros(pl_output.cpx_phase.shape[-2:], dtype=np.int16)
-        else:
-            shp_counts = np.sum(neighbor_arrays, axis=(-2, -1))
-
         # Compress the ministack using only the non-compressed SLCs
         cur_comp_slc = compress(
             # Get the inner portion of the full-res SLC data
@@ -256,10 +250,10 @@ def run_wrapped_phase_single(
         # All other outputs are strided (smaller in size)
         out_datas = [
             pl_output.temp_coh,
+            pl_output.shp_counts,
             pl_output.eigenvalues,
             pl_output.estimator,
             pl_output.avg_coh,
-            shp_counts,
         ]
         for data, output_file in zip(out_datas, output_files[1:]):
             if data is None:  # May choose to skip some outputs, e.g. "avg_coh"
diff --git a/src/dolphin/workflows/stitching_bursts.py b/src/dolphin/workflows/stitching_bursts.py
index 75cb64d4..5ed6e4a3 100644
--- a/src/dolphin/workflows/stitching_bursts.py
+++ b/src/dolphin/workflows/stitching_bursts.py
@@ -22,11 +22,12 @@ def run(
     temp_coh_file_list: Sequence[Path],
     ps_file_list: Sequence[Path],
     amp_dispersion_list: Sequence[Path],
+    shp_count_file_list: Sequence[Path],
     stitched_ifg_dir: Path,
     output_options: OutputOptions,
     file_date_fmt: str = "%Y%m%d",
     corr_window_size: tuple[int, int] = (11, 11),
-) -> tuple[list[Path], list[Path], Path, Path, Path]:
+) -> tuple[list[Path], list[Path], Path, Path, Path, Path]:
     """Run the displacement workflow on a stack of SLCs.
 
     Parameters
@@ -39,6 +40,8 @@ def run(
         Sequence of paths to the (looked) burst-wise ps mask files.
     amp_dispersion_list : Sequence[Path]
         Sequence of paths to the (looked) burst-wise amplitude dispersion files.
+    shp_count_file_list : Sequence[Path]
+        Sequence of paths to the burst-wise SHP counts files.
     stitched_ifg_dir : Path
         Location to store the output stitched ifgs and correlations
     output_options : OutputOptions
@@ -63,6 +66,8 @@ def run(
         Path to ps mask file created.
     stitched_amp_disp_file : Path
         Path to amplitude dispersion file created.
+    stitched_shp_count_file : Path
+        Path to SHP count file created.
 
     """
     stitched_ifg_dir.mkdir(exist_ok=True, parents=True)
@@ -117,6 +122,15 @@ def run(
         out_bounds_epsg=output_options.bounds_epsg,
     )
 
+    stitched_shp_count_file = stitched_ifg_dir / "shp_counts.tif"
+    stitching.merge_images(
+        shp_count_file_list,
+        outfile=stitched_shp_count_file,
+        driver="GTiff",
+        out_bounds=out_bounds,
+        out_bounds_epsg=output_options.bounds_epsg,
+    )
+
     if output_options.add_overviews:
         logger.info("Creating overviews for stitched images")
         create_overviews(stitched_ifg_paths, image_type=ImageType.INTERFEROGRAM)
@@ -124,6 +138,7 @@ def run(
         create_image_overviews(stitched_ps_file, image_type=ImageType.PS)
         create_image_overviews(stitched_temp_coh_file, image_type=ImageType.CORRELATION)
         create_image_overviews(stitched_amp_disp_file, image_type=ImageType.CORRELATION)
+        create_image_overviews(stitched_shp_count_file, image_type=ImageType.PS)
 
     return (
         stitched_ifg_paths,
@@ -131,4 +146,5 @@ def run(
         stitched_temp_coh_file,
         stitched_ps_file,
         stitched_amp_disp_file,
+        stitched_shp_count_file,
     )
diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py
index 4d1725a1..19721ad6 100644
--- a/src/dolphin/workflows/wrapped_phase.py
+++ b/src/dolphin/workflows/wrapped_phase.py
@@ -18,7 +18,7 @@
 @log_runtime
 def run(
     cfg: DisplacementWorkflow, debug: bool = False, tqdm_kwargs=None
-) -> tuple[list[Path], list[Path], Path, Path, Path]:
+) -> tuple[list[Path], list[Path], Path, Path, Path, Path]:
     """Run the displacement workflow on a stack of SLCs.
 
     Parameters
@@ -46,6 +46,8 @@ def run(
         The multilooked boolean persistent scatterer file.
     amp_disp_looked_file : Path
         The multilooked amplitude dispersion file.
+    shp_count_file : Path
+        Path to the created SHP counts file.
 
     """
     if tqdm_kwargs is None:
@@ -146,6 +148,7 @@ def run(
         logger.info(f"Skipping EVD step, {len(phase_linked_slcs)} files already exist")
         comp_slc_list = sorted(pl_path.glob("compressed*tif"))
         temp_coh_file = next(pl_path.glob("temporal_coherence*tif"))
+        shp_count_file = next(pl_path.glob("shp_count*tif"))
     else:
         logger.info(f"Running sequential EMI step in {pl_path}")
         kwargs = tqdm_kwargs | {"desc": f"Phase linking ({pl_path})"}
@@ -154,27 +157,25 @@ def run(
         # If we pre-compute it from some big stack, we need to use that for SHP
         # finding, not use the size of `slc_vrt_file`
         shp_nslc = None
-        (
-            phase_linked_slcs,
-            comp_slc_list,
-            temp_coh_file,
-        ) = sequential.run_wrapped_phase_sequential(
-            slc_vrt_file=vrt_stack.outfile,
-            ministack_planner=ministack_planner,
-            ministack_size=cfg.phase_linking.ministack_size,
-            half_window=cfg.phase_linking.half_window.model_dump(),
-            strides=strides,
-            use_evd=cfg.phase_linking.use_evd,
-            beta=cfg.phase_linking.beta,
-            mask_file=nodata_mask_file,
-            ps_mask_file=ps_output,
-            amp_mean_file=cfg.ps_options._amp_mean_file,
-            amp_dispersion_file=cfg.ps_options._amp_dispersion_file,
-            shp_method=cfg.phase_linking.shp_method,
-            shp_alpha=cfg.phase_linking.shp_alpha,
-            shp_nslc=shp_nslc,
-            block_shape=cfg.worker_settings.block_shape,
-            **kwargs,
+        (phase_linked_slcs, comp_slc_list, temp_coh_file, shp_count_file) = (
+            sequential.run_wrapped_phase_sequential(
+                slc_vrt_file=vrt_stack.outfile,
+                ministack_planner=ministack_planner,
+                ministack_size=cfg.phase_linking.ministack_size,
+                half_window=cfg.phase_linking.half_window.model_dump(),
+                strides=strides,
+                use_evd=cfg.phase_linking.use_evd,
+                beta=cfg.phase_linking.beta,
+                mask_file=nodata_mask_file,
+                ps_mask_file=ps_output,
+                amp_mean_file=cfg.ps_options._amp_mean_file,
+                amp_dispersion_file=cfg.ps_options._amp_dispersion_file,
+                shp_method=cfg.phase_linking.shp_method,
+                shp_alpha=cfg.phase_linking.shp_alpha,
+                shp_nslc=shp_nslc,
+                block_shape=cfg.worker_settings.block_shape,
+                **kwargs,
+            )
         )
 
     # ###################################################
@@ -191,6 +192,7 @@ def run(
             temp_coh_file,
             ps_looked_file,
             amp_disp_looked_file,
+            shp_count_file,
         )
 
     logger.info(f"Creating virtual interferograms from {len(phase_linked_slcs)} files")
@@ -203,6 +205,7 @@ def run(
         temp_coh_file,
         ps_looked_file,
         amp_disp_looked_file,
+        shp_count_file,
     )
 
 


Test patch:
diff --git a/tests/test_shp.py b/tests/test_shp.py
index 4850d99d..b3743459 100644
--- a/tests/test_shp.py
+++ b/tests/test_shp.py
@@ -9,7 +9,7 @@
 simulate._seed(1234)
 
 
-@pytest.mark.parametrize("method", ["glrt", "ks", "kld"])
+@pytest.mark.parametrize("method", ["glrt", "ks"])
 def test_shp_glrt_tf_smoketest(method):
     shape = (5, 50, 50)
     slcs = 20 * (np.random.rand(*shape) + 1j * np.random.rand(*shape))
@@ -29,27 +29,37 @@ def test_shp_glrt_tf_smoketest(method):
     )
 
 
+NUM_SLCS = 30
+
+
 @pytest.fixture(scope="module")
-def slcs(shape=(30, 11, 11)):
+def slcs(shape=(NUM_SLCS, 11, 11)):
     return 20 * (np.random.rand(*shape) + 1j * np.random.rand(*shape))
 
 
-# the KLD and GLRT methods have same API for estimate_neighbors
-@pytest.mark.parametrize("method", ["kld", "glrt"])
-def test_shp_glrt_kld(slcs, method):
+@pytest.fixture(scope="module")
+def mean(slcs):
     amp_stack = np.abs(slcs)
-    mean = np.mean(amp_stack, axis=0)
-    var = np.var(amp_stack, axis=0)
+    return np.mean(amp_stack, axis=0)
+
+
+@pytest.fixture(scope="module")
+def var(slcs):
+    amp_stack = np.abs(slcs)
+    return np.var(amp_stack, axis=0)
+
+
+def test_shp_glrt(mean, var):
+    method = "glrt"
 
     halfwin_rowcol = (5, 5)  # Looking at the entire stack
-    nslc = slcs.shape[0]
 
     # First try a small alpha
     neighbors = shp.estimate_neighbors(
         mean=mean,
         var=var,
         halfwin_rowcol=halfwin_rowcol,
-        nslc=nslc,
+        nslc=NUM_SLCS,
         alpha=0.005,
         method=method,
     )
@@ -59,26 +69,25 @@ def test_shp_glrt_kld(slcs, method):
 
     # Check the edges are cut off
     top_left = neighbors[0, 0]
-    assert top_left[:5, :5].all()
+    assert top_left[:5, :5].sum() == 5 * 5 - 1
     assert top_left[6:, :].sum() == top_left[:, 6:].sum() == 0
 
 
 def test_shp_ks(slcs):
     amp_stack = np.abs(slcs)
     halfwin_rowcol = (5, 5)  # Looking at the entire stack
-    nslc = slcs.shape[0]
 
     # First try a tiny alpha
     neighbors = shp.estimate_neighbors(
         amp_stack=amp_stack,
         halfwin_rowcol=halfwin_rowcol,
-        nslc=nslc,
+        nslc=NUM_SLCS,
         alpha=1e-6,
         method="ks",
     )
     shps_mid_pixel = neighbors[5, 5]
     # Check that everything is counted as a neighbor
-    assert shps_mid_pixel.all()
+    assert shps_mid_pixel.sum() == shps_mid_pixel.size - 1
 
     # Check the edges are cut off and all zeros
     # NOTE: this is different than TF test since KS is slower
@@ -90,22 +99,17 @@ def test_shp_ks(slcs):
     neighbors = shp.estimate_neighbors(
         amp_stack=amp_stack,
         halfwin_rowcol=halfwin_rowcol,
-        nslc=nslc,
+        nslc=NUM_SLCS,
         alpha=1.0,
         method="ks",
     )
     shps_mid_pixel = neighbors[5, 5]
-    assert shps_mid_pixel.sum() == 1  # only itself
-    assert shps_mid_pixel[5, 5] == 1
+    assert shps_mid_pixel.sum() == 0  # only itself
 
 
-@pytest.mark.parametrize("method", ["kld", "glrt"])
-def test_shp_half_mean_different(slcs, method):
+def test_shp_half_mean_different(mean, var):
     """Run a test where half the image has different mean"""
-    amp_stack = np.abs(slcs)
-    mean = np.mean(amp_stack, axis=0)
-    var = np.var(amp_stack, axis=0)
-    nslc = slcs.shape[0]
+    method = "glrt"
 
     halfwin_rowcol = (5, 5)
     # make the top half different amplitude
@@ -119,23 +123,20 @@ def test_shp_half_mean_different(slcs, method):
         mean=mean2,
         var=var,
         halfwin_rowcol=halfwin_rowcol,
-        nslc=nslc,
+        nslc=NUM_SLCS,
         alpha=0.01,
         method=method,
     )
     shps_mid_pixel = neighbors[5, 5]
     # Check that everything is counted as a neighbor
-    assert shps_mid_pixel[5:, :].all()
+    top_block = shps_mid_pixel[5:, :]
+    assert top_block.sum() == top_block.size - 1
     assert not shps_mid_pixel[:5, :].any()
 
 
-@pytest.mark.parametrize("method", ["kld", "glrt"])
-def test_shp_half_var_different(slcs, method):
+def test_shp_half_var_different(mean, var):
     """Run a test where half the image has different variance"""
-    amp_stack = np.abs(slcs)
-    mean = np.mean(amp_stack, axis=0)
-    var = np.var(amp_stack, axis=0)
-    nslc = slcs.shape[0]
+    method = "glrt"
 
     halfwin_rowcol = (5, 5)
     # make the top half different amplitude
@@ -149,17 +150,40 @@ def test_shp_half_var_different(slcs, method):
         mean=mean,
         var=var2,
         halfwin_rowcol=halfwin_rowcol,
-        nslc=nslc,
+        nslc=NUM_SLCS,
         alpha=0.01,
         method=method,
     )
     shps_mid_pixel = neighbors[5, 5]
-    # Check that everything is counted as a neighbor
-    assert shps_mid_pixel[5:, :].all()
+    top_block = shps_mid_pixel[5:, :]
+    assert top_block.sum() == top_block.size - 1
     assert not shps_mid_pixel[:5, :].any()
 
 
-@pytest.mark.parametrize("method", ["glrt", "kld", "ks"])
+@pytest.mark.parametrize("strides", [{"x": 1, "y": 1}, {"x": 2, "y": 2}])
+def test_shp_glrt_nodata_0(mean, var, strides):
+    """Ensure"""
+    method = "glrt"
+
+    halfwin_rowcol = (5, 5)  # Looking at the entire stack
+
+    mean = mean.copy()
+    mean[:4, :4] = 0
+    # First try a small alpha
+    neighbors = shp.estimate_neighbors(
+        mean=mean,
+        var=var,
+        halfwin_rowcol=halfwin_rowcol,
+        strides=strides,
+        nslc=NUM_SLCS,
+        alpha=0.005,
+        method=method,
+    )
+    out_col, out_row = 2 // strides["x"], 2 // strides["x"]
+    assert neighbors[:out_row, :out_col, :, :].sum() == 0
+
+
+@pytest.mark.parametrize("method", ["glrt", "ks"])
 @pytest.mark.parametrize("alpha", [0.01, 0.05])
 @pytest.mark.parametrize("strides", [{"x": 1, "y": 1}, {"x": 2, "y": 2}])
 def test_shp_statistics(method, alpha, strides):
@@ -175,13 +199,13 @@ def test_shp_statistics(method, alpha, strides):
         amp_stack = amp_stack_sims[i]
         mean = np.mean(amp_stack, axis=0)
         var = np.var(amp_stack, axis=0)
-        nslc = amp_stack.shape[0]
+        amp_stack.shape[0]
 
         neighbors = shp.estimate_neighbors(
             mean=mean,
             var=var,
             halfwin_rowcol=halfwin_rowcol,
-            nslc=nslc,
+            nslc=NUM_SLCS,
             strides=strides,
             alpha=alpha,
             method=method,