Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .pre-commit-config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions goals/exfit_upgrade_summary.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
20 changes: 11 additions & 9 deletions goals/lmfit_sparse_jac_bug.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -50,22 +50,24 @@ 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)
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)
```

Expand Down Expand Up @@ -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)
Expand Down
31 changes: 18 additions & 13 deletions goals/perf_optimization_session_2026-03-06.md
Original file line number Diff line number Diff line change
Expand Up @@ -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()
}
```

Expand Down Expand Up @@ -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))
```

Expand All @@ -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
Expand All @@ -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)."""
...
```
Expand Down Expand Up @@ -212,24 +212,29 @@ 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,)
```

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",
)
```

Expand All @@ -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
Expand Down