Allow for list of int. Only project active times #1689

domfournier · 2025-08-06T17:48:01Z

(cherry picked from commit 1dfafa5)

Summary

This proposed change streamlines the compute of predicted data from time domain fields. Instead of looping over all dts, only loop over the times that needed by the projection of the receivers.

PR Checklist

If this is a work in progress PR, set as a Draft PR
Linted my code according to the style guides.
Added tests to verify changes to the code.
Added necessary documentation to any new functions/classes following the
expect style.
Marked as ready for review (if this is was a draft PR), and converted
to a Pull Request
Tagged @simpeg/simpeg-developers when ready for review.

Reference issue

What does this implement/fix?

Additional information

(cherry picked from commit 1dfafa5)

codecov · 2025-08-06T18:34:55Z

Codecov Report

❌ Patch coverage is 77.77778% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.24%. Comparing base (cba34e9) to head (94aafe1).

Files with missing lines	Patch %	Lines
simpeg/electromagnetics/time_domain/receivers.py	75.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1689      +/-   ##
==========================================
- Coverage   85.79%   85.24%   -0.56%     
==========================================
  Files         427      427              
  Lines       55506    55509       +3     
  Branches     5189     5188       -1     
==========================================
- Hits        47623    47320     -303     
- Misses       6466     6736     +270     
- Partials     1417     1453      +36

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

lheagy

Thanks for this @domfournier! This is a nice update. I ran it for a simple example and am seeing about a 5-10% speed-up, which is meaningful!

**Before: **

After:

It looks like some of the EM tests are failing. I can take a look. I see the pr is still a draft. Is there something else you would like to update before bringing it in? or just the tests?

lheagy · 2025-09-01T16:38:07Z

It looks like it is falling down inside of evalDeriv. If instead of calling the P from the receiver, we construct the full P, that should do the trick. Alternatively (and this will likely be faster), we can just grab the correct indices in v for the forward and project into a larger vector with zeros in the adjoint. We would need to get the active_times which we should be able to do without constructing the full Pt matrix

simpeg/simpeg/electromagnetics/time_domain/receivers.py

Lines 178 to 208 in d58b57f

    
               def evalDeriv(self, src, mesh, time_mesh, f, v, adjoint=False): 
        
                   """Derivative of projected fields with respect to the inversion model times a vector. 
        
                   Parameters 
        
                   ---------- 
        
                   src : simpeg.electromagnetics.frequency_domain.sources.BaseTDEMSrc 
        
                       A time-domain EM source 
        
                   mesh : discretize.base.BaseMesh 
        
                       The mesh on which the discrete set of equations is solved 
        
                   time_mesh : discretize.TensorMesh 
        
                       A 1D ``TensorMesh`` defining the time discretization 
        
                   f : simpeg.electromagnetic.time_domain.fields.FieldsTDEM 
        
                       The solution for the fields defined on the mesh 
        
                   v : numpy.ndarray 
        
                       A vector 
        
                   adjoint : bool, default = ``False`` 
        
                       If ``True``, return the adjoint 
        
                   Returns 
        
                   ------- 
        
                   numpy.ndarray 
        
                       derivative of fields times a vector projected to the receiver(s) 
        
                   """ 
        
                   P = self.getP(mesh, time_mesh, f) 
        
                   if not adjoint: 
        
                       return P * v 
        
                   elif adjoint: 
        
                       # dP_dF_T = P.T * v #[src, self] 
        
                       # newshape = (len(dP_dF_T)/time_mesh.nN, time_mesh.nN ) 
        
                       return P.T * v  # np.reshape(dP_dF_T, newshape, order='F')

Ps = self.getSpatialP(mesh, f)
Pt = self.getTimeP(time_mesh, f)
P = sp.kron(Pt, Ps)

lheagy · 2025-09-01T20:55:49Z

simpeg/electromagnetics/time_domain/receivers.py

+        numpy.ndarray
+            Active times for the receiver.
+        """
+        return np.unique(sp.find(projection)[1])


I think we should be able to do the following since it is a sparse matrix

Suggested change

return np.unique(sp.find(projection)[1])

return projection.indices

I suspect that using .indices is going to be faster. Just two caveats:

We still need to do np.unique(projection.indices). If we don't get the unique ones, the Pt matrix after doing Pt = Pt[:, self.active_times(Pt)] will be actually larger than the original one (in case we have repeated column indices).

Using indices or the find function is equivalent only if projection is a csr_matrix. For example, if it's a csc_matrix, then the .indices will be row indices instead of the column indices. Since getP always return a csr_matrix, this should not be a problem, But I would clarify this in the docstring, and maybe checking isinstance(projection, csr_matrix) would be a good idea.

santisoler

Leaving a few comments below. I'll try to experiment with my second thought and share results.

santisoler · 2025-09-02T18:16:29Z

simpeg/electromagnetics/time_domain/receivers.py


        Ps = self.getSpatialP(mesh, f)
        Pt = self.getTimeP(time_mesh, f)
+        Pt = Pt[:, self.active_times(Pt)]


The one issue I see by adding this line in getP is that any time we want to compute the dot product between P and any vector, the vector should also be "reduced" to only the active times. This makes the getP method less standalone: we need to call active_times again if we want to compute dot products of P with any vector.

Would it be possible to apply the optimization directly in eval and evalDeriv, without modifying the output of getP? ~~I guess I'm asking if filtering out non-active times in the P matrix (P = sp.kron(Pt, Ps)) is possible and would achieve same performance as doing so in Pt.~~ Update: no, it wouldn't be possible because the output of kron is a coo_matrix that cannot be subscripted.

santisoler · 2025-09-03T20:10:43Z

Sharing here a few ideas we discussed today in the meeting:

It might be a good idea to make the P-related methods private. We can tackle it in another PR, no need to do it here.
We should fix the evalDeriv methods so we they can compute the dot product of P with any vector without running into the shape mismatch error.
One alternative we discussed is making the output of getP a LinearOperator, so through the matvec and rmatvec it can handle the optimization. The methods themselves can apply the reduction both to the Ps matrix and to the vector, and then cast it into an output vector, filling the zeros.

Allow for list of int. Only project active times

fe3c824

(cherry picked from commit 1dfafa5)

domfournier requested review from lheagy and santisoler August 6, 2025 17:49

Fix reference to timeP

94aafe1

domfournier marked this pull request as draft August 6, 2025 19:30

lheagy approved these changes Sep 1, 2025 •

edited

Loading

View reviewed changes

lheagy reviewed Sep 1, 2025

View reviewed changes

santisoler reviewed Sep 2, 2025

View reviewed changes

santisoler mentioned this pull request Sep 25, 2025

Make ComplexMap.deriv to return sparse matrix or 1D array #1687

Closed

6 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Allow for list of int. Only project active times #1689

Allow for list of int. Only project active times #1689

Uh oh!

domfournier commented Aug 6, 2025

Uh oh!

codecov bot commented Aug 6, 2025 •

edited

Loading

Uh oh!

lheagy left a comment

Uh oh!

lheagy commented Sep 1, 2025 •

edited

Loading

Uh oh!

lheagy Sep 1, 2025 •

edited

Loading

Uh oh!

santisoler Sep 2, 2025

Uh oh!

santisoler left a comment

Uh oh!

santisoler Sep 2, 2025 •

edited

Loading

Uh oh!

santisoler commented Sep 3, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

	return np.unique(sp.find(projection)[1])
	return projection.indices

Allow for list of int. Only project active times #1689

Are you sure you want to change the base?

Allow for list of int. Only project active times #1689

Uh oh!

Conversation

domfournier commented Aug 6, 2025

Summary

PR Checklist

Reference issue

What does this implement/fix?

Additional information

Uh oh!

codecov bot commented Aug 6, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Uh oh!

lheagy left a comment

Choose a reason for hiding this comment

Uh oh!

lheagy commented Sep 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lheagy Sep 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

santisoler Sep 2, 2025

Choose a reason for hiding this comment

Uh oh!

santisoler left a comment

Choose a reason for hiding this comment

Uh oh!

santisoler Sep 2, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

santisoler commented Sep 3, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

codecov bot commented Aug 6, 2025 •

edited

Loading

lheagy commented Sep 1, 2025 •

edited

Loading

lheagy Sep 1, 2025 •

edited

Loading

santisoler Sep 2, 2025 •

edited

Loading