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

Issue:
Progress bar of wrapped phase is broken

E.g.

![Image](https://user-images.githubusercontent.com/8291800/216123629-17ef6ac2-c869-4723-a919-00c5d67be194.png)


Ignoring the `mean of empty slice` error, it ended at 86/195 progress. 
This means that this part is probably wrong: https://github.com/opera-adt/dolphin/blob/cc74bd516e0cb2ba94bbea9e277ce74fc15fc447/src/dolphin/sequential.py#L151-L153
[Bug]: Maximum-brightness selection for PS in the strided output is broken

The phase linking result with no ps_mask looks like this for one day:
<img width="604" alt="image" src="https://user-images.githubusercontent.com/8291800/222523709-34fb0dd5-c19c-49cc-b74c-a0fcb8b9aa7f.png">

The code where the PS pixels are filled by taking the brightest magnitude pixel in each stride window has bad artifacts:
<img width="753" alt="image" src="https://user-images.githubusercontent.com/8291800/222523891-aff92545-d6d5-4ef2-9ba2-edda90472785.png">

For a temporary fix, the new PR will just average the phase of PS pixels within the stride window (which is much easier to get the implementation right):
<img width="695" alt="image" src="https://user-images.githubusercontent.com/8291800/222524475-f624925f-2291-4f3d-8712-dc605640d974.png">


Two side notes

1. the PS mask, even at 0.25 $D_A$ threshold, looks like this, where most of Mauna Loa is  a PS:
<img width="534" alt="image" src="https://user-images.githubusercontent.com/8291800/222524262-4b3922aa-b6ed-486c-bbd2-1b30ba524e6d.png">
2. It looks like the PS pixels are still somewhat noisy... (zooming in)
<img width="374" alt="image" src="https://user-images.githubusercontent.com/8291800/222524896-42953f76-5ce9-433c-964d-e369e538de96.png">
this might require the addition of Sara's PS pruning methods, even at the lower 0.25 threshold.




Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index fabbd7f9..dcca90b5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,68 @@
-# Unreleased, 0.0.2
+# Unreleased
+
+# v0.0.4
+
+**Added**
+
+- Created first version of the single-update workflow, usable with `dolphin config --single`
+- `_background.py` module as the abstract classes for background input/output with `EagerLoader` and `Writer`.
+- `from_vrt_file` for the `VRTInterferogram` class.
+- Arbitrary interferogram index selection in `Network` class.
+- Parallel CPU eigenvector finding using `scipy.linalg.eigh`.
+- PS selection for strided outputs using the average PS phase within a window (that contains multiple PS).
+- Comments in the YAML file output by the `dolphin config` command.
+
+
+**Changed**
+
+- The main workflow has been renamed to `s1_disp.py` to better reflect the workflow, since it can handle both single and stack workflows.
+    - The `sequential.py` and `single.py` are where these differences are handled.
+- More uniform naming in `io.get_raster_<X>` functions.
+- The SLC compression is now done in `_compress.py` to declutter the `mle.py` module.
+- Replace `tqdm` with `rich` for progress bars.
+
+- Docs are now using the mkdocs `material` theme.
+
+**Removed**
+
+- `utils.parse_slc_strings` in favor of always using `utils.get_dates`.
+- `io.get_stack_nodata_mask`. This will be done using the nodata polygon, or not at all.
+
+
+**Dependencies**
+
+Added requirements:
+
+- pyproj>=3.2
+- rich>=12.0
+- threadpoolctl>=3.0
+
+For docs:
+- mkdocs-material
+- pymdown-extensions
+
+Removed requirements:
+
+- tqdm
+
+
+# 0.0.3
+
+**Added**
+
+- Ability for `VRTStack` to handle HDF5 files with subdatasets.
+    - The OPERA specific HDF5 files are now supported without extra configuration.
+- First stitching of interferograms in `stitching.py`.
+    - Users can pass multiple SLC burst (like COMPASS bursts) per date, and the workflow will process per stack then stitch per date.
+- More features for `load_gdal` to load in blocks.
+
+**Changed**
+
+- A small amount of regularization on the coherence matrix is done before inversion during phase linking to avoid singular matrices.
+- Renamed module to `_log.py`
+- `workflows/wrapped_phase.py` absorbed much logic formerly in `s1_disp_stack.py`.
+
+# 0.0.2
 
 **Added**
 
@@ -14,15 +78,6 @@
     - `vrt.py`
 - Created the `phase_link` subpackage for wrapped phase estimation.
 
-**Changed**
-
-**Deprecated**
-
-**Removed**
-
-**Fixed**
-
-**Dependencies**
 
 Added requirements:
 
@@ -30,7 +85,7 @@ Added requirements:
 - tqdm>=4.60
 
 
-# Unreleased, 0.0.1
+# 0.0.1
 
 **Added**
 
@@ -38,13 +93,6 @@ Added requirements:
 - Command line interface for running the workflows
 - Outline of project structure and utilities
 
-**Changed**
-
-**Deprecated**
-
-**Removed**
-
-**Fixed**
 
 **Dependencies**
 
diff --git a/README.md b/README.md
index cc50d8c3..1215efe5 100644
--- a/README.md
+++ b/README.md
@@ -43,8 +43,9 @@ Example usage:
 
 ```bash
 $ dolphin config --slc-files /path/to/slcs/*tif
+$ dolphin run dolphin_config.yaml
 ```
-This will create a YAML file (by default `dolphin_config.yaml` in the current directory).
+The `config` command creates a YAML file (by default `dolphin_config.yaml` in the current directory).
 
 The only required inputs for the workflow is a list of coregistered SLC files (in either geographic or radar coordinates).
 If the SLC files are spread over multiple files, you can either
diff --git a/conda-env-gpu-extras.yml b/conda-env-gpu-extras.yml
index fd6286af..c55ba9cd 100644
--- a/conda-env-gpu-extras.yml
+++ b/conda-env-gpu-extras.yml
@@ -9,3 +9,4 @@ channels:
 dependencies:
   - cupy>=11.0
   - cudatoolkit
+  - pynvml>=11.0
diff --git a/conda-env.yml b/conda-env.yml
index 8af06247..68e3d241 100644
--- a/conda-env.yml
+++ b/conda-env.yml
@@ -12,6 +12,7 @@ dependencies:
   - pydantic>=1.10
   - pymp-pypi>=0.4.5
   - pyproj>=3.2
+  - rich>=12.0
   - ruamel.yaml>=0.15
   - scipy>=1.5 # "scipy 0.16+ is required for linear algebra", numba. 1.5 is the oldest version that supports Python 3.7
-  - tqdm>=4.60
+  - threadpoolctl>=3.0
diff --git a/docker/create-lockfile.sh b/docker/create-lockfile.sh
index d642c0d5..643a9425 100755
--- a/docker/create-lockfile.sh
+++ b/docker/create-lockfile.sh
@@ -24,7 +24,7 @@ main() {
 
     # Get concretized package list.
     local PKGLIST
-    PKGLIST=$(docker run \
+    PKGLIST=$(docker run --network=host \
         -v "$ENVFILE:/tmp/environment.yml:ro" --rm \
         mambaorg/micromamba:1.1.0 bash -c '\
             micromamba install -y -n base -f /tmp/environment.yml > /dev/null && \
diff --git a/docker/specfile.txt b/docker/specfile.txt
index dd54b5c7..7b5e7926 100644
--- a/docker/specfile.txt
+++ b/docker/specfile.txt
@@ -11,18 +11,18 @@ https://conda.anaconda.org/conda-forge/linux-64/cairo-1.16.0-ha61ee94_1014.tar.b
 https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.18.1-h7f98852_0.tar.bz2#f26ef8098fab1f719c91eb760d63381a
 https://conda.anaconda.org/conda-forge/linux-64/cffi-1.15.1-py310h255011f_3.conda#800596144bb613cd7ac58b80900ce835
 https://conda.anaconda.org/conda-forge/linux-64/cfitsio-4.2.0-hd9d235c_0.conda#8c57a9adbafd87f5eff842abde599cb4
-https://conda.anaconda.org/conda-forge/linux-64/cryptography-39.0.0-py310h34c0648_0.conda#af4b0c22dc4006ce3c095e840cb2efd7
-https://conda.anaconda.org/conda-forge/linux-64/curl-7.87.0-hdc1c0ab_0.conda#b14123ca479b9473d7f7395b0fd25c97
+https://conda.anaconda.org/conda-forge/linux-64/cryptography-39.0.2-py310h34c0648_0.conda#99dc5a02a8b16cd88ca9a12354496e78
+https://conda.anaconda.org/conda-forge/linux-64/curl-7.88.1-hdc1c0ab_0.conda#1968e4fef727858ac04746560e820928
 https://conda.anaconda.org/conda-forge/linux-64/expat-2.5.0-h27087fc_0.tar.bz2#c4fbad8d4bddeb3c085f18cbf97fbfad
-https://conda.anaconda.org/conda-forge/linux-64/fontconfig-2.14.1-hc2a2eb6_0.tar.bz2#78415f0180a8d9c5bcc47889e00d5fb1
+https://conda.anaconda.org/conda-forge/linux-64/fontconfig-2.14.2-h14ed4e7_0.conda#0f69b688f52ff6da70bccb7ff7001d1d
 https://conda.anaconda.org/conda-forge/linux-64/freetype-2.12.1-hca18f0e_1.conda#e1232042de76d24539a436d37597eb06
 https://conda.anaconda.org/conda-forge/linux-64/freexl-1.0.6-h166bdaf_1.tar.bz2#897e772a157faf3330d72dd291486f62
-https://conda.anaconda.org/conda-forge/linux-64/gdal-3.6.2-py310hc1b7723_5.conda#2caf4a7da07394c02af296a19afed558
+https://conda.anaconda.org/conda-forge/linux-64/gdal-3.6.2-py310hc1b7723_7.conda#c0703fb11ff48f0e8408db47bcfdbd8a
 https://conda.anaconda.org/conda-forge/linux-64/geos-3.11.1-h27087fc_0.tar.bz2#917b9a50001fffdd89b321b5dba31e55
 https://conda.anaconda.org/conda-forge/linux-64/geotiff-1.7.1-h7a142b4_6.conda#b7963c107ed1f6a95cadc244f95cd3cd
 https://conda.anaconda.org/conda-forge/linux-64/gettext-0.21.1-h27087fc_0.tar.bz2#14947d8770185e5153fdd04d4673ed37
 https://conda.anaconda.org/conda-forge/linux-64/giflib-5.2.1-h36c2ea0_2.tar.bz2#626e68ae9cc5912d6adb79d318cf962d
-https://conda.anaconda.org/conda-forge/linux-64/git-2.39.1-pl5321h693f4a3_0.conda#12f9ef434e479d7ec9dcd3ab9799a49d
+https://conda.anaconda.org/conda-forge/linux-64/git-2.39.2-pl5321h693f4a3_0.conda#7311557191b94fb57e6f58a8d996b694
 https://conda.anaconda.org/conda-forge/linux-64/h5py-3.8.0-nompi_py310h0311031_100.conda#b1dfa8111fff8d0cc63091ce7b6dac38
 https://conda.anaconda.org/conda-forge/linux-64/hdf4-4.2.15-h9772cbc_5.tar.bz2#ee08782aff2ff9b3291c967fa6bc7336
 https://conda.anaconda.org/conda-forge/linux-64/hdf5-1.12.2-nompi_h4df4325_101.conda#162a25904af6586b234b2dd52ee99c61
@@ -33,39 +33,40 @@ https://conda.anaconda.org/conda-forge/linux-64/kealib-1.5.0-ha7026e8_0.conda#c9
 https://conda.anaconda.org/conda-forge/linux-64/keyutils-1.6.1-h166bdaf_0.tar.bz2#30186d27e2c9fa62b45fb1476b7200e3
 https://conda.anaconda.org/conda-forge/linux-64/krb5-1.20.1-h81ceb04_0.conda#89a41adce7106749573d883b2f657d78
 https://conda.anaconda.org/conda-forge/linux-64/lcms2-2.14-hfd0df8a_1.conda#c2566c2ea5f153ddd6bf4acaf7547d97
-https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.39-hcc3a1bd_1.conda#737be0d34c22d24432049ab7a3214de4
+https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.40-h41732ed_0.conda#7aca3059a1729aa76c597603f10b0dd3
 https://conda.anaconda.org/conda-forge/linux-64/lerc-4.0.0-h27087fc_0.tar.bz2#76bbff344f0134279f225174e9064c8f
 https://conda.anaconda.org/conda-forge/linux-64/libaec-1.0.6-hcb278e6_1.conda#0f683578378cddb223e7fd24f785ab2a
 https://conda.anaconda.org/conda-forge/linux-64/libblas-3.9.0-16_linux64_openblas.tar.bz2#d9b7a8639171f6c6fa0a983edabcfe2b
 https://conda.anaconda.org/conda-forge/linux-64/libcblas-3.9.0-16_linux64_openblas.tar.bz2#20bae26d0a1db73f758fc3754cab4719
-https://conda.anaconda.org/conda-forge/linux-64/libcurl-7.87.0-hdc1c0ab_0.conda#bc302fa1cf8eda15c60f669b7524a320
+https://conda.anaconda.org/conda-forge/linux-64/libcurl-7.88.1-hdc1c0ab_0.conda#81eaeb3b35163c8e90e57532bc93754d
 https://conda.anaconda.org/conda-forge/linux-64/libdeflate-1.17-h0b41bf4_0.conda#5cc781fd91968b11a8a7fdbee0982676
 https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20191231-he28a2e2_2.tar.bz2#4d331e44109e3f0e19b4cb8f9b82f3e1
 https://conda.anaconda.org/conda-forge/linux-64/libev-4.33-h516909a_1.tar.bz2#6f8720dff19e17ce5d48cfe7f3d2f0a3
 https://conda.anaconda.org/conda-forge/linux-64/libffi-3.4.2-h7f98852_5.tar.bz2#d645c6d2ac96843a2bfaccd2d62b3ac3
 https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2#d7c89558ba9fa0495403155b64376d81
 https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-12.2.0-h65d4601_19.tar.bz2#e4c94f80aef025c17ab0828cd85ef535
-https://conda.anaconda.org/conda-forge/linux-64/libgdal-3.6.2-h8c90c07_5.conda#b09f2d4b3b87a66b70df5778eed056b6
+https://conda.anaconda.org/conda-forge/linux-64/libgdal-3.6.2-hea5766e_7.conda#99afb7804b3c0b0e0be1cf0beb1fb1f1
 https://conda.anaconda.org/conda-forge/linux-64/libgfortran5-12.2.0-h337968e_19.tar.bz2#164b4b1acaedc47ee7e658ae6b308ca3
 https://conda.anaconda.org/conda-forge/linux-64/libgfortran-ng-12.2.0-h69a702a_19.tar.bz2#cd7a806282c16e1f2d39a7e80d3a3e0d
 https://conda.anaconda.org/conda-forge/linux-64/libglib-2.74.1-h606061b_1.tar.bz2#ed5349aa96776e00b34eccecf4a948fe
 https://conda.anaconda.org/conda-forge/linux-64/libgomp-12.2.0-h65d4601_19.tar.bz2#cedcee7c064c01c403f962c9e8d3c373
 https://conda.anaconda.org/conda-forge/linux-64/libiconv-1.17-h166bdaf_0.tar.bz2#b62b52da46c39ee2bc3c162ac7f1804d
+https://conda.anaconda.org/conda-forge/linux-64/libjpeg-turbo-2.1.4-h166bdaf_0.tar.bz2#b4f717df2d377410b462328bf0e8fb7d
 https://conda.anaconda.org/conda-forge/linux-64/libkml-1.3.0-h37653c0_1015.tar.bz2#37d3747dd24d604f63d2610910576e63
 https://conda.anaconda.org/conda-forge/linux-64/liblapack-3.9.0-16_linux64_openblas.tar.bz2#955d993f41f9354bf753d29864ea20ad
 https://conda.anaconda.org/conda-forge/linux-64/libllvm11-11.1.0-he0ac6c6_5.tar.bz2#cae79c6fd61cc6823cbebdbb2c16c60e
-https://conda.anaconda.org/conda-forge/linux-64/libnetcdf-4.8.1-nompi_h261ec11_106.tar.bz2#9b25de670ce5753a33c18b1090d1d3bf
-https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.51.0-hff17c54_0.conda#dd682f0b6d65e75b2bc868fc8e93d87e
+https://conda.anaconda.org/conda-forge/linux-64/libnetcdf-4.9.1-nompi_h34a3ff0_101.conda#9ce7f5bde721be2f82da9f861c57f1b3
+https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.52.0-h61bc06f_0.conda#613955a50485812985c059e7b269f42e
 https://conda.anaconda.org/conda-forge/linux-64/libnsl-2.0.0-h7f98852_0.tar.bz2#39b1328babf85c7c3a61636d9cd50206
 https://conda.anaconda.org/conda-forge/linux-64/libopenblas-0.3.21-pthreads_h78a6416_3.tar.bz2#8c5963a49b6035c40646a763293fbb35
 https://conda.anaconda.org/conda-forge/linux-64/libpng-1.6.39-h753d276_0.conda#e1c890aebdebbfbf87e2c917187b4416
-https://conda.anaconda.org/conda-forge/linux-64/libpq-15.1-hb675445_3.conda#9873ab80ec8fab4a2c26c7580e0d7f58
+https://conda.anaconda.org/conda-forge/linux-64/libpq-15.2-hb675445_0.conda#4654b17eccaba55b8581d6b9c77f53cc
 https://conda.anaconda.org/conda-forge/linux-64/librttopo-1.1.0-ha49c73b_12.tar.bz2#d2047c6de84b07a1db9cbe1683939956
 https://conda.anaconda.org/conda-forge/linux-64/libspatialite-5.0.1-h221c8f1_23.conda#859297085081cfbc123dc60015864f6b
 https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.40.0-h753d276_0.tar.bz2#2e5f9a37d487e1019fd4d8113adb2f9f
 https://conda.anaconda.org/conda-forge/linux-64/libssh2-1.10.0-hf14f497_3.tar.bz2#d85acad4b47dff4e3def14a769a97906
 https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-12.2.0-h46fd767_19.tar.bz2#1030b1f38c129f2634eae026f704fe60
-https://conda.anaconda.org/conda-forge/linux-64/libtiff-4.5.0-h6adf6a1_2.conda#2e648a34072eb39d7c4fc2a9981c5f0c
+https://conda.anaconda.org/conda-forge/linux-64/libtiff-4.5.0-h61b78d3_3.conda#8156bb7a598a505c9cbe4b7c0cf66117
 https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.32.1-h7f98852_1000.tar.bz2#772d69f030955d9646d3d0eaf21d859d
 https://conda.anaconda.org/conda-forge/linux-64/libwebp-base-1.2.4-h166bdaf_0.tar.bz2#ac2ccf7323d21f2994e4d1f5da664f37
 https://conda.anaconda.org/conda-forge/linux-64/libxcb-1.13-h7f98852_1004.tar.bz2#b3653fdc58d03face9724f602218a904
@@ -73,30 +74,30 @@ https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.10.3-h7463322_0.tar.bz
 https://conda.anaconda.org/conda-forge/linux-64/libzip-1.9.2-hc929e4a_1.tar.bz2#5b122b50e738c4be5c3f2899f010d7cf
 https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.2.13-h166bdaf_4.tar.bz2#f3f9de449d32ca9b9c66a22863c96f41
 https://conda.anaconda.org/conda-forge/linux-64/llvmlite-0.39.1-py310h58363a5_1.tar.bz2#18c28e036ae9c5366c56ecd0c875a108
-https://conda.anaconda.org/conda-forge/linux-64/lz4-c-1.9.3-h9c3ff4c_1.tar.bz2#fbe97e8fa6f275d7c76a09e795adc3e6
+https://conda.anaconda.org/conda-forge/linux-64/lz4-c-1.9.4-hcb278e6_0.conda#318b08df404f9c9be5712aaa5a6f0bb0
 https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.3-h27087fc_1.tar.bz2#4acfc691e64342b9dae57cf2adc63238
 https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda#da0ec11a6454ae19bff5b02ed881a2b1
-https://conda.anaconda.org/conda-forge/linux-64/nss-3.82-he02c5a1_0.conda#f8d7f11d19e4cb2207eab159fd4c0152
+https://conda.anaconda.org/conda-forge/linux-64/nss-3.88-he45b914_0.conda#d7a81dfb99ad8fbb88872fb7ec646e6c
 https://conda.anaconda.org/conda-forge/linux-64/numba-0.56.4-py310ha5257ce_0.conda#5a4904c5c57a3d1f9a5ae6533401f652
 https://conda.anaconda.org/conda-forge/linux-64/numpy-1.23.5-py310h53a5b5f_0.conda#3b114b1559def8bad228fec544ac1812
 https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.0-hfec8fc6_2.conda#5ce6a42505c6e9e6151c54c3ec8d68ea
 https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-2_gnu.tar.bz2#73aaf86a425cc6e73fcf236a5a46396d
-https://conda.anaconda.org/conda-forge/linux-64/openssl-3.0.7-h0b41bf4_1.conda#7adaac6ff98219bcb99b45e408b80f4e
+https://conda.anaconda.org/conda-forge/linux-64/openssl-3.0.8-h0b41bf4_0.conda#e043403cd18faf815bf7705ab6c1e092
 https://conda.anaconda.org/conda-forge/linux-64/pcre2-10.40-hc3806b6_0.tar.bz2#69e2c796349cd9b273890bee0febfe1b
 https://conda.anaconda.org/conda-forge/linux-64/perl-5.32.1-2_h7f98852_perl5.tar.bz2#09ba115862623f00962e9809ea248f1a
 https://conda.anaconda.org/conda-forge/linux-64/pixman-0.40.0-h36c2ea0_0.tar.bz2#660e72c82f2e75a6b3fe6a6e75c79f19
 https://conda.anaconda.org/conda-forge/linux-64/poppler-23.01.0-h091648b_0.conda#33ba6d8025df115bcbe50c69e9b808ed
-https://conda.anaconda.org/conda-forge/linux-64/postgresql-15.1-h3248436_3.conda#4f686d5d582a3e3b32a3cb010fa12103
+https://conda.anaconda.org/conda-forge/linux-64/postgresql-15.2-h3248436_0.conda#4dbd6c4bc33369751a4d728b392943ba
 https://conda.anaconda.org/conda-forge/linux-64/proj-9.1.1-h8ffa02c_2.conda#c264aea0e16bba26afa0a0940e954492
 https://conda.anaconda.org/conda-forge/linux-64/pthread-stubs-0.4-h36c2ea0_1001.tar.bz2#22dad4df6e8630e8dff2428f6f6a7036
-https://conda.anaconda.org/conda-forge/linux-64/pydantic-1.10.4-py310h1fa729e_1.conda#9141df307497cf8de501bfd0205c6d1b
+https://conda.anaconda.org/conda-forge/linux-64/pydantic-1.10.5-py310h1fa729e_0.conda#1af0b4a028ef4812010304e94ff04f92
 https://conda.anaconda.org/conda-forge/linux-64/pyproj-3.4.1-py310h15e2413_1.conda#5be35366687def87437d210fd673100c
-https://conda.anaconda.org/conda-forge/linux-64/python-3.10.8-h4a9ceb5_0_cpython.conda#be2a6d78752c2ab85f360ce37d2c64e2
+https://conda.anaconda.org/conda-forge/linux-64/python-3.10.9-he550d4f_0_cpython.conda#3cb3e91b3fe66baa68a12c85f39b9b40
 https://conda.anaconda.org/conda-forge/linux-64/python_abi-3.10-3_cp310.conda#4eb33d14d794b0f4be116443ffed3853
 https://conda.anaconda.org/conda-forge/linux-64/readline-8.1.2-h0f457ee_0.tar.bz2#db2ebbe2943aae81ed051a6a9af8e0fa
 https://conda.anaconda.org/conda-forge/linux-64/ruamel.yaml-0.17.21-py310h5764c6d_2.tar.bz2#93643151fba5c31c676effdd7661f075
 https://conda.anaconda.org/conda-forge/linux-64/ruamel.yaml.clib-0.2.7-py310h1fa729e_1.conda#2f9b517412af46255cef5e53a22c264e
-https://conda.anaconda.org/conda-forge/linux-64/scipy-1.10.0-py310h8deb116_0.conda#ef72eeddf5316330730b11907c6c07d8
+https://conda.anaconda.org/conda-forge/linux-64/scipy-1.10.1-py310h8deb116_0.conda#4c9604c5ec179c21f8f0a09e3c164480
 https://conda.anaconda.org/conda-forge/linux-64/snappy-1.1.9-hbd366e4_2.tar.bz2#48018e187dacc6002d3ede9c824238ac
 https://conda.anaconda.org/conda-forge/linux-64/sqlite-3.40.0-h4ff8645_0.tar.bz2#bb11803129cbbb53ed56f9506ff74145
 https://conda.anaconda.org/conda-forge/linux-64/tiledb-2.13.2-hd532e3d_0.conda#6d97164f19dbd27575ef1899b02dc1e0
@@ -106,23 +107,21 @@ https://conda.anaconda.org/conda-forge/linux-64/xerces-c-3.2.4-h55805fa_1.tar.bz
 https://conda.anaconda.org/conda-forge/linux-64/xorg-kbproto-1.0.7-h7f98852_1002.tar.bz2#4b230e8381279d76131116660f5a241a
 https://conda.anaconda.org/conda-forge/linux-64/xorg-libice-1.0.10-h7f98852_0.tar.bz2#d6b0b50b49eccfe0be0373be628be0f3
 https://conda.anaconda.org/conda-forge/linux-64/xorg-libsm-1.2.3-hd9c2040_1000.tar.bz2#9e856f78d5c80d5a78f61e72d1d473a3
-https://conda.anaconda.org/conda-forge/linux-64/xorg-libx11-1.7.2-h7f98852_0.tar.bz2#12a61e640b8894504326aadafccbb790
+https://conda.anaconda.org/conda-forge/linux-64/xorg-libx11-1.8.4-h0b41bf4_0.conda#ea8fbfeb976ac49cbeb594e985393514
 https://conda.anaconda.org/conda-forge/linux-64/xorg-libxau-1.0.9-h7f98852_0.tar.bz2#bf6f803a544f26ebbdc3bfff272eb179
 https://conda.anaconda.org/conda-forge/linux-64/xorg-libxdmcp-1.1.3-h7f98852_0.tar.bz2#be93aabceefa2fac576e971aef407908
-https://conda.anaconda.org/conda-forge/linux-64/xorg-libxext-1.3.4-h7f98852_1.tar.bz2#536cc5db4d0a3ba0630541aec064b5e4
+https://conda.anaconda.org/conda-forge/linux-64/xorg-libxext-1.3.4-h0b41bf4_2.conda#82b6df12252e6f32402b96dacc656fec
 https://conda.anaconda.org/conda-forge/linux-64/xorg-libxrender-0.9.10-h7f98852_1003.tar.bz2#f59c1242cc1dd93e72c2ee2b360979eb
 https://conda.anaconda.org/conda-forge/linux-64/xorg-renderproto-0.11.1-h7f98852_1002.tar.bz2#06feff3d2634e3097ce2fe681474b534
-https://conda.anaconda.org/conda-forge/linux-64/xorg-xextproto-7.3.0-h7f98852_1002.tar.bz2#1e15f6ad85a7d743a2ac68dae6c82b98
+https://conda.anaconda.org/conda-forge/linux-64/xorg-xextproto-7.3.0-h0b41bf4_1003.conda#bce9f945da8ad2ae9b1d7165a64d0f87
 https://conda.anaconda.org/conda-forge/linux-64/xorg-xproto-7.0.31-h7f98852_1007.tar.bz2#b4a4381d54784606820704f7b5f05a15
 https://conda.anaconda.org/conda-forge/linux-64/xz-5.2.6-h166bdaf_0.tar.bz2#2161070d867d1b1204ea749c8eec4ef0
 https://conda.anaconda.org/conda-forge/linux-64/zlib-1.2.13-h166bdaf_4.tar.bz2#4b11e365c0275b808be78b30f904e295
 https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.2-h3eb15da_6.conda#6b63daed8feeca47be78f323e793d555
-https://conda.anaconda.org/conda-forge/noarch/appdirs-1.4.4-pyh9f0ad1d_0.tar.bz2#5f095bc6454094e96f146491fd03633b
 https://conda.anaconda.org/conda-forge/noarch/cached-property-1.5.2-hd8ed1ab_1.tar.bz2#9b347a7ec10940d3f7941ff6c460b551
 https://conda.anaconda.org/conda-forge/noarch/cached_property-1.5.2-pyha770c72_1.tar.bz2#576d629e47797577ab0f1b351297ef4a
 https://conda.anaconda.org/conda-forge/noarch/certifi-2022.12.7-pyhd8ed1ab_0.conda#fb9addc3db06e56abe03e0e9f21a63e6
 https://conda.anaconda.org/conda-forge/noarch/charset-normalizer-2.1.1-pyhd8ed1ab_0.tar.bz2#c1d5b294fbf9a795dec349a6f4d8be8e
-https://conda.anaconda.org/conda-forge/noarch/colorama-0.4.6-pyhd8ed1ab_0.tar.bz2#3faab06a954c2a04039983f2c4a50d99
 https://conda.anaconda.org/conda-forge/noarch/fonts-conda-ecosystem-1-0.tar.bz2#fee5683a3f04bd15cbd8318b096a27ab
 https://conda.anaconda.org/conda-forge/noarch/fonts-conda-forge-1-0.tar.bz2#f766549260d6815b0c52253f1fb1bb29
 https://conda.anaconda.org/conda-forge/noarch/font-ttf-dejavu-sans-mono-2.37-hab24e00_0.tar.bz2#0c96522c6bdaed4b1566d11387caaf45
@@ -130,17 +129,22 @@ https://conda.anaconda.org/conda-forge/noarch/font-ttf-inconsolata-3.000-h77eed3
 https://conda.anaconda.org/conda-forge/noarch/font-ttf-source-code-pro-2.038-h77eed37_0.tar.bz2#4d59c254e01d9cde7957100457e2d5fb
 https://conda.anaconda.org/conda-forge/noarch/font-ttf-ubuntu-0.83-hab24e00_0.tar.bz2#19410c3df09dfb12d1206132a1d357c5
 https://conda.anaconda.org/conda-forge/noarch/idna-3.4-pyhd8ed1ab_0.tar.bz2#34272b248891bddccc64479f9a7fffed
+https://conda.anaconda.org/conda-forge/noarch/markdown-it-py-2.2.0-pyhd8ed1ab_0.conda#b2928a6c6d52d7e3562b4a59c3214e3a
+https://conda.anaconda.org/conda-forge/noarch/mdurl-0.1.0-pyhd8ed1ab_0.tar.bz2#f8dab71fdc13b1bf29a01248b156d268
 https://conda.anaconda.org/conda-forge/noarch/packaging-23.0-pyhd8ed1ab_0.conda#1ff2e3ca41f0ce16afec7190db28288b
-https://conda.anaconda.org/conda-forge/noarch/pip-22.3.1-pyhd8ed1ab_0.tar.bz2#da66f2851b9836d3a7c5190082a45f7d
-https://conda.anaconda.org/conda-forge/noarch/pooch-1.6.0-pyhd8ed1ab_0.tar.bz2#6429e1d1091c51f626b5dcfdd38bf429
-https://conda.anaconda.org/conda-forge/noarch/poppler-data-0.4.11-hd8ed1ab_0.tar.bz2#abc27381c4f005da588cffa1f76403ee
+https://conda.anaconda.org/conda-forge/noarch/pip-23.0.1-pyhd8ed1ab_0.conda#8025ca83b8ba5430b640b83917c2a6f7
+https://conda.anaconda.org/conda-forge/noarch/platformdirs-3.1.0-pyhd8ed1ab_0.conda#0b8fbdfd52918bc2f1b76feccd95c919
+https://conda.anaconda.org/conda-forge/noarch/pooch-1.7.0-pyhd8ed1ab_0.conda#eb2e0fc33ad0d04b51f9280360c13c1e
+https://conda.anaconda.org/conda-forge/noarch/poppler-data-0.4.12-hd8ed1ab_0.conda#d8d7293c5b37f39b2ac32940621c6592
 https://conda.anaconda.org/conda-forge/noarch/pycparser-2.21-pyhd8ed1ab_0.tar.bz2#076becd9e05608f8dc72757d5f3a91ff
+https://conda.anaconda.org/conda-forge/noarch/pygments-2.14.0-pyhd8ed1ab_0.conda#c78cd16b11cd6a295484bd6c8f24bea1
 https://conda.anaconda.org/conda-forge/noarch/pymp-pypi-0.4.5-pyhd8ed1ab_0.tar.bz2#3d56b4c5a162223b9f20cef1fdc8a89a
 https://conda.anaconda.org/conda-forge/noarch/pyopenssl-23.0.0-pyhd8ed1ab_0.conda#d41957700e83bbb925928764cb7f8878
 https://conda.anaconda.org/conda-forge/noarch/pysocks-1.7.1-pyha2e5f31_6.tar.bz2#2a7de29fb590ca14b5243c4c812c8025
 https://conda.anaconda.org/conda-forge/noarch/requests-2.28.2-pyhd8ed1ab_0.conda#11d178fc55199482ee48d6812ea83983
-https://conda.anaconda.org/conda-forge/noarch/setuptools-66.1.1-pyhd8ed1ab_0.conda#9467d520d1457018e055bbbfdf9b7567
-https://conda.anaconda.org/conda-forge/noarch/tqdm-4.64.1-pyhd8ed1ab_0.tar.bz2#5526ff3f88f9db87bb0924b9ce575345
+https://conda.anaconda.org/conda-forge/noarch/rich-13.3.2-pyhd8ed1ab_0.conda#0aec87e5ad4d9170d58d3f0525ff5825
+https://conda.anaconda.org/conda-forge/noarch/setuptools-67.5.1-pyhd8ed1ab_0.conda#5aec57dd42104f4916c397fe291f9697
+https://conda.anaconda.org/conda-forge/noarch/threadpoolctl-3.1.0-pyh8a188c0_0.tar.bz2#a2995ee828f65687ac5b1e71a2ab1e0c
 https://conda.anaconda.org/conda-forge/noarch/typing-extensions-4.4.0-hd8ed1ab_0.tar.bz2#be969210b61b897775a0de63cd9e9026
 https://conda.anaconda.org/conda-forge/noarch/typing_extensions-4.4.0-pyha770c72_0.tar.bz2#2d93b130d148d7fc77e583677792fc6a
 https://conda.anaconda.org/conda-forge/noarch/tzdata-2022g-h191b570_0.conda#51fc4fcfb19f5d95ffc8c339db5068e8
diff --git a/pyproject.toml b/pyproject.toml
index 689e4812..9ae8167b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -66,6 +66,7 @@ ignore_missing_imports = true
 ignore = "D100,D102,D104,D105,D106,D107,D203,D204,D213,D413"
 
 [tool.pytest.ini_options]
+doctest_optionflags = "NORMALIZE_WHITESPACE"
 addopts = "  --cov=dolphin  -n auto --maxprocesses=8 --doctest-modules --ignore=scripts --ignore=docs --ignore=data"
 filterwarnings = [
   "error",
diff --git a/requirements.txt b/requirements.txt
index 46f8517a..15f03b23 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,6 +5,7 @@ numpy>=1.20
 pydantic>=1.10
 pymp-pypi>=0.4.5
 pyproj>=3.2
+rich>=12.0
 ruamel_yaml>=0.15
 scipy>=1.5
-tqdm>=4.60
+threadpoolctl>=3.0
diff --git a/scripts/benchmarking-cpu.sh b/scripts/benchmarking-cpu.sh
new file mode 100755
index 00000000..3287e329
--- /dev/null
+++ b/scripts/benchmarking-cpu.sh
@@ -0,0 +1,37 @@
+#/usr/bin/env bash
+
+set -e
+
+SLCDIR="/u/aurora-r0/staniewi/dev/hawaii/frame97/stack-no-lut-copy"
+NSLC="15 27"
+for blocksize_gb in 1 4; do
+    for s in 2 3; do
+        for nworkers in 4 16; do
+            for tpw in 1 8 32; do
+                for n in $NSLC; do
+                    slcfiles=$(find $SLCDIR -name "*185684_iw2*.h5" | head -n $n | paste -sd " ")
+                    strides="$((2 * s)) $s"
+                    outfile="dolphin_config_cpu_block${blocksize_gb}GB_strides${s}_tpw${tpw}_nslc${n}_nworkers${nworkers}.yaml"
+
+                    dolphin config --no-gpu --amp-dispersion-threshold 0.25 --single-update \
+                        --slc-files $slcfiles \
+                        --block-size-gb $blocksize_gb \
+                        --strides $strides \
+                        --threads-per-worker $tpw \
+                        --n-workers $nworkers \
+                        --outfile $outfile
+                done
+            done
+        done
+    done
+done
+
+for f in *.yaml; do
+    echo "Running $f"
+    echo "Removing the scratch/linked_phase directory"
+    rm -rf scratch/linked_phase scratch/slc_stack.vrt
+    echo "Removing the SLCs from the cache"
+    vmtouch -e $SLCDIR
+
+    dolphin run $f &>${f/.yaml/.log}
+done
diff --git a/scripts/benchmarking-gpu.sh b/scripts/benchmarking-gpu.sh
new file mode 100755
index 00000000..558b27ea
--- /dev/null
+++ b/scripts/benchmarking-gpu.sh
@@ -0,0 +1,34 @@
+#/usr/bin/env bash
+
+set -e
+
+SLCDIR="/u/aurora-r0/staniewi/dev/hawaii/frame97/stack-no-lut"
+NSLC="15 27"
+for blocksize_gb in 1 4; do
+    for s in 2 3; do # strides
+        for tpw in 4 32; do
+            for n in $NSLC; do
+                slcfiles=$(find $SLCDIR -name "*185684_iw2*.h5" | head -n $n | paste -sd " ")
+                strides="$((2 * s)) $s"
+
+                outfile="dolphin_config_gpu_block${blocksize_gb}GB_strides${s}_tpw${tpw}_nslc${n}.yaml"
+                dolphin config --amp-dispersion-threshold 0.25 --single-update \
+                    --slc-files $slcfiles \
+                    --block-size-gb $blocksize_gb \
+                    --strides $strides \
+                    --threads-per-worker $tpw \
+                    --outfile $outfile
+            done
+        done
+    done
+done
+
+for f in *.yaml; do
+    echo "Running $f"
+    echo "Removing the scratch/linked_phase directory"
+    rm -rf scratch/linked_phase scratch/slc_stack.vrt
+    echo "Removing the SLCs from the cache"
+    vmtouch -e $SLCDIR
+
+    dolphin run $f &>${f/.yaml/.log}
+done
diff --git a/scripts/parse_benchmark_logs.py b/scripts/parse_benchmark_logs.py
new file mode 100644
index 00000000..f364aa3c
--- /dev/null
+++ b/scripts/parse_benchmark_logs.py
@@ -0,0 +1,59 @@
+import re
+from pathlib import Path
+
+import pandas as pd
+
+
+def get_df(path):
+    """Create a dataframe from a directory of log files."""
+    rows = []
+    for f in Path(path).glob("*.log"):
+        result = _get_memory(f)
+        attr = _parse_filename(f)
+        result.update(attr)
+        rows.append(result)
+    return pd.DataFrame(rows)
+
+
+def _get_memory(logfile):
+    out = {"file": str(Path(logfile).resolve())}
+    mempat = r"Maximum memory usage: (\d\.\d{2}) GB"
+    timepat = (
+        r"Total elapsed time for dolphin.workflows.s1_disp.run : (\d*\.\d{2}) minutes"
+        r" \((\d*\.\d{2}) seconds\)"
+    )
+    for line in open(logfile).readlines():
+        m = re.search(mempat, line)
+        if m:
+            out["memory"] = float(m.groups()[0])
+            continue
+        m = re.search(timepat, line)
+        if m:
+            out["runtime"] = float(m.groups()[1])
+    return out
+
+
+def _parse_filename(logfile):
+    # dolphin_config_cpu_block1GB_strides2_tpw8_nslc27_nworkers4.log
+    # dolphin_config_gpu_block1GB_strides2_tpw32_nslc15.log
+    # Create the regular expression pattern with named groups
+    pat = (
+        r"dolphin_config_"
+        r"(?P<device>\w+)_"
+        r"block(?P<block>\d+)GB_"
+        r"strides(?P<strides>\d+)_"
+        r"tpw(?P<tpw>\d+)_"
+        r"nslc(?P<nslc>\d+)"
+        # nworkers is optional
+        r"(_nworkers(?P<nworkers>\d+))?"
+    )
+    m = re.match(pat, Path(logfile).stem)
+    if not m:
+        raise ValueError(f"Could not parse filename: {logfile}")
+    out_dict = m.groupdict()
+    # convert numbers to ints
+    for k, v in out_dict.items():
+        if v is not None and k != "device":
+            out_dict[k] = int(v)
+
+    return out_dict
diff --git a/src/dolphin/_background.py b/src/dolphin/_background.py
index 131252a3..063b7ff1 100644
--- a/src/dolphin/_background.py
+++ b/src/dolphin/_background.py
@@ -1,4 +1,6 @@
 import abc
+import os
+import time
 from queue import Empty, Full, Queue
 from threading import Event, Thread
 from threading import enumerate as threading_enumerate
@@ -59,32 +61,56 @@ def __init__(
         timeout=_DEFAULT_TIMEOUT,
         name="BackgroundWorker",
     ):
+        self.name = name
         self.store_results = store_results
         self.timeout = timeout
         self._finished_event = Event()
         self._work_queue = Queue(num_work_queue)
-        self._results_queue = Queue(num_results_queue)
+        if self.store_results:
+            self._results_queue = Queue(num_results_queue)
         self._thread = Thread(target=self._consume_work_queue, name=name)
         self._thread.start()
         self._drop_unfinished_results = drop_unfinished_results
 
     def _consume_work_queue(self):
-        while not self._finished_event.is_set() and is_main_thread_active():
-            logger.debug("getting work")
+        while True:
+            if not is_main_thread_active():
+                break
+
+            logger.debug(f"{self.name} getting work")
+            if self._finished_event.is_set():
+                do_exit = self._drop_unfinished_results or (
+                    self._work_queue.unfinished_tasks == 0
+                )
+                if do_exit:
+                    break
+                else:
+                    # Keep going even if finished event is set
+                    logger.debug(
+                        f"{self.name} Finished... but waiting for work queue to empty,"
+                        f" {self._work_queue.qsize()} items left,"
+                        f" {self._work_queue.unfinished_tasks} unfinished"
+                    )
             try:
                 args, kw = self._work_queue.get(timeout=self.timeout)
+                logger.debug(f"{self.name} processing")
+                result = self.process(*args, **kw)
+                self._work_queue.task_done()
+                # Notify the queue that processing is done
+                logger.debug(f"{self.name} got result")
             except Empty:
-                logger.debug("timed out, checking if done")
+                logger.debug(f"{self.name} timed out, checking if done")
                 continue
-            logger.debug("processing")
-            result = self.process(*args, **kw)
-            logger.debug("got result")
+
             if self.store_results:
-                logger.debug("saving result in queue")
-                try:
-                    self._results_queue.put(result, timeout=2)
-                except Full:
-                    logger.debug("result queue full, waiting...")
+                logger.debug(f"{self.name} saving result in queue")
+                while True:
+                    try:
+                        self._results_queue.put(result, timeout=2)
+                        break
+                    except Full:
+                        logger.debug(f"{self.name} result queue full, waiting...")
+                        continue
 
     @abc.abstractmethod
     def process(self, *args, **kw):
@@ -107,8 +133,16 @@ def get_result(self):
         Blocks until a result is available.
         Same output interface as `process`.
         """
-        result = self._results_queue.get()
-        self._results_queue.task_done()
+        while True:
+            try:
+                result = self._results_queue.get(timeout=self.timeout)
+                self._results_queue.task_done()
+                break
+            except Empty:
+                logger.debug(f"{self.name} get_result timed out, checking if done")
+                if self._finished_event.is_set():
+                    raise RuntimeError("Attempted to get_result after notify_finished!")
+                continue
         return result
 
     def notify_finished(self, timeout=None):
@@ -119,7 +153,7 @@ def notify_finished(self, timeout=None):
         If `store_results=True` also block until all results have been retrieved.
         """
         self._finished_event.set()
-        if not self._drop_unfinished_results:
+        if self.store_results and not self._drop_unfinished_results:
             self._results_queue.join()
         self._thread.join(timeout)
 
@@ -227,3 +261,74 @@ def process(self, *args, **kw):
     def read(self, *args, **kw):
         """User-defined method for reading a chunk of data."""
         pass
+
+
+class NvidiaMemoryWatcher(Thread):
+    """Watch the memory usage of the GPU and log it to a file.
+
+    Parameters
+    ----------
+    log_file : str
+        The file to write the memory usage to.
+    refresh_rate : float, optional
+        The refresh_rate in seconds to check the memory usage, by default 1.0
+    """
+
+    def __init__(
+        self,
+        log_file: str = "nvidia_memory.log",
+        refresh_rate: float = 1.0,
+        gpu_id: int = 0,
+    ):
+        try:
+            from pynvml.smi import nvidia_smi  # noqa: F401
+        except ImportError:
+            raise ImportError("Please install pynvml through pip or conda")
+
+        super().__init__(name="NvidiaMemoryWatcher")
+        self.log_file = log_file
+        self.pid = os.getpid()
+        self.t0 = time.time()
+        self.refresh_rate = refresh_rate
+        self.gpu_id = gpu_id
+        # The query lag is the time it takes to query the GPU memory
+        # This is used to try and refresh close to the refresh rate
+        self._query_lag = 0.5
+        self._finished_event = Event()
+        self._thread = Thread(target=self.run)
+
+        self._thread.start()
+
+    def run(self):
+        """Run the background task."""
+        logger.info(
+            f"Logging GPU memory usage to {self.log_file} every {self.refresh_rate} s"
+        )
+        with open(self.log_file, "w") as f:
+            # Write the header
+            f.write("time(s),memory(GB)\n")
+
+        while not self._finished_event.is_set() and is_main_thread_active():
+            mem = self._get_gpu_memory()
+            t_cur = time.time() - self.t0
+            with open(self.log_file, "a") as f:
+                row = f"{t_cur:.3f},{mem:.2f}\n"
+                f.write(row)
+
+            # Sleep until the next refresh
+            time.sleep(max(0, self.refresh_rate - self._query_lag))
+
+    def join(self):
+        """Wait for the thread to finish."""
+        self._thread.join()
+
+    def notify_finished(self):
+        """Notify the thread that it should finish."""
+        self._finished_event.set()
+        self._thread.join()
+
+    def _get_gpu_memory(self) -> float:
+        """Get the memory usage (in GiB) of the GPU for the current pid."""
+        from dolphin.utils import get_gpu_memory
+
+        return get_gpu_memory(pid=self.pid, gpu_id=self.gpu_id)
diff --git a/src/dolphin/_log.py b/src/dolphin/_log.py
index 509bb198..1b298cc5 100644
--- a/src/dolphin/_log.py
+++ b/src/dolphin/_log.py
@@ -23,6 +23,8 @@
 from logging import Formatter
 from typing import Optional
 
+from rich.logging import RichHandler
+
 __all__ = ["get_log", "log_runtime"]
 
 
@@ -48,47 +50,36 @@ def get_log(
     logging.Logger
     """
     logger = logging.getLogger(name)
-    return format_log(logger, debug=debug, filename=filename)
+    if not logger.hasHandlers():
+        setup_logging(debug=debug)
+    if debug:
+        logger.setLevel(logging.DEBUG)
+
+    # In addition to stderr, log to a file if requested
+    if filename:
+        file_handler = logging.FileHandler(filename)
+        file_handler.setLevel(logging.DEBUG)
+        formatter = Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s")
+        file_handler.setFormatter(formatter)
+        logger.addHandler(file_handler)
+    return logger
 
 
-def format_log(
-    logger: logging.Logger, debug: bool = False, filename: Optional[str] = None
-) -> logging.Logger:
+def setup_logging(debug: bool = False) -> None:
     """Make the logging output pretty and colored with times.
 
     Parameters
     ----------
-    logger : logging.Logger
-        The logger to format
     debug : bool (Default value = False)
         If true, sets logging level to DEBUG
-    filename : str, optional
-        If provided, will log to this file in addition to stderr.
 
-    Returns
-    -------
-    logging.Logger
     """
+    # Set for all dolphin modules
+    logger = logging.getLogger("dolphin")
+    h = RichHandler(rich_tracebacks=True, log_time_format="[%Y-%m-%d %H:%M:%S]")
+    logger.addHandler(h)
     log_level = logging.DEBUG if debug else logging.INFO
-    format_ = "[%(asctime)s] [%(levelname)s %(filename)s] %(message)s"
-    formatter = Formatter(format_, datefmt="%m/%d %H:%M:%S")
-
-    handler = logging.StreamHandler()
-    handler.setFormatter(formatter)
-
-    if not logger.handlers:
-        logger.addHandler(handler)
-        logger.setLevel(log_level)
-
-    if filename is not None:
-        file_handler = logging.FileHandler(filename)
-        file_handler.setFormatter(formatter)
-        logger.addHandler(file_handler)
-
-    if debug:
-        logger.setLevel(debug)
-
-    return logger
+    logger.setLevel(log_level)
 
 
 def log_runtime(f: Callable) -> Callable:
@@ -109,14 +100,15 @@ def wrapper(*args, **kwargs):
         result = f(*args, **kwargs)
 
         t2 = time.time()
-        elapsed_time = t2 - t1
-        time_string = "Total elapsed time for {} : {} minutes ({} seconds)".format(
-            f.__name__,
-            "{0:.2f}".format(elapsed_time / 60.0),
-            "{0:.2f}".format(elapsed_time),
+        elapsed_seconds = t2 - t1
+        elapsed_minutes = elapsed_seconds / 60.0
+        time_string = (
+            f"Total elapsed time for {f.__module__}.{f.__name__} : "
+            f"{elapsed_minutes:.2f} minutes ({elapsed_seconds:.2f} seconds)"
         )
 
         logger.info(time_string)
+
         return result
 
     return wrapper
diff --git a/src/dolphin/io.py b/src/dolphin/io.py
index 5ac4551d..25c0e8d1 100644
--- a/src/dolphin/io.py
+++ b/src/dolphin/io.py
@@ -13,12 +13,11 @@
 from numpy.typing import ArrayLike, DTypeLike
 from osgeo import gdal
 from pyproj import CRS
-from tqdm.auto import tqdm
 
 from dolphin._background import _DEFAULT_TIMEOUT, BackgroundReader, BackgroundWriter
 from dolphin._log import get_log
 from dolphin._types import Filename
-from dolphin.utils import gdal_to_numpy_type, numpy_to_gdal_type
+from dolphin.utils import gdal_to_numpy_type, numpy_to_gdal_type, progress
 
 gdal.UseExceptions()
 
@@ -598,6 +597,7 @@ def __init__(
         self._queue_size = queue_size
         self._skip_empty = skip_empty
         self._nodata_mask = nodata_mask
+        self._block_shape = block_shape
         self._nodata = get_raster_nodata(filename)
         if self._nodata is None:
             self._nodata = np.nan
@@ -614,23 +614,26 @@ def iter_blocks(
         for rows, cols in self.slices:
             self.queue_read(rows, cols)
 
-        for _ in tqdm(range(len(self.slices))):
-            cur_block, (rows, cols) = self.get_data()
-            logger.debug(f"got data for {rows, cols}: {cur_block.shape}")
-
-            if self._skip_empty and self._nodata_mask is not None:
-                if self._nodata_mask[rows, cols].all():
-                    continue
-
-            if self._skip_empty:
-                # Otherwise look at the actual block we loaded
-                if np.isnan(self._nodata):
-                    block_nodata = np.isnan(cur_block)
-                else:
-                    block_nodata = cur_block == self._nodata
-                if np.all(block_nodata):
-                    continue
-            yield cur_block, (rows, cols)
+        s_iter = range(len(self.slices))
+        desc = f"Processing {self._block_shape} sized blocks..."
+        with progress() as p:
+            for _ in p.track(s_iter, description=desc):
+                cur_block, (rows, cols) = self.get_data()
+                logger.debug(f"got data for {rows, cols}: {cur_block.shape}")
+
+                if self._skip_empty and self._nodata_mask is not None:
+                    if self._nodata_mask[rows, cols].all():
+                        continue
+
+                if self._skip_empty:
+                    # Otherwise look at the actual block we loaded
+                    if np.isnan(self._nodata):
+                        block_nodata = np.isnan(cur_block)
+                    else:
+                        block_nodata = cur_block == self._nodata
+                    if np.all(block_nodata):
+                        continue
+                yield cur_block, (rows, cols)
 
         self.notify_finished()
 
diff --git a/src/dolphin/phase_link/_mle_gpu.py b/src/dolphin/phase_link/_mle_gpu.py
index a2daa57d..27ebdc94 100644
--- a/src/dolphin/phase_link/_mle_gpu.py
+++ b/src/dolphin/phase_link/_mle_gpu.py
@@ -98,6 +98,9 @@ def run_gpu(
     if free_mem:
         del d_slc_stack
         del d_C_arrays
+        del d_output_phase
+        del d_cpx_phase
+        cp.get_default_memory_pool().free_all_blocks()
         cp.get_default_pinned_memory_pool().free_all_blocks()
 
     if use_slc_amp:
diff --git a/src/dolphin/phase_link/mle.py b/src/dolphin/phase_link/mle.py
index c038061a..43a9ab8f 100644
--- a/src/dolphin/phase_link/mle.py
+++ b/src/dolphin/phase_link/mle.py
@@ -7,12 +7,7 @@
 from scipy.linalg import eigh
 
 from dolphin._types import Filename
-from dolphin.utils import (
-    get_array_module,
-    gpu_is_available,
-    take_looks,
-    upsample_nearest,
-)
+from dolphin.utils import get_array_module, gpu_is_available, take_looks
 
 logger = logging.getLogger(__name__)
 
@@ -284,7 +279,9 @@ def _check_all_nans(slc_stack):
     # Check that there are no SLCS which are all nans:
     bad_slc_idxs = np.where(np.all(nans, axis=(1, 2)))[0]
     if bad_slc_idxs.size > 0:
-        raise PhaseLinkRuntimeError(f"SLC stack[{bad_slc_idxs}] has are all NaNs.")
+        raise PhaseLinkRuntimeError(
+            f"slc_stack[{bad_slc_idxs}] out of {len(slc_stack)} are all NaNs."
+        )
 
 
 def _fill_ps_pixels(
@@ -332,20 +329,19 @@ def _fill_ps_pixels(
     ps_mask_looked = ps_mask_looked[: mle_est.shape[1], : mle_est.shape[2]]
 
     if use_max_ps:
+        print("Using max PS pixel to fill in MLE estimate")
         # Get the indices of the brightest pixels within each look window
-        slc_r_idxs, slc_c_idxs = _get_maxes(mag, strides["y"], strides["x"])
+        slc_r_idxs, slc_c_idxs = _get_max_idxs(mag, strides["y"], strides["x"])
         # we're only filling where there are PS pixels
         ref = np.conj(slc_stack[0][slc_r_idxs, slc_c_idxs])
         for i in range(len(slc_stack)):
             mle_est[i][ps_mask_looked] = slc_stack[i][slc_r_idxs, slc_c_idxs] * ref
     else:
         # Get the average of all PS pixels within each look window
-
+        # The referencing to SLC 0 is done in _get_avg_ps
         avg_ps = _get_avg_ps(slc_stack, ps_mask, strides)[
             :, : mle_est.shape[1], : mle_est.shape[2]
         ]
-        ref = np.conj(avg_ps[0])[None, :, :]
-        avg_ps *= ref
         mle_est[:, ps_mask_looked] = avg_ps[:, ps_mask_looked]
 
     # Force PS pixels to have high temporal coherence
@@ -358,6 +354,8 @@ def _get_avg_ps(
     # First, set all non-PS pixels to NaN
     slc_stack_nanned = slc_stack.copy()
     slc_stack_nanned[:, ~ps_mask] = np.nan
+    # Reference all ps pixels in the SLC stack to the first SLC
+    slc_stack_nanned[:, ps_mask] *= np.conj(slc_stack_nanned[0, ps_mask])[None]
     # Then, take the average of all PS pixels within each look window
     return take_looks(
         slc_stack_nanned,
@@ -368,23 +366,31 @@ def _get_avg_ps(
     )
 
 
-def _get_maxes(arr, row_looks, col_looks):
+def _get_max_idxs(arr, row_looks, col_looks):
+    """Get the indices of the maximum value in each look window."""
     if row_looks == 1 and col_looks == 1:
         # No need to pad if we're not looking
         return np.where(arr == arr)
-    # Get the max value in each look window
-    # Start by dithering the array to avoid ties in the max
-    arr = arr.copy() + 1e-5 * np.random.rand(*arr.shape)
-    max_nums = take_looks(
-        arr, row_looks, col_looks, func_type="nanmax", edge_strategy="pad"
-    )
-    out_rows, out_cols = np.array(arr.shape) // [row_looks, col_looks]
-    # Make sure it's the output size we want:
-    max_nums = max_nums[:out_rows, :out_cols]
-
-    # Repeat the max values to back to the original size
-    maxes_filled = upsample_nearest(
-        max_nums, arr.shape[-2:], looks=(row_looks, col_looks)
-    )
-    # Find the indices of the max values in the original image
-    return np.where(maxes_filled == arr)
+    # Adjusted from this answer to not take every moving window
+    # https://stackoverflow.com/a/72742009/4174466
+    windows = np.lib.stride_tricks.sliding_window_view(arr, (row_looks, col_looks))[
+        ::row_looks, ::col_looks
+    ]
+    maxvals = np.nanmax(windows, axis=(2, 3))
+    indx = np.array((windows == np.expand_dims(maxvals, axis=(2, 3))).nonzero())
+
+    # In [82]: (windows == np.expand_dims(maxvals, axis = (2, 3))).nonzero()
+    # This gives 4 arrays:
+    # First two are the window indices
+    # (array([0, 0, 0, 1, 1, 1]),
+    # array([0, 1, 2, 0, 1, 2]),
+    # last two are the relative indices (within each window)
+    # array([0, 0, 1, 1, 1, 1]),
+    # array([1, 1, 1, 1, 1, 0]))
+    window_positions, relative_positions = indx.reshape((2, 2, -1))
+    # Multiply the first two by the window size to get the absolute indices
+    # of the top lefts of the windows
+    window_offsets = np.array([row_looks, col_looks]).reshape((2, 1))
+    # Then add the last two to get the relative indices
+    rows, cols = relative_positions + window_positions * window_offsets
+    return rows, cols
diff --git a/src/dolphin/ps.py b/src/dolphin/ps.py
index 1ffcd453..6099cf17 100644
--- a/src/dolphin/ps.py
+++ b/src/dolphin/ps.py
@@ -24,8 +24,8 @@ def create_ps(
     output_file: Filename,
     amp_mean_file: Filename,
     amp_dispersion_file: Filename,
-    amp_dispersion_threshold: float = 0.35,
-    max_ram_gb: float = 1.0,
+    amp_dispersion_threshold: float = 0.25,
+    block_size_gb: float = 1.0,
 ):
     """Create the amplitude dispersion, mean, and PS files.
 
@@ -40,8 +40,8 @@ def create_ps(
     amp_mean_file : Filename
         The output mean amplitude file.
     amp_dispersion_threshold : float, optional
-        The threshold for the amplitude dispersion. Default is 0.35.
-    max_ram_gb : int, optional
+        The threshold for the amplitude dispersion. Default is 0.25.
+    block_size_gb : int, optional
         The maximum amount of data to read at a time (in GB).
         Default is 1.0 GB.
     """
@@ -60,7 +60,7 @@ def create_ps(
         )
 
     vrt_stack = VRTStack.from_vrt_file(slc_vrt_file)
-    max_bytes = 1e9 * max_ram_gb
+    max_bytes = 1e9 * block_size_gb
     block_shape = vrt_stack._get_block_shape(max_bytes=max_bytes)
 
     # Initialize the intermediate arrays for the calculation
@@ -100,7 +100,7 @@ def create_ps(
 
 def calc_ps_block(
     stack_mag: ArrayLike,
-    amp_dispersion_threshold: float = 0.35,
+    amp_dispersion_threshold: float = 0.25,
     min_count: Optional[int] = None,
 ):
     r"""Calculate the amplitude dispersion for a block of data.
@@ -121,7 +121,7 @@ def calc_ps_block(
     amp_dispersion_threshold : float, optional
         The threshold for the amplitude dispersion to label a pixel as a PS:
             ps = amp_disp < amp_dispersion_threshold
-        Default is 0.35.
+        Default is 0.25.
     min_count : int, optional
         The minimum number of valid pixels to calculate the mean and standard
         deviation. If the number of valid pixels is less than `min_count`,
diff --git a/src/dolphin/stack.py b/src/dolphin/stack.py
index 4249d907..c0fc6c13 100755
--- a/src/dolphin/stack.py
+++ b/src/dolphin/stack.py
@@ -416,12 +416,8 @@ def iter_blocks(
         yield from self._loader.iter_blocks()
 
     def _get_block_shape(self, max_bytes=DEFAULT_BLOCK_BYTES):
-        test_file = self._get_non_vrt_file(self._gdal_file_strings[0])
-
         return io.get_max_block_shape(
-            # Note that we're using the actual first file, not the VRT
-            # since the VRT always has the same block size.
-            test_file,
+            self._gdal_file_strings[0],
             len(self),
             max_bytes=max_bytes,
         )
@@ -439,20 +435,6 @@ def _get_nodata_mask(self, nodata=nan, buffer_pixels=100):
                 buffer_pixels=buffer_pixels,
             )
 
-    @staticmethod
-    def _get_non_vrt_file(filename: Filename):
-        """Get one of the files within a VRT.
-
-        If the file is not a VRT, return the file itself.
-        Will traverse nested VRTs.
-        """
-        if Path(filename).suffix == ".vrt":
-            file_list = gdal.Info(fspath(filename), format="json")["files"]
-            if len(file_list) <= 1:
-                raise ValueError(f"VRT file {filename} contains no files")
-            return VRTStack._get_non_vrt_file(file_list[1])
-        return filename
-
     @property
     def shape(self):
         """Get the 3D shape of the stack."""
diff --git a/src/dolphin/template_derived.vrt b/src/dolphin/template_derived.vrt
deleted file mode 100644
index c368692e..00000000
--- a/src/dolphin/template_derived.vrt
+++ /dev/null
@@ -1,20 +0,0 @@
-<VRTDataset rasterXSize="{xsize}" rasterYSize="{ysize}">
-    <VRTRasterBand dataType="CFloat32" band="1" subClass="VRTDerivedRasterBand">
-        <PixelFunctionType>adjust</PixelFunctionType>
-        <PixelFunctionLanguage>Python</PixelFunctionLanguage>
-        <PixelFunctionCode><![CDATA[
-            import numpy as np
-            def adjust(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
-                       raster_ysize, buf_radius, gt, **kwargs):
-
-                out_ar[:] = np.abs(in_ar[0]) * np.exp(1j * (np.angle(in_ar[0]) + np.angle(in_ar[1])))
-            ]]>
-        </PixelFunctionCode>
-        <SimpleSource>
-            <SourceFilename relativeToVRT="1">{file1}</SourceFilename>
-        </SimpleSource>
-        <SimpleSource>
-            <SourceFilename relativeToVRT="1">{file2}</SourceFilename>
-        </SimpleSource>
-    </VRTRasterBand>
-</VRTDataset>
diff --git a/src/dolphin/unwrap.py b/src/dolphin/unwrap.py
index 5bc367a4..e9980c1a 100644
--- a/src/dolphin/unwrap.py
+++ b/src/dolphin/unwrap.py
@@ -7,11 +7,11 @@
 
 import numpy as np
 from osgeo import gdal
-from tqdm import tqdm
 
 from dolphin._log import get_log, log_runtime
 from dolphin._types import Filename
 from dolphin.io import copy_projection, get_raster_xysize
+from dolphin.utils import progress
 
 logger = get_log(__name__)
 
@@ -278,8 +278,8 @@ def run(
             )
             for inf, outf in zip(in_files, out_files)
         ]
-        for idx, fut in enumerate(tqdm(as_completed(futures)), start=1):
-            fut.result()
-            tqdm.write("Done with {} / {}".format(idx, len(futures)))
+        with progress() as p:
+            for fut in p.track(as_completed(futures)):
+                fut.result()
 
     return all_out_files
diff --git a/src/dolphin/utils.py b/src/dolphin/utils.py
index addad27e..1e60bdc0 100644
--- a/src/dolphin/utils.py
+++ b/src/dolphin/utils.py
@@ -1,5 +1,7 @@
 import datetime
 import re
+import resource
+import sys
 import warnings
 from os import fspath
 from pathlib import Path
@@ -8,6 +10,7 @@
 import numpy as np
 from numpy.typing import DTypeLike
 from osgeo import gdal, gdal_array, gdalconst
+from rich.progress import MofNCompleteColumn, Progress, SpinnerColumn, TimeElapsedColumn
 
 from dolphin._log import get_log
 from dolphin._types import Filename
@@ -16,6 +19,24 @@
 logger = get_log(__name__)
 
 
+def progress():
+    """Create a Progress bar context manager.
+
+    Usage
+    -----
+    >>> with progress() as p:
+    ...     for i in p.track(range(10)):
+    ...         pass
+    10/10 Working... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00
+    """
+    return Progress(
+        SpinnerColumn(),
+        MofNCompleteColumn(),
+        *Progress.get_default_columns()[:-1],  # Skip the ETA column
+        TimeElapsedColumn(),
+    )
+
+
 def numpy_to_gdal_type(np_dtype: DTypeLike) -> int:
     """Convert numpy dtype to gdal type.
 
@@ -82,8 +103,6 @@ def get_dates(filename: Filename, fmt: str = "%Y%m%d") -> List[datetime.date]:
     pattern = _date_format_to_regex(fmt)
     date_list = re.findall(pattern, path.stem)
     if not date_list:
-        msg = f"{filename} does not contain date like {fmt}"
-        logger.warning(msg)
         return []
     return [_parse_date(d, fmt) for d in date_list]
 
@@ -486,3 +505,78 @@ def upsample_nearest(
     arr_out = xp.zeros(shape=shape, dtype=arr.dtype)
     arr_out[..., :out_r, :out_c] = arr_up[..., :out_r, :out_c]
     return arr_out
+
+
+def get_max_memory_usage(units: str = "GB", children: bool = True) -> float:
+    """Get the maximum memory usage of the current process.
+
+    Parameters
+    ----------
+    units : str, optional, choices=["GB", "MB", "KB", "byte"]
+        The units to return, by default "GB".
+    children : bool, optional
+        Whether to include the memory usage of child processes, by default True
+
+    Returns
+    -------
+    float
+        The maximum memory usage in the specified units.
+
+    Raises
+    ------
+    ValueError
+        If the units are not recognized.
+
+    References
+    ----------
+    1. https://stackoverflow.com/a/7669279/4174466
+    2. https://unix.stackexchange.com/a/30941/295194
+    3. https://manpages.debian.org/bullseye/manpages-dev/getrusage.2.en.html
+
+    """
+    max_mem = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
+    if children:
+        max_mem += resource.getrusage(resource.RUSAGE_CHILDREN).ru_maxrss
+    if units.lower().startswith("g"):
+        factor = 1e9
+    elif units.lower().startswith("m"):
+        factor = 1e6
+    elif units.lower().startswith("k"):
+        factor = 1e3
+    elif units.lower().startswith("byte"):
+        factor = 1.0
+    else:
+        raise ValueError(f"Unknown units: {units}")
+    if sys.platform.startswith("linux"):
+        # on linux, ru_maxrss is in kilobytes, while on mac, ru_maxrss is in bytes
+        factor /= 1e3
+
+    return max_mem / factor
+
+
+def get_gpu_memory(pid: Optional[int] = None, gpu_id: int = 0) -> float:
+    """Get the memory usage (in GiB) of the GPU for the current pid."""
+    try:
+        from pynvml.smi import nvidia_smi
+    except ImportError:
+        raise ImportError("Please install pynvml through pip or conda")
+
+    def get_mem(process):
+        used_mem = process["used_memory"] if process else 0
+        if process["unit"] == "MiB":
+            multiplier = 1 / 1024
+        else:
+            logger.warning(f"Unknown unit: {process['unit']}")
+        return used_mem * multiplier
+
+    nvsmi = nvidia_smi.getInstance()
+    processes = nvsmi.DeviceQuery()["gpu"][gpu_id]["processes"]
+    if not processes:
+        return 0.0
+
+    if pid is None:
+        # Return sum of all processes
+        return sum(get_mem(p) for p in processes)
+    else:
+        procs = [p for p in processes if p["pid"] == pid]
+        return get_mem(procs[0]) if procs else 0.0
diff --git a/src/dolphin/workflows/_config_cli.py b/src/dolphin/workflows/_config_cli.py
index 7fcc46bc..33b57c69 100644
--- a/src/dolphin/workflows/_config_cli.py
+++ b/src/dolphin/workflows/_config_cli.py
@@ -14,10 +14,11 @@ def create_config(
     subdataset: Optional[str] = None,
     mask_files: Optional[List[str]] = None,
     ministack_size: Optional[int] = 15,
-    amp_dispersion_threshold: float = 0.35,
+    amp_dispersion_threshold: float = 0.25,
     strides: Tuple[int, int],
-    max_ram_gb: float = 1,
+    block_size_gb: float = 1,
     n_workers: int = 16,
+    threads_per_worker: int = 1,
     no_gpu: bool = False,
     single_update: bool = False,
 ):
@@ -53,8 +54,9 @@ def create_config(
             amp_dispersion_threshold=amp_dispersion_threshold,
         ),
         worker_settings=dict(
-            max_ram_gb=max_ram_gb,
+            block_size_gb=block_size_gb,
             n_workers=n_workers,
+            threads_per_worker=threads_per_worker,
             gpu_enabled=(not no_gpu),
         ),
     )
@@ -119,7 +121,7 @@ def get_parser(subparser=None, subcommand_name="run"):
     ps_group.add_argument(
         "--amp-dispersion-threshold",
         type=float,
-        default=0.35,
+        default=0.25,
         help="Threshold for the amplitude dispersion.",
     )
 
@@ -150,10 +152,10 @@ def get_parser(subparser=None, subcommand_name="run"):
         help="Disable the GPU (if using a machine that has one available).",
     )
     worker_group.add_argument(
-        "--max-ram-gb",
+        "--block-size-gb",
         type=float,
         default=1,
-        help="Maximum amount of RAM to use per worker.",
+        help="Size (in GB) of blocks of data to load at once time.",
     )
     worker_group.add_argument(
         "--n-workers",
@@ -161,6 +163,12 @@ def get_parser(subparser=None, subcommand_name="run"):
         default=16,
         help="Number of workers to use.",
     )
+    worker_group.add_argument(
+        "--threads-per-worker",
+        type=int,
+        default=1,
+        help="Number of threads to use per worker.",
+    )
     # parser.add_argument(
     #     "--mask-files",
     #     nargs=argparse.ZERO_OR_MORE,
diff --git a/src/dolphin/workflows/_run_cli.py b/src/dolphin/workflows/_run_cli.py
index f790b597..ec003fce 100644
--- a/src/dolphin/workflows/_run_cli.py
+++ b/src/dolphin/workflows/_run_cli.py
@@ -2,13 +2,9 @@
 import argparse
 from typing import Optional
 
-from dolphin._log import log_runtime
+from dolphin._log import get_log
 
-from . import s1_disp
-from .config import Workflow
 
-
-@log_runtime
 def run(config_file: str, debug: bool = False, log_file: Optional[str] = None):
     """Run the displacement workflow.
 
@@ -21,11 +17,30 @@ def run(config_file: str, debug: bool = False, log_file: Optional[str] = None):
     log_file : str, optional
         If provided, will log to this file in addition to stderr.
     """
+    from threadpoolctl import ThreadpoolController
+
+    from dolphin.utils import get_max_memory_usage
+
+    from . import s1_disp
+    from .config import Workflow
+
+    # Set the logging level for all `dolphin.` modules
+    logger = get_log("dolphin", debug=debug)
+
+    # Set the environment variables for the workers
+    # TODO: Is this the best place to do this?
     cfg = Workflow.from_yaml(config_file)
     cfg.create_dir_tree(debug=debug)
 
+    controller = ThreadpoolController()
+    controller.limit(limits=cfg.worker_settings.threads_per_worker)
+
     s1_disp.run(cfg, debug=debug, log_file=log_file)
 
+    # Print the maximum memory usage for each worker
+    max_mem = get_max_memory_usage(units="GB")
+    logger.info(f"Maximum memory usage: {max_mem:.2f} GB")
+
 
 def get_parser(subparser=None, subcommand_name="run"):
     """Set up the command line interface."""
@@ -60,6 +75,7 @@ def main(args=None):
     """Get the command line arguments and run the workflow."""
     parser = get_parser()
     parsed_args = parser.parse_args(args)
+
     run(parsed_args.config_file, debug=parsed_args.debug)
 
 
diff --git a/src/dolphin/workflows/_utils.py b/src/dolphin/workflows/_utils.py
index 45036c49..9a4906da 100644
--- a/src/dolphin/workflows/_utils.py
+++ b/src/dolphin/workflows/_utils.py
@@ -90,6 +90,7 @@ def setup_output_folder(
     start_idx: int = 0,
     strides: Dict[str, int] = {"y": 1, "x": 1},
     creation_options: Optional[List] = None,
+    nodata: Optional[float] = 0,
     output_folder: Optional[Path] = None,
 ) -> List[Path]:
     """Create empty output files for each band after `start_idx` in `vrt_stack`.
@@ -115,6 +116,8 @@ def setup_output_folder(
         [dolphin.io.compute_out_shape][]
     creation_options : list, optional
         List of options to pass to the GDAL driver, by default None
+    nodata : float, optional
+        Nodata value to use for the output files, by default 0.
     output_folder : Path, optional
         Path to output folder, by default None
         If None, will use the same folder as the first SLC in `vrt_stack`
@@ -150,6 +153,7 @@ def setup_output_folder(
             nbands=1,
             dtype=dtype,
             strides=strides,
+            nodata=nodata,
             options=creation_options,
         )
 
diff --git a/src/dolphin/workflows/config.py b/src/dolphin/workflows/config.py
index ca568e10..ba5d5f63 100644
--- a/src/dolphin/workflows/config.py
+++ b/src/dolphin/workflows/config.py
@@ -78,7 +78,7 @@ class PsOptions(BaseModel):
     )
 
     amp_dispersion_threshold: float = Field(
-        0.35,
+        0.25,
         description="Amplitude dispersion threshold to consider a pixel a PS.",
         gt=0.0,
     )
@@ -240,9 +240,17 @@ class WorkerSettings(BaseSettings):
     n_workers: int = Field(
         16, ge=1, description="Number of cpu cores to use for processing (if CPU)"
     )
-    max_ram_gb: float = Field(
+    threads_per_worker: int = Field(
+        4,
+        ge=1,
+        description=(
+            "Number of threads to use per worker. This sets the OMP_NUM_THREADS"
+            " environment variable."
+        ),
+    )
+    block_size_gb: float = Field(
         1.0,
-        description="Maximum RAM (in GB) to use for processing",
+        description="Size (in GB) of blocks of data to load at once time.",
         gt=0.1,
     )
 
diff --git a/src/dolphin/workflows/sequential.py b/src/dolphin/workflows/sequential.py
index 88fa2eb0..a8c0d02c 100644
--- a/src/dolphin/workflows/sequential.py
+++ b/src/dolphin/workflows/sequential.py
@@ -148,7 +148,10 @@ def run_evd_sequential(
 
     logger.info(f"Running EVD on compressed files: {adjustment_vrt_stack}")
     adjusted_comp_slc_files = setup_output_folder(
-        adjustment_vrt_stack, driver="GTiff", strides=strides
+        adjustment_vrt_stack,
+        driver="GTiff",
+        strides=strides,
+        nodata=0,
     )
 
     writer = io.Writer()
@@ -165,7 +168,7 @@ def run_evd_sequential(
         logger.debug(msg)
 
         # Run the phase linking process on the current adjustment stack
-        cur_mle_stack, tcorr = run_mle(
+        cur_mle_stack, _ = run_mle(
             cur_data,
             half_window=half_window,
             strides=strides,
@@ -177,6 +180,7 @@ def run_evd_sequential(
             n_workers=n_workers,
             gpu_enabled=gpu_enabled,
         )
+        np.nan_to_num(cur_mle_stack, copy=False)
 
         # Get the location within the output file, shrinking down the slices
         out_row_start = rows.start // ys
diff --git a/src/dolphin/workflows/single.py b/src/dolphin/workflows/single.py
index 4856560d..0382ccbe 100644
--- a/src/dolphin/workflows/single.py
+++ b/src/dolphin/workflows/single.py
@@ -91,6 +91,7 @@ def run_evd_single(
         f"{vrt}: from {Path(vrt.file_list[first_non_comp_idx]).name} to"
         f" {Path(vrt.file_list[-1]).name}"
     )
+    logger.info(f"Total stack size (in pixels): {vrt.shape}")
     # Set up the output folder with empty files to write into
     output_slc_files = setup_output_folder(
         vrt,
@@ -98,6 +99,7 @@ def run_evd_single(
         start_idx=first_non_comp_idx,
         strides=strides,
         output_folder=output_folder,
+        nodata=0,
     )
 
     # Create the empty compressed SLC file
@@ -107,6 +109,7 @@ def run_evd_single(
         like_filename=vrt.outfile,
         output_name=comp_slc_file,
         nbands=1,
+        nodata=0,
         # Note that the compressed SLC is the same size as the original SLC
     )
 
@@ -119,6 +122,7 @@ def run_evd_single(
         nbands=1,
         dtype=np.float32,
         strides=strides,
+        nodata=0,
     )
 
     # Iterate over the stack in blocks
@@ -159,9 +163,20 @@ def run_evd_single(
         except PhaseLinkRuntimeError as e:
             # note: this is a warning instead of info, since it should
             # get caught at the "skip_empty" step
-            logger.warning(f"Exception at ({rows}, {cols}): {e}")
+            msg = f"At block {rows.start}, {cols.start}: {e}"
+            if "are all NaNs" in e.args[0]:
+                # Some SLCs in the ministack are all NaNs
+                # This happens from a shifting burst window near the edges,
+                # and seems to cause no issues
+                logger.info(msg)
+            else:
+                logger.warning(msg)
             continue
 
+        # Fill in the nan values with 0
+        np.nan_to_num(cur_mle_stack, copy=False)
+        np.nan_to_num(tcorr, copy=False)
+
         # Save each of the MLE estimates (ignoring the compressed SLCs)
         assert len(cur_mle_stack[first_non_comp_idx:]) == len(output_slc_files)
         # Get the location within the output file, shrinking down the slices
@@ -182,7 +197,6 @@ def run_evd_single(
         # TODO: make a flag? We don't always need to save the compressed SLCs
         writer.queue_write(cur_comp_slc, comp_slc_file, out_row_start, out_col_start)
         # logger.debug(f"Saved compressed block SLC to {cur_comp_slc_file}")
-        # tqdm.write(" Finished block, loading next block.")
 
     # Block until all the writers for this ministack have finished
     logger.info(f"Waiting to write {writer.num_queued} blocks of data.")
diff --git a/src/dolphin/workflows/wrapped_phase.py b/src/dolphin/workflows/wrapped_phase.py
index e5808641..b6c25026 100644
--- a/src/dolphin/workflows/wrapped_phase.py
+++ b/src/dolphin/workflows/wrapped_phase.py
@@ -55,7 +55,7 @@ def run(cfg: Workflow, debug: bool = False) -> List[VRTInterferogram]:
             amp_mean_file=cfg.ps_options.amp_mean_file,
             amp_dispersion_file=cfg.ps_options.amp_dispersion_file,
             amp_dispersion_threshold=cfg.ps_options.amp_dispersion_threshold,
-            max_ram_gb=cfg.worker_settings.max_ram_gb,
+            block_size_gb=cfg.worker_settings.block_size_gb,
         )
 
     # #########################
@@ -63,6 +63,9 @@ def run(cfg: Workflow, debug: bool = False) -> List[VRTInterferogram]:
     # #########################
     pl_path = cfg.phase_linking.directory
 
+    from dolphin._background import NvidiaMemoryWatcher
+
+    watcher = NvidiaMemoryWatcher() if cfg.worker_settings.gpu_enabled else None
     existing_files = list(pl_path.glob("*.tif"))
     if len(existing_files) > 0:
         logger.info(f"Skipping EVD step, {len(existing_files)} files already exist")
@@ -77,7 +80,7 @@ def run(cfg: Workflow, debug: bool = False) -> List[VRTInterferogram]:
                 reference_idx=0,
                 # mask_file=cfg.inputs.mask_file,
                 ps_mask_file=cfg.ps_options.output_file,
-                max_bytes=cfg.worker_settings.max_ram_gb * 1e9,
+                max_bytes=cfg.worker_settings.block_size_gb * 1e9,
                 n_workers=cfg.worker_settings.n_workers,
                 gpu_enabled=cfg.worker_settings.gpu_enabled,
                 beta=cfg.phase_linking.beta,
@@ -91,11 +94,13 @@ def run(cfg: Workflow, debug: bool = False) -> List[VRTInterferogram]:
                 ministack_size=cfg.phase_linking.ministack_size,
                 # mask_file=cfg.inputs.mask_file,
                 ps_mask_file=cfg.ps_options.output_file,
-                max_bytes=cfg.worker_settings.max_ram_gb * 1e9,
+                max_bytes=cfg.worker_settings.block_size_gb * 1e9,
                 n_workers=cfg.worker_settings.n_workers,
                 gpu_enabled=cfg.worker_settings.gpu_enabled,
                 beta=cfg.phase_linking.beta,
             )
+    if watcher:
+        watcher.notify_finished()
 
     # ###################################################
     # 4. Form interferograms from estimated wrapped phase


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 9c651d18..8a3aec28 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -33,9 +33,10 @@ jobs:
               pydantic=1.10
               pymp-pypi=0.4.5
               pyproj=3.2
+              rich=12.0
               gdal=3.3
               scipy=1.5
-              tqdm=4.60
+              threadpoolctl>=3.0
 
       fail-fast: false
     name: ${{ matrix.os.label }} • ${{ matrix.deps.label }}
diff --git a/tests/test_workflows_config.py b/tests/test_workflows_config.py
index 8bebd409..48ba462e 100644
--- a/tests/test_workflows_config.py
+++ b/tests/test_workflows_config.py
@@ -34,7 +34,7 @@ def test_ps_options_defaults(tmpdir):
     # Change directory so the creation of the default files doesn't fail
     with tmpdir.as_cwd():
         pso = config.PsOptions()
-        assert pso.amp_dispersion_threshold == 0.35
+        assert pso.amp_dispersion_threshold == 0.25
         assert pso.directory == Path("PS")
         assert pso.output_file == Path("PS/ps_pixels.tif")
         assert pso.amp_dispersion_file == Path("PS/amp_dispersion.tif")
@@ -94,7 +94,7 @@ def test_worker_settings_defaults():
     assert ws.gpu_enabled is True
     assert ws.gpu_id == 0
     assert ws.n_workers == 16
-    assert ws.max_ram_gb == 1.0
+    assert ws.block_size_gb == 1.0
 
 
 def test_worker_env_defaults(monkeypatch):
@@ -124,9 +124,9 @@ def test_worker_env_defaults(monkeypatch):
     ws = config.WorkerSettings()
     assert ws.n_workers == 8
 
-    monkeypatch.setenv("DOLPHIN_MAX_RAM_GB", "4.5")
+    monkeypatch.setenv("DOLPHIN_BLOCK_SIZE_GB", "4.5")
     ws = config.WorkerSettings()
-    assert ws.max_ram_gb == 4.5
+    assert ws.block_size_gb == 4.5
 
 
 @pytest.fixture()