I will have a look later whether this performs like what I wrote yesterday (which is how @sarangnemo implemented it in NutMEG back then) - I guess @larsoner got so excited by the plots I shared that he had to try by himself! 😄

For what it's worth, the code I wrote up yesterday performs close to a scalar beamformer with optimal orientation picking - I have some preliminary code to also allow orientations to vary across time, which might make things more interesting, which I'd be happy to contribute.

I guess @larsoner got so excited by the plots I shared that he had to try by himself!

Yes sorry I got excited about figuring out how to take care of the complex-valued case... :)

I have some preliminary code to also allow orientations to vary across time, which might make things more interesting, which I'd be happy to contribute.

Eventually we can allow this by allowing project_onto take directions of shape n_vertices, 3, n_times as well. I'm not sure if we want to allow it by default, though. At least for the extreme case of doing SVD per time point, one direction will capture 100% of the variance, and the result will be equivalent to just np.linalg.norm(stc.data, axis=1) (i.e., stc.magnitude()) at that point.

In between these two extremes -- i.e., not doing SVD across all time points, or SVD for each individual time point -- things will vary a bit, but I'd be tempted to have people .crop and project_onto and recombine if they want to try this sort of thing. Before adding it in MNE we'd want to look to see the extent to which it's useful in practice. The SVD-across-all-time approach here already makes sense because it enables things like envelope correlation computation with free-orientation (typically volumetric) minimum norm inverse solutions, which is not really doable currently in master (and we also know we'll need it for some beamformer implementations, as you say).

So our methods agree - left image is my version (minus polarity flipping), right image is your version !

(Data is cropped toy data, so not necessarily meaningful output, just to compare numbers.)

Excellent -- in that case this is good to go @agramfort . Once this is merged I'll rebase the beamformer refactoring PR and we can properly test the two vector "unit-noise-gain" methods with and without SVD afterward

should directions be a public parameter?

It was easy enough to add and validation was simple. People could want to for example project onto the max power orientations they got from some other filter. So I don't see a problem/maintenance cost from keeping it.

we have stc.normal so should we rather of stc.principal_component?

Starting from scratch I'd actually prefer the API to be stc.project_onto('normal' | 'pca') or so. I can rename it to principal_component but if we think of other ways to do the projections our namespace gets more annoying -- I'd rather keep and potentially in the future add additional options to a generic project_onto function...

ok I propose to deprecate VectorSourceEstimate.normal and introduce

VectorSourceEstimate.project(directions, src=None, use_cps=True, return_directions=False)

where directions can be
'normal' -> requires src
'pca' -> does not require src? and use return_directions if you want to get the estimated directions
array of shape (n_dipoles, 3)

sounds like a plan?

I'd rather just always return the directions, return_ kwargs I'd rather leave as a mechanism for backward compatibility. It's very easy to do [0] to get the first output of a function.

src is mandatory in both cases to make the sign flip directions not arbitrary

(deal with sign flip ambiguity in pca)

Yep, I'll push shortly

Thx @larsoner