Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

JP-3653: Nsclean speedup #8547

Merged
merged 5 commits into from
Jun 13, 2024
Merged
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
7 changes: 7 additions & 0 deletions CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,13 @@ master_background_mos
``slit.source_name`` now contains the string "BKG" instead of
"background". [#8533]

nsclean
-------

- Improved run time of NSClean.fit() by using a vector rather than a large,
sparse matrix to perform linear algebra operations with a diagonal weight
matrix. [#8547]

outlier_detection
-----------------

Expand Down
9 changes: 6 additions & 3 deletions jwst/nsclean/lib.py
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,10 @@

# Get data and weights for this line
d = data[y][self.mask[y]] # unmasked (useable) data
p = np.diag(self.P[y][self.mask[y]]) # Weights
# The line below uses a vector to represent a diagonal weight matrix.
# Multiplications by this vector later on may be viewed as
# equivalent formulations to multiplication by diag(p).
p = self.P[y][self.mask[y]] # Weights

Check warning on line 140 in jwst/nsclean/lib.py

View check run for this annotation

Codecov / codecov/patch

jwst/nsclean/lib.py#L140

Added line #L140 was not covered by tests

# If none of the pixels in this line is useable (all masked out),
# skip and move on to the next line.
Expand Down Expand Up @@ -164,15 +167,15 @@

# Compute the Moore-Penrose inverse of A = P*B.
# $A^+ = (A^H A)^{-1} A^H$
A = np.matmul(p, B)
A = B*p[:, np.newaxis]

Check warning on line 170 in jwst/nsclean/lib.py

View check run for this annotation

Codecov / codecov/patch

jwst/nsclean/lib.py#L170

Added line #L170 was not covered by tests
AH = np.conjugate(A.transpose()) # Hermitian transpose of A
pinv_PB = np.matmul(np.linalg.inv(np.matmul(AH, A)), AH)

# Solve for the Fourier transform of this line's background samples.
# The way that we have done it, this multiplies the input data by the
# number of samples used for the fit.
rfft = np.zeros(self.nx//2 + 1, dtype=np.complex64)
rfft[:k.shape[1]] = np.matmul(np.matmul(pinv_PB, p), d)
rfft[:k.shape[1]] = np.matmul(pinv_PB, p*d)

Check warning on line 178 in jwst/nsclean/lib.py

View check run for this annotation

Codecov / codecov/patch

jwst/nsclean/lib.py#L178

Added line #L178 was not covered by tests

# Numpy requires that the forward transform multiply
# the data by n. Correct normalization.
Expand Down