-
Notifications
You must be signed in to change notification settings - Fork 23
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
Update channel_combination.py #48
base: master
Are you sure you want to change the base?
Conversation
Added module RAMRS_SVD_weighting, which performs SVD channel combination on the raw data (prior to averaging). Useful in cases where there might be drift across the averages, which can be corrected after channel combination more easily, since channel combination increases the SNR. Also added a few additional outputs to the svd_weighting module to return just the channel weights (without phase correction). This is useful for cases where channel weights are derived from the unsuppressed H2O signal, and then applied to the suppressed H2O signal. The phase shift term won't be correct for the suppressed H2O signal.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi Laura,
It seems like we have some significant differences of opinion about how the library should be structured. Do we need to have a tcon to discuss? I hate to always be criticising your PRs because I do really appreciate your work on the project. I will be submitting my PRs to you for review as well now that you are more active on the project again.
Best,
Ben
Hi Ben - No worries. You are probably right about us having a bit of a difference of opinion on this : ) That's OK, though, and since you are the primary author on this and have far more time to spend on taking care of it, I'm glad to defer to you. Realistically, I only have about 3 weeks worth of hours to work on this before this project completely shuts down forever, so I'm just trying to put as much of the stuff I've developed over the years out there before then, in case we never get more funding. Anyway, take a look at some of my comments so you can get a sense of my thought process and why I structured things this way, and let me know how to proceed.
-Laura
|
||
return channel_weights * numpy.exp(-1j * phase_shift) / numpy.sum(numpy.abs(channel_weights)) | ||
return channel_weights_phase_corrected, channel_weights, phase_shift |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think that this is unnecessary - the phase shift here is completely insignificant as there is going to be some proper phase correction happening later on anyway, this is just a very rough and ready attempt to put the signal in phase to make it easier to see what is going on when plotting e.g. to check for frequency drift. The point is that the channel weights are only defined up to an arbitrary phase, so we have to choose something and this makes the most sense - it is actually better than doing nothing in all cases.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, I'm fine with changing that. I realize that there is additional phase correction coming down the line. The only issue is that I like to pull my channel weights off of the unsuppressed H2O signal since the SNR is so much higher, and then apply them to the H2O suppressed signal. The phase correction for one doesn't really work for the other, so it ends up making the H2O suppressed signal's phase a bit worse in some cases. This will be easily corrected later on since it's only 0th order anyway, so it's not that big of a deal. I can change it back if that would be easier.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am surprised that the unsuppressed water phase is not correct for the suppressed signal, they should have an exactly matched phase, except for the suppressed water peak, and who cares about that anyway. It should also be extremely simple to modify the phase of the suppressed signal after using the wrongly phased weights, simply by making the first point of the FID real, which is all this does anyway. I would really prefer to keep this as simple as possible as a function, which means returning just one thing.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also there is no reason that the "unphased" weights from the wref are any better as an option for the phase of the suppressed signal.
@@ -19,3 +24,125 @@ def combine_channels(data, weights=None): | |||
weighted_data = weights.reshape((len(weights), 1)) * data | |||
combined_data = weighted_data.sum(axis=0) | |||
return combined_data | |||
|
|||
|
|||
def RAMRS_SVD_weighting(data, phase_correct = True, whiten = False, W = None, noise_range = None): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Whitening the channel noise is something that I have been meaning to do for a long time, so very pleased to see that has made an appearance. However, I am not sure this is the best way to set up the interface. My general approach throughout has been to provide small, modular functions to perform each step of the processing, rather than to give a few monolithic functions with lots of options in the parameter string, like whiten, phase_correct etc. This is easier to understand, easier to test and easier to customise. I would much prefer code like
whitened_wref = whiten(wref_data)
channel_weights = svd_weighting(whitened_wref)
combined_data = np.average(wsup_data, axis=1, weights=channel_weights)
rather than
combined_data = RAMRS_SVD_weighting(wsup_data, phase_correct=True, whiten=True, W = channel_weights)
""" | ||
This channel combination approach operates on raw MRS (spectral) data, | ||
prior to averaging repetitions from multiple acquisitions. It is useful in | ||
cases where there is instability in frequency and phase shifts across |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is actually the same method of calculating the channel weights, as you call the existing function. Essentially the major difference is that this function expects a 3d array instead of 2d. I think a far better solution would be to modify the existing svd_weighting()
function so that it accepts an axis
parameter to define which axis the channels run over, then internally it reshapes all other axes into a single line and calculates the channel weights that way. That would then be a more seamless improvement of the existing code.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
My original idea was to do it exactly like you described - but I thought instead it might be a good idea to actually create a function that specifically references this paper, and makes it explicit that it is replicating the method described therein. It becomes sort of self-explanatory in that way, and people can follow along with the paper and see exactly what is going on. So often I'll see an algorithm implemented in a paper, try to find some code that replicates it somewhere, end up cobbling something together from bits and pieces of other libraries, but I'm never quite sure if it's done exactly right. I've worked out the kinks on this approach and tested it so extensively that I thought it would be good to present it in this kind of monolithic way (which it definitely is!) and use the other functions from Suspect internally as needed. An alternative would be to put all of this into an example script somewhere, but then people will have to slog through it, if they find it at all. Anyway, that was my thought process on the whole thing. I'm glad to change it to the method you described, though. I already have a version of that somewhere in my code base.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Interesting that we both thought about this in the same way and came to opposite conclusions. I am trying quite hard to maintain an easy to understand set of example notebooks on readthedocs, there is already a notebook on channel combination that looks at different ways of weighting the channels including the SVD method. My plan was to either add another notebook later in the series to cover noise whitening, or to go back and incorporate that into the existing notebook. The great thing about the notebooks is that they can combine explanation and plots inline with the code, so they help with understanding the concepts much better than simple comments. Although your docstring is excellent, most end users aren't going to dig through the codebase to understand how the algorithm works, they just want some functions they can call.
I think that a chunk of the code inside this function could do with being extracted into a whiten()
function which can be called separately. Ultimately I think that the implementation of this function should be essentially just calling about 3 other functions from the module, and that may make it superfluous. If you don't have the time to work on that, then I am happy to take it on.
Added module RAMRS_SVD_weighting (Receive Array MRS SVD weighting, from the reference in the comments). Performs SVD channel combination on the raw data (prior to averaging). Useful in cases where there might be drift across the averages, which can be corrected after channel combination more easily, since channel combination increases the SNR of the averages.
Added a few additional outputs to the svd_weighting module to also return just the channel weights (without phase correction). This is useful for cases where channel weights are derived from the unsuppressed H2O signal, and then applied to the suppressed H2O signal. The phase shift term won't be correct for the suppressed H2O signal.
Note: I know the updates to the svd_weighting algorithm will break the way it is used in the channel combination Jupyter notebook, though, so we can discuss alternative ways to implement this if it is too cumbersome. It would be helpful to be able to return the channel weights, minus the phase correction term, though.