Profiling example

profiling_example/README.md

@@ -0,0 +1,55 @@
+# Profiling example for Network unwrapping
+
+This is an example of how to perform profiling on the Arc Unwrapping algorithms. In this example we used the so-called "LAMBDA method" developed by TUDelft, and use a wrapper function in DePSI to call the unwrapping algorithm.
+
+- LAMBDA method: [LAMBDA.py](https://github.com/TUDelftGeodesy/DePSI_group/blob/dev/depsi/LAMBDA.py) (private repo)
+- Wrapper funtion `lambda_estimation`: [lambda_estimation](https://github.com/TUDelftGeodesy/DePSI_group/blob/dev/depsi/LAMBDA.py) (private repo)
+- Data needed for the example: [stm_amsterdam_173p.zarr](https://zenodo.org/records/15324181/files/stm_amsterdam_173p.zarr.zip?download=1)
+
+
+## File structure
+
+- lambda_unwrap.ipynb: Jupyter notebook for debugging and visualization
+- lambda_unwrap.py: Python script for running the unwrapping algorithm in a recursive loop
+- lambda_unwrap_dask.ipynb: Jupyter notebook for debugging dask method
+- lambda_unwrap_dask.py: Python script for running the same unwrapping algorithm with dask databags
+
+## Profiling the unwrapping algorithm
+
+We use the `py-spy` package to profile the unwrapping algorithm. For example, we can use the following command to profile the unwrapping algorithm without dask
+
+### Loop method
+
+```sh
+py-spy record --output profile_loop_60pnts --idle --rate 5 --subprocesses --format speedscope python lambda_unwrap.py
+```
+
+### Dask with `processes` schedular:
+
+In `lambda_unwrap_dask.py` configure: 
+
+```py
+# Configure dask scheduler
+dask.config.set(scheduler="processes") 
+```
+
+```sh
+py-spy record --output profile_dask_60pnts_processes --idle --rate 5 --subprocesses --format speedscope python lambda_unwrap_dask.py
+```
+
+### Dask with `threads` schedular:
+
+```py
+# Configure dask scheduler
+dask.config.set(scheduler="threads") 
+```
+
+```sh
+py-spy record --output profile_dask_60pnts_threads --idle --rate 5 --subprocesses --format speedscope python lambda_unwrap_dask.py
+```
+
+Then you can visualize the profile using the [`speedscope` web tool](https://www.speedscope.app/)
+
+## Profiling results
+
+The results for 5 points, 24 points and 60 points can be found at: [this link](https://zenodo.org/records/15393928/files/experiments.zip)

profiling_example/lambda.ipynb

@@ -0,0 +1,317 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# An example to generate STM from dynamic estimation\n",
+    "\n",
+    "adapted from https://github.com/TUDelftGeodesy/DePSI_group/blob/dev/examples/notebooks/demo_dynamic_estimation.ipynb"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# from depsi.LAMBDA import *\n",
+    "from depsi.dynamic_estimation import lambda_estimation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Set input parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# set wavelength of Sentinel-1 (C-band)\n",
+    "wavelength = 0.055465763 # m\n",
+    "init_len = 50\n",
+    "sigma_acc = 0.003 # m/yr^2\n",
+    "L = 30/365\n",
+    "\n",
+    "# Initial gauss for the sigma of the unknown parameters\n",
+    "sigma_offset = 0.001 #m\n",
+    "sigma_vel =  0.0001 #m/yr\n",
+    "sigma_h = 5 #m\n",
+    "sigma_ther = 0.00005 #m/°C\n",
+    "\n",
+    "# the option for the LAMBDA method\n",
+    "# method == 1: ILS with shrinking search\n",
+    "# method == 2: Integer rounding\n",
+    "# method == 3: Integer bootstrapping\n",
+    "# method == 4: PAR\n",
+    "# method == 5: ILS with Ratio Test\n",
+    "method = 3"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Load the STM derived from PS selection"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load all to memory\n",
+    "stm = xr.open_zarr('../../data/stm_amsterdam_173p.zarr').compute()\n",
+    "stm"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Choose the reference point"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Select the point with the smallest NMAD for initialization as the reference point\n",
+    "nmad_init = np.array(stm['nmad_init'])\n",
+    "ref_pnt_idx= int(np.argmin(nmad_init))\n",
+    "ref_pnt_idx"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Define functions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def polynomial_fitting(x, y, degree=1):\n",
+    "    \"\"\"\n",
+    "    Perform polynomial fitting on the input data (x, y).\n",
+    "\n",
+    "    Parameters:\n",
+    "    \n",
+    "    x : array_like\n",
+    "        Independent variable data points.\n",
+    "    y : array_like \n",
+    "        Dependent variable data points.\n",
+    "    degree: int\n",
+    "        Degree of the polynomial fit (default is 1 for linear).\n",
+    "\n",
+    "    Returns:\n",
+    "    y_fitted : ndarray\n",
+    "        Fitted y values for the given x based on the polynomial fit.\n",
+    "    poly_fn : np.poly1d\n",
+    "        Polynomial function that represents the fit.\n",
+    "    coeffs : ndarray \n",
+    "        Coefficients of the fitted polynomial.\n",
+    "    \"\"\"\n",
+    "    # Validate inputs\n",
+    "    if len(x) != len(y):\n",
+    "        raise ValueError(\"Input arrays 'x' and 'y' must have the same length.\")\n",
+    "    if degree < 1:\n",
+    "        raise ValueError(\"Degree must be at least 1.\")\n",
+    "\n",
+    "    # Perform polynomial fitting\n",
+    "    coeffs = np.polyfit(x, y, degree)\n",
+    "    poly_fn = np.poly1d(coeffs)\n",
+    "    y_fitted = poly_fn(x)\n",
+    "\n",
+    "    return y_fitted, poly_fn, coeffs\n",
+    "\n",
+    "\n",
+    "\n",
+    "def NMAD_to_sigma_phase(nmad, method):\n",
+    "\n",
+    "    \"\"\"\n",
+    "    Converts Normalized Median Absolute Deviation (NMAD) to the sigma of phase observations.\n",
+    "\n",
+    "    This function uses an empirical cubic approximation (and is not based on physics):\n",
+    "    sigma = a + b*nmad + c*nmad**2 + d*nmad**3\n",
+    "\n",
+    "    Parameters:\n",
+    "    nmad : array_like\n",
+    "        Normalized Median Absolute Deviation value.\n",
+    "    method : {'mean', 'mean_2_sigma'}\n",
+    "        Method used to compute sigma.\n",
+    "    \n",
+    "    Returns:\n",
+    "    sigma: ndarray\n",
+    "        Estimated sigma value.\n",
+    "    \"\"\"\n",
+    "    if method == 'mean':\n",
+    "        a, b, c, d = -0.0144869469, 2.00028682, -5.23271341, 21.1111801\n",
+    "    elif method == 'mean_2_sigma':\n",
+    "        a, b, c, d = 0.01907808, 1.2852969, 1.90052824, 11.60677721\n",
+    "    else:\n",
+    "        raise ValueError(\"Invalid method. Choose 'mean' or 'mean_2_sigma'.\")\n",
+    "\n",
+    "    sigma = a + b*nmad + c*nmad**2 + d*nmad**3\n",
+    "    \n",
+    "    return sigma"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Perform unwrapping with lambda estimation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# select the reference point\n",
+    "stm_refpnt = stm.isel(space=ref_pnt_idx)\n",
+    "\n",
+    "# subset the data for debug\n",
+    "NUM_POINTS = 5\n",
+    "stm = stm.isel(space=range(5))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Initiate empty arrays to store results\n",
+    "x_hat = np.zeros((stm.sizes[\"space\"], 4))\n",
+    "Q_xhat = np.zeros(( 4, 4, stm.sizes[\"space\"]))\n",
+    "y_hats = np.zeros((stm.sizes[\"space\"], stm.sizes[\"time\"]))\n",
+    "phs_unw_init = np.zeros((stm.sizes[\"space\"], stm.sizes[\"time\"]))\n",
+    "\n",
+    "\n",
+    "# Loop over all points and perform phase unwrapping\n",
+    "for pnt_id in range(stm.sizes[\"space\"]):\n",
+    "    # Select current point\n",
+    "    stm_1pnt = stm.isel(space=pnt_id)\n",
+    "\n",
+    "    print(f\"{pnt_id}/{stm.sizes['space']}\")\n",
+    "\n",
+    "    # Get the sigma for the arc with the NMAD from the incremental time series\n",
+    "    sigma_nmad_inc_i = NMAD_to_sigma_phase(\n",
+    "        stm_1pnt[\"nmad_inc_stm\"].data, \"mean_2_sigma\"\n",
+    "    )\n",
+    "    sigma_nmad_inc_j = NMAD_to_sigma_phase(\n",
+    "        stm_refpnt[\"nmad_inc_stm\"].data, \"mean_2_sigma\"\n",
+    "    )\n",
+    "    sigma_nmad_inc_arc = np.sqrt(\n",
+    "        np.square(sigma_nmad_inc_i) + np.square(sigma_nmad_inc_j)\n",
+    "    )\n",
+    "\n",
+    "    # Compute arc phase\n",
+    "    phs_wrapped_arc = np.angle(\n",
+    "        stm_refpnt[\"sd_complex\"].data * stm_1pnt[\"sd_complex\"].data.conj()\n",
+    "    )\n",
+    "\n",
+    "    # Compute 'mean' h2ph value for the arc (which we currently model as the average of the two time series)\n",
+    "    h2ph_arc = (stm_refpnt[\"h2ph_values\"].data + stm_1pnt[\"h2ph_values\"].data) / 2\n",
+    "\n",
+    "    results = lambda_estimation(\n",
+    "        wavelength=wavelength,\n",
+    "        phs_wrapped=phs_wrapped_arc,\n",
+    "        sigma_phs_apri=sigma_nmad_inc_arc,\n",
+    "        years=stm_1pnt[\"years\"].data,\n",
+    "        h2ph_arc=h2ph_arc,\n",
+    "        temp=stm_1pnt[\"temperature\"].data,\n",
+    "        sigma_offset=sigma_offset,\n",
+    "        sigma_vel=sigma_vel,\n",
+    "        sigma_h=sigma_h,\n",
+    "        sigma_ther=sigma_ther,\n",
+    "        method=method,\n",
+    "    )\n",
+    "\n",
+    "    # Store results\n",
+    "    x_hat[pnt_id, :] = results[0]\n",
+    "    Q_xhat[:, :, pnt_id] = results[1]\n",
+    "    y_hats[pnt_id, :] = results[2]\n",
+    "    phs_unw_init[pnt_id, :] = results[3]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# For now, write x_hat and Q_xhat to npz files, since they so not fit in dimensions of an STM\n",
+    "np.savez(\n",
+    "    \"./x_hat_Q_xhat.npz\",\n",
+    "    x_hat=x_hat,\n",
+    "    Q_xhat=Q_xhat,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Attach y_hats and phs_unw_init to the STM\n",
+    "stm[\"y_hats\"] = ((\"space\", \"time\"), y_hats)\n",
+    "stm[\"phs_unw_init\"] = ((\"space\", \"time\"), phs_unw_init)\n",
+    "stm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stm.to_zarr(\n",
+    "    \"./stm_amsterdam_173p_init_unw.zarr\",\n",
+    "    mode=\"w\",\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "mobyle",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}