diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index c362deeb..557e46a2 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -25,7 +25,7 @@ repos:
             #      - id: check-docstring-first
       - id: check-case-conflict # Check for files with names that would conflict on a case-insensitive filesystem
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.16.2
+    rev: v0.16.5
     hooks:
       - id: ruff-format
       - id: ruff-check
diff --git a/goals/exfit_upgrade_summary.md b/goals/exfit_upgrade_summary.md
index 88d7d262..9d2f7d98 100644
--- a/goals/exfit_upgrade_summary.md
+++ b/goals/exfit_upgrade_summary.md
@@ -48,9 +48,9 @@ Pyy built once per pixel.  Total complexity O(n) per pixel.
 `run()` now accepts a `workers` keyword (mirrors `LineRatioFit` API):
 
 ```python
-fit.run(components=2, workers=-1)   # all CPUs
-fit.run(components=2, workers=4)    # 4 processes
-fit.run(components=2)               # serial (default)
+fit.run(components=2, workers=-1)  # all CPUs
+fit.run(components=2, workers=4)  # 4 processes
+fit.run(components=2)  # serial (default)
 ```
 
 Implementation uses `ProcessPoolExecutor` with:
diff --git a/goals/lmfit_sparse_jac_bug.md b/goals/lmfit_sparse_jac_bug.md
index db0e5619..7255ccd8 100644
--- a/goals/lmfit_sparse_jac_bug.md
+++ b/goals/lmfit_sparse_jac_bug.md
@@ -26,7 +26,7 @@ element-wise multiplication on a non-square sparse matrix:
 ```python
 # minimizer.py line 1595–1596
 if issparse(ret.jac):
-    hess = (ret.jac.T * ret.jac).toarray()   # BUG: * is element-wise in scipy >= 1.9
+    hess = (ret.jac.T * ret.jac).toarray()  # BUG: * is element-wise in scipy >= 1.9
 ```
 
 `ret.jac` has shape `(m, n)` where `m` = number of residuals and `n` = number of
@@ -50,14 +50,16 @@ from lmfit import Minimizer, Parameters
 x_data = np.array([1.0, 2.0, 3.0, 4.0])
 y_data = np.array([2.1, 4.0, 5.9, 8.1])
 
+
 def residual(params):
-    a = params['a'].value
-    b = params['b'].value
-    return a * x_data + b - y_data   # 4 residuals
+    a = params["a"].value
+    b = params["b"].value
+    return a * x_data + b - y_data  # 4 residuals
+
 
 params = Parameters()
-params.add('a', value=1.0)
-params.add('b', value=0.0)
+params.add("a", value=1.0)
+params.add("b", value=0.0)
 
 # Jacobian sparsity: 4 residuals × 2 parameters, all entries non-zero
 sparsity = lil_matrix((4, 2), dtype=np.int8)
@@ -65,7 +67,7 @@ sparsity[:, :] = 1
 sparsity = sparsity.tocsr()
 
 mini = Minimizer(residual, params)
-result = mini.minimize(method='least_squares', jac_sparsity=sparsity)
+result = mini.minimize(method="least_squares", jac_sparsity=sparsity)
 # ValueError: inconsistent shapes (2, 4) and (4, 2)
 ```
 
@@ -133,11 +135,11 @@ from scipy.optimize import least_squares
 from scipy.sparse import lil_matrix
 import numpy as np
 
-result = least_squares(fun, x0, bounds=bounds, jac_sparsity=sparsity, method='trf')
+result = least_squares(fun, x0, bounds=bounds, jac_sparsity=sparsity, method="trf")
 
 # Manual covariance from Jacobian (works for both sparse and dense)
 J = result.jac
-if hasattr(J, 'toarray'):
+if hasattr(J, "toarray"):
     J = J.toarray()
 try:
     cov = np.linalg.inv(J.T @ J)
diff --git a/goals/perf_optimization_session_2026-03-06.md b/goals/perf_optimization_session_2026-03-06.md
index 4f92eea2..507b1153 100644
--- a/goals/perf_optimization_session_2026-03-06.md
+++ b/goals/perf_optimization_session_2026-03-06.md
@@ -79,8 +79,7 @@ Four Horsehead Nebula FITS files from `pdrtpy/testdata/`, read with `Measurement
 
 ```python
 self._observedratios_flat = {
-    k: (v.data.flatten(), v.uncertainty.array.flatten())
-    for k, v in self._observedratios.items()
+    k: (v.data.flatten(), v.uncertainty.array.flatten()) for k, v in self._observedratios.items()
 }
 ```
 
@@ -114,9 +113,7 @@ _qq = np.squeeze(np.reshape(residuals, newshape))
 
 # After
 modelpix_exp = modelpix.reshape((-1,) + (1,) * mdata.ndim)
-residuals_arr = ma.masked_invalid(
-    (mdata[np.newaxis, ...] - modelpix_exp) / merror[np.newaxis, ...]
-)
+residuals_arr = ma.masked_invalid((mdata[np.newaxis, ...] - modelpix_exp) / merror[np.newaxis, ...])
 _qq = np.squeeze(np.reshape(residuals_arr, newshape))
 ```
 
@@ -135,6 +132,7 @@ Note: test suite runtime also dropped from 367 s to 143 s.
 ```python
 _worker_model_interps = None  # per-process cache
 
+
 def _init_worker(model_points, model_values):
     """Build RegularGridInterpolators once per worker process."""
     global _worker_model_interps
@@ -143,8 +141,10 @@ def _init_worker(model_points, model_values):
         for pts, vals in zip(model_points, model_values)
     ]
 
-def _fit_pixel_worker(j, obs_data_j, obs_err_j, init_density, init_rf,
-                      minn, maxn, minfuv, maxfuv, nan_policy, minimize_kwargs):
+
+def _fit_pixel_worker(
+    j, obs_data_j, obs_err_j, init_density, init_rf, minn, maxn, minfuv, maxfuv, nan_policy, minimize_kwargs
+):
     """Fit a single pixel in a worker process. Returns (j, MinimizerResult)."""
     ...
 ```
@@ -212,7 +212,7 @@ Fits all N valid pixels in a **single** `scipy.optimize.least_squares` call:
 2. **Vectorized residual** using batched `RegularGridInterpolator` queries:
    ```python
    def _joint_residual(x):
-       pts = np.column_stack([x[::2], x[1::2]])          # (n_valid, 2)
+       pts = np.column_stack([x[::2], x[1::2]])  # (n_valid, 2)
        mvalues = np.array([interp(pts) for interp in interps])  # (n_ratios, n_valid)
        return ((obs_data - mvalues) / obs_err).flatten()  # (n_valid * n_ratios,)
    ```
@@ -220,16 +220,21 @@ Fits all N valid pixels in a **single** `scipy.optimize.least_squares` call:
 3. **Block-diagonal `jac_sparsity`**: pixel `j`'s parameters only affect pixel `j`'s residuals.
    ```python
    from scipy.sparse import lil_matrix
+
    sparsity = lil_matrix((N * n_ratios, 2 * N), dtype=np.int8)
    for j in range(N):
-       sparsity[j*n_ratios:(j+1)*n_ratios, 2*j:2*j+2] = 1
+       sparsity[j * n_ratios : (j + 1) * n_ratios, 2 * j : 2 * j + 2] = 1
    ```
 
 4. **Direct scipy call**:
    ```python
    result = scipy.optimize.least_squares(
-       _joint_residual, x0, bounds=(lb, ub),
-       jac_sparsity=sparsity.tocsr(), tr_solver='lsmr', method='trf',
+       _joint_residual,
+       x0,
+       bounds=(lb, ub),
+       jac_sparsity=sparsity.tocsr(),
+       tr_solver="lsmr",
+       method="trf",
    )
    ```
 
@@ -249,10 +254,10 @@ With ~5200 valid pixels, the global TRF `ftol` is satisfied in ~2 outer iteratio
 
 Extracted from the block-diagonal joint Jacobian:
 ```python
-J_i = jac[i*n_ratios:(i+1)*n_ratios, 2*i:2*i+2]  # (n_ratios, 2)
+J_i = jac[i * n_ratios : (i + 1) * n_ratios, 2 * i : 2 * i + 2]  # (n_ratios, 2)
 cov_i = inv(J_i.T @ J_i)
 stderr_density = sqrt(cov_i[0, 0])
-stderr_rf      = sqrt(cov_i[1, 1])
+stderr_rf = sqrt(cov_i[1, 1])
 ```
 
 ### Support for future regularization