initial implementation of vce

depsi/arc_estimation.py

@@ -1438,7 +1438,7 @@ def periodogram(
         Key for the temporal baseline in the STM.
         The value should be in decimal years.
     std_obs : float, optional
-        A-poriori standard deviation of the observations in rads, by default 1.0.
+        A-priori standard deviation of the observations in rads, by default 1.0.
         This value is used to construct the stochastic model (Qyy) of the observations.
     std_height : float, optional
         A-priori standard deviation of the height in meters, by default 50.0.

depsi/constants.py

@@ -6,3 +6,6 @@
 
 # Speed of light in vacuum
 SPEED_OF_LIGHT = 299792458.0  # m/s
+
+# Sentinel-1 wavelength
+WAVELENGTH_S1 = 0.055465763  # m

depsi/network.py

@@ -2007,3 +2007,48 @@ def _network_relation_matrix(idx_source, idx_target, n_points, idx_refpnt, spars
         A = np.array(A.todense())
 
     return A
+
+
+def _independent_arcs(arcs: np.ndarray) -> np.ndarray:
+    """Select independent arcs from a list of arcs.
+
+    An arc is independent if its starting and ending points do not exist in any other arc's
+    starting or ending points.
+
+    Parameters
+    ----------
+    arcs : np.ndarray
+        A 2D array of shape (n_points, 2) where each row represents indices of the starting and ending points
+        of an arc.
+
+    Returns
+    -------
+    np.ndarray
+        A 2D array of independent arcs, where each row represents indices of the starting and ending points
+        of an arc.
+    """
+    # Select arcs with unique starting points
+    _, unique_idx_start = np.unique(arcs[:, 0], return_index=True)
+    arcs = arcs[unique_idx_start, :]
+
+    # Select arcs with unique ending points
+    _, unique_idx_end = np.unique(arcs[:, 1], return_index=True)
+    arcs = arcs[unique_idx_end, :]
+
+    # After previous two steps, no arcs will share starting or ending points.
+    # However, there starting points may be the ending points of other arcs, and vice versa.
+    # To ensure independency, we loop through the rest arcs and add arc one by one
+    # In each interation, remove arcs that
+    # 1) start with the ending point of this arc, or
+    # 2) end with the starting point of this arc
+    arcs_selected = np.empty((0, 2), dtype=int)
+    while arcs.shape[0] > 0:
+        arc_current = arcs[0, :]
+        arcs_selected = np.append(arcs_selected, [arc_current], axis=0)
+        # Remove arcs which contain the starting point or ending point of the current arc
+        idx_remove = np.where((arcs[:, 1] == arc_current[0]) | (arcs[:, 0] == arc_current[1]))[0]
+        # add the index of the current arc to idx_remove
+        idx_remove = np.append(idx_remove, 0)
+        arcs = np.delete(arcs, idx_remove, axis=0)
+
+    return arcs_selected

depsi/stochastic.py

@@ -0,0 +1,160 @@
+"""stochastic model related functions."""
+
+import numpy as np
+import xarray as xr
+
+from depsi.arc_estimation import periodogram
+from depsi.network import _independent_arcs, form_network
+from depsi.utils import get_m2ph
+
+SIGMA_MOTHER_ATMO = 15.0  # std of mother atmosphere when its included in the stochastic model
+SIGMA_OTHERS = 20.0  # std of other interferograms when master
+SIGMA_OVERALL = 30  # std when mother atmosphere is in functional model instead of stochastic model
+SIGMA_MINIMUM = 10  # minimum value for variance components, used to avoid negative values, in degrees
+
+
+def vce_temporal(
+    stm: xr.Dataset,
+    key_phase: str,
+    key_Btemporal: str,
+    key_h2ph: str,
+    key_x="lon",
+    key_y="lat",
+    max_length: float = 0.01,  # TODO: change distance from degree to meters
+    include_mother_atmo: bool = False,
+) -> np.ndarray:
+    """Estimate variance components per epoch.
+
+    This estimation is performed on an Space-Time Matrix (STM) of points.
+    Independent arcs are formed from the STM and unwrapped as redundancies of the estimation.
+
+    Parameters
+    ----------
+    stm : xr.Dataset
+        Input space-time matrix (STM).
+    key_phase : str
+        key for the phase data variable in the STM.
+    key_Btemporal : str
+        key for the Btemp data variable in the STM.
+    key_h2ph : str
+        key for the h2ph data variable in the STM.
+    key_x : str, optional
+        x coordinate for network formation, by default "lon"
+    key_y : str, optional
+        y coordinate for network formation, by default "lat"
+    max_length : float, optional
+        maximum length of the arcs, in degrees, by default 0.01
+    include_mother_atmo : bool, optional
+        whether to include mother atmosphere in the stochastic model.
+        when False, it is assumed that the mother atmosphere is included in the functional model,
+        by default False
+
+    Returns
+    -------
+    np.ndarray
+        Estimated variance components per epoch, in radians squared.
+        This is an array of shape (Nifgs + 1,) where Nifgs is the number of interferograms.
+        The extra element is for the mother atmosphere.
+        The first element is the variance component for the mother atmosphere.
+        If `include_mother_atmo` is False, the first element is 0.0.
+    """
+    # Generate a delaunay network
+    arcs = form_network(
+        stm, key_phase=key_phase, key_h2ph=key_h2ph, key_Btemporal=key_Btemporal, network_method="delaunay"
+    )
+
+    # Select independent arcs
+    arcs_source_target = _independent_arcs(np.stack([arcs["source"].data, arcs["target"].data], axis=1))
+    arcs_source_target_set = set(map(tuple, arcs_source_target))
+    pairs = np.stack([arcs["source"].data, arcs["target"].data], axis=1)
+    mask = np.array([tuple(x) in arcs_source_target_set for x in pairs])
+    arcs = arcs.isel(space=mask)
+
+    # Unwrap the arcs, arcs stm has standard data vars 'd_phase', 'h2ph', 'Btemp'
+    phase_unwrapped, _, _, _, _ = periodogram(arcs, key_dphase="d_phase", key_h2ph="h2ph", key_Btemporal="Btemp")
+
+    # Intiate variance components
+    Nifgs = stm.sizes["time"]
+    if include_mother_atmo:
+        # include mother atmosphere in the stochastic model
+        Qy1, Qy = _q_with_mother_atmo(Nifgs)
+    else:
+        # estimate variance components in the functional model
+        Qy1, Qy = _q_no_mother_atmo(Nifgs)
+
+    Qyinv = np.linalg.inv(Qy)
+
+    # Compute Pao and QP
+    Btemp = stm[key_Btemporal].values
+    h2ph_approx = stm[key_h2ph].mean(dim="space").values  # Mean h2ph of all arcs
+    m2ph = get_m2ph()  # Convert meters to phase
+    B = np.stack([h2ph_approx * m2ph, Btemp * m2ph]).T  # Design matrix of the functional model
+    Pao = np.eye(Nifgs) - B @ np.linalg.inv(B.T @ Qyinv @ B) @ B.T @ Qyinv
+    QP = Qyinv @ Pao
+
+    # Compute QPQy1QP and N, optimized using einsum
+    # The following code is equivalent to this nested loop:
+    # Nsig = Qy1.shape[2]  # Number of sigmas, i.e. the number of components to estimate
+    # Narcs_vce = phase_unwrapped.shape[0]  # Number of arcs for VCE
+    # QPQy1QP = np.full((Nifgs, Nifgs, Nsig), np.nan)
+    # N = np.full((Nsig, Nsig), np.nan)
+    # for k in range(Nsig):
+    #     QPQy1QP[:, :, k] = QP @ Qy1[:, :, k] @ QP
+    #     for j in range(Nsig):
+    #         N[k, j] = np.trace(QPQy1QP[:, :, k] @ Qy1[:, :, j])
+    QPQy1QP = np.einsum("ij,jlk,lm ->imk", QP, Qy1, QP, optimize=True)
+    N = np.einsum("abk,baj->kj", QPQy1QP, Qy1, optimize=True)
+
+    Ninv = np.linalg.inv(N)
+
+    # Estimate variance components from all independent arcs
+    # The following code is equivalent to this nested loop:
+    # sig2 = np.full((Nsig, Narcs_vce), np.nan)
+    # l = np.full((Nsig, 1), np.nan)
+    # for v in range(Narcs_vce):
+    #     y = phase_unwrapped[v, :].reshape(-1, 1)
+    #     for k in range(Nsig):
+    #         l[k, 0] = (y.T @ QPQy1QP[:, :, k] @ y).squeeze()
+    #     sig2[:, v] = (Ninv @ l).flatten()
+    l_vec = np.einsum("ij,jmk,mi ->ki", phase_unwrapped, QPQy1QP, phase_unwrapped.T, optimize=True)  # l vector
+    sig2_all_arcs = Ninv @ l_vec
+    sig2_est = np.mean(sig2_all_arcs, axis=1)
+
+    # Apply threshold to avoid small and negative values
+    threshold = (np.pi * SIGMA_MINIMUM / 180) ** 2
+    sig2_est[sig2_est < threshold] = threshold
+
+    # if atmosphere is not included in the stochastic model, add a zero at the beginning
+    if not include_mother_atmo:
+        sig2_est = np.insert(sig2_est, 0, 0.0)
+
+    return sig2_est
+
+
+def _q_with_mother_atmo(Nifgs: int) -> tuple:
+    """Build Qy1 and Qy matrices with mother atmosphere."""
+    sig0_mother = (np.pi * SIGMA_MOTHER_ATMO / 180) ** 2
+    sig0_others = (np.pi * SIGMA_OTHERS / 180) ** 2
+
+    # Build "Design matrix" for variance components estimation
+    Qy1 = np.zeros((Nifgs, Nifgs, Nifgs + 1), dtype=np.int16)  # Build Qy1 matrix, has 0 or 2
+    Qy1[:, :, 0] = 2  # Mother epoch, all 2
+    for v in range(0, Nifgs):
+        Qy1[v, v, v + 1] = 2  # Other epochs, on location set 2
+    Qy = np.full((Nifgs, Nifgs), 2 * sig0_mother)  # Build Qy matrix has 2 * sig0_mother in background
+    for v in range(Nifgs):
+        Qy[v, v] += 2 * sig0_others  # Per epoch add 2 * sig0_others
+
+    return Qy1, Qy
+
+
+def _q_no_mother_atmo(Nifgs: int) -> tuple:
+    """Build Qy1 and Qy matrices without mother atmosphere."""
+    sig0 = (np.pi * SIGMA_OVERALL / 180) ** 2
+
+    # Build "Design matrix" for variance components estimation
+    Qy1 = np.zeros((Nifgs, Nifgs, Nifgs), dtype=np.int16)  # Build Qy1 matrix, has 0 or 2
+    for v in range(Nifgs):
+        Qy1[v, v, v] = 2  # Other epochs, on location set 2
+    Qy = np.diag(np.full(Nifgs, 2 * sig0))  # Build Qy matrix has 2 * sig0 in background
+    return Qy1, Qy

depsi/utils.py

@@ -31,7 +31,7 @@
 import pytz
 import xarray as xr
 
-from depsi.constants import EARTH_RADIUS
+from depsi.constants import EARTH_RADIUS, WAVELENGTH_S1
 
 logger = logging.getLogger(__name__)
 
@@ -834,3 +834,8 @@ def concatenate_stms(
     stm_dens_pnts_output = stm_dens_pnts_output.reset_coords(names=coords_to_reset, drop=False)
 
     return stm_dens_pnts_output
+
+
+def get_m2ph(wavelength: float = WAVELENGTH_S1):
+    """Get the conversion factor from meters to phase."""
+    return -4 * np.pi / wavelength

tests/test_network.py

@@ -7,6 +7,7 @@
 from depsi.network import (
     _ensure_network_min_connections,
     _ensure_single_network,
+    _independent_arcs,
     _network_relation_matrix,
     _remove_network_points_min_connections,
     form_network,
@@ -500,3 +501,31 @@ def test_init_network_relation_matrix_sparse(
 
         assert A.shape == A_exp.shape
         assert np.all(A.todense() == A_exp)
+
+
+class TestArcsUtils:
+    @pytest.mark.timeout(10)  # Each should finish in 10 seconds
+    @pytest.mark.parametrize(
+        "npoints, narcs",
+        [
+            (103, 1000),
+            (1923, 10000),
+            (12, 30),
+        ],
+    )
+    def test_independent_arcs(self, npoints, narcs):
+        # Simulate random arcs
+        rng = np.random.default_rng(42)
+        arcs = rng.integers(0, npoints, size=(narcs, 2))
+        # remove arcs which has the same start and end point
+        arcs = arcs[arcs[:, 0] != arcs[:, 1]]
+        # Remove duplicate arcs
+        arcs = np.unique(np.sort(arcs, axis=1), axis=0)
+
+        # Test that the arcs are independent.
+        independent_arcs = _independent_arcs(arcs)
+
+        # A point index should only appear once
+        # either as a start or end point of an arc.
+        all_idx = independent_arcs.flatten()
+        assert all_idx.shape == np.unique(all_idx).shape

tests/test_stochastic.py

@@ -0,0 +1,81 @@
+import numpy as np
+import pytest
+import xarray as xr
+
+from depsi.stochastic import _q_no_mother_atmo, _q_with_mother_atmo, vce_temporal
+
+
+def simulated_stm(n_ifg, n_points):
+    """Function to simulate arcs stm for testing."""
+    rng = np.random.default_rng(31)
+
+    # Simulate coordinates
+    lat = np.linspace(51.14, 51.15, n_points)
+    lon = rng.uniform(6.9, 7.0, n_points)
+
+    # Simulate Btemp
+    TIME_STEP = 15  # Time step in days
+    years = np.linspace(0, (n_ifg - 1) * TIME_STEP / 365.25, n_ifg)  # years from 0 to n_ifg-1
+
+    # Simulated phase values, linear+ noise
+    phase = np.tile(np.linspace(-np.pi, np.pi, n_ifg), (n_points, 1)) + rng.normal(0, 0.1, (n_points, n_ifg))
+
+    # Simulated h2ph values
+    h2ph = rng.random((n_points, n_ifg)) * 1e-3  # fixed h2ph values for each arc
+
+    arcs = xr.Dataset(
+        data_vars={
+            "lon": (("space",), lon),
+            "lat": (("space",), lat),
+            "phase": (("space", "time"), phase),
+            "h2ph": (("space", "time"), h2ph),
+            "Btemporal": (("time",), years),
+        }
+    )
+
+    arcs.attrs["wavelength"] = 0.056  # example wavelength in meters
+
+    return arcs
+
+
+@pytest.mark.parametrize(
+    ["n_ifg", "n_points", "include_mother_atmo"],
+    [
+        (12, 41, False),
+        (19, 107, False),
+        (42, 127, True),
+    ],
+)
+def test_vce_temporal(n_ifg, n_points, include_mother_atmo):
+    arcs = simulated_stm(n_ifg, n_points)
+    sigma2 = vce_temporal(
+        arcs, key_phase="phase", key_Btemporal="Btemporal", key_h2ph="h2ph", include_mother_atmo=include_mother_atmo
+    )
+
+    assert sigma2.shape[0] == n_ifg + 1
+    if include_mother_atmo:
+        assert sigma2[0] != 0
+    else:
+        assert sigma2[0] == 0
+
+
+@pytest.mark.parametrize("Nifgs", [4, 11, 23])
+def test_q_no_mother_atmo(Nifgs):
+    """Test the _q_with_mother_atmo function."""
+    Qy1, Qy = _q_no_mother_atmo(Nifgs)
+
+    assert Qy1.shape == (Nifgs, Nifgs, Nifgs)
+    assert Qy.shape == (Nifgs, Nifgs)
+    assert np.all(np.isin(np.unique(Qy1.flatten()), [0, 2]))
+    assert np.all(Qy == np.diag(np.diag(Qy)))  # check Qy is a diagonal matrix
+
+
+@pytest.mark.parametrize("Nifgs", [4, 11, 23])
+def test_q_with_mother_atmo(Nifgs):
+    """Test the _q_with_mother_atmo function."""
+    Qy1, Qy = _q_with_mother_atmo(Nifgs)
+
+    assert Qy1.shape == (Nifgs, Nifgs, Nifgs + 1)
+    assert Qy.shape == (Nifgs, Nifgs)
+    assert np.all(np.isin(np.unique(Qy1.flatten()), [0, 2]))
+    assert np.all(np.isin(np.unique(Qy1[:, :, 0].flatten()), [2]))  # check first slice is all 2