[ENH, FIX] Add wavelet denoising and fix t2smap workflow #90

tsalo · 2018-07-16T00:02:57Z

Closes #80.

Changes proposed in this pull request:

Re-implements wavelet denoising from an older version of MEICA, with minor changes to work with pywt's current API.
- Adds wvpca argument to tedana.decomposition.eigendecomp.tedpca, tedana.workflows.tedana, and the tedana workflow CLI.
Fixes the t2smap workflow, which hadn't been updated to work with the recent changes to the package.
- Also adds a new argument, fitmode, which can fit the mono exponential model per timepoint if desired.

I have *no* clue if this works properly or not. Taken from https://github.com/ME-ICA/me-ica/blob/da6ac4c23a59145530a16bd6f25c676460 fe9436/meica.libs/tedana.py#L514-L518 and adjusted to work with pywt.idwt’s arguments (i.e., dropping correct_size argument).

Hopefully it works.

codecov · 2018-07-16T17:08:22Z

Codecov Report

Merging #90 into master will increase coverage by 14.79%.
The diff coverage is 86.55%.

@@             Coverage Diff             @@
##           master      #90       +/-   ##
===========================================
+ Coverage   27.95%   42.74%   +14.79%     
===========================================
  Files          25       27        +2     
  Lines        1381     1516      +135     
===========================================
+ Hits          386      648      +262     
+ Misses        995      868      -127

Impacted Files	Coverage Δ
tedana/info.py	`100% <ø> (ø)`	⬆️
tedana/tests/test_tedana.py	`0% <ø> (ø)`	⬆️
tedana/model/monoexponential.py	`95.74% <ø> (ø)`	⬆️
tedana/cli/run_tedana.py	`0% <0%> (ø)`	⬆️
tedana/model/fit.py	`8.08% <0%> (ø)`	⬆️
tedana/utils/io.py	`13.59% <0%> (ø)`	⬆️
tedana/workflows/tedana.py	`11.39% <0%> (+11.39%)`	⬆️
tedana/tests/test_model_monoexponential.py	`100% <100%> (ø)`	⬆️
tedana/tests/test_decomposition_utils.py	`100% <100%> (ø)`
tedana/tests/test_t2smap.py	`100% <100%> (ø)`
... and 18 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ec80e16...5e51b26. Read the comment docs.

emdupre

This is great, thank you ! I'm assuming this is still WIP from the travis errors ?

I had a few questions about variable names, but just wanted to say how great the addition is -- and the improved docstrings are a nice touch !

I just patched the circle error (finally) -- any chance you could merge in master, again ? Sorry for the trouble !

emdupre · 2018-07-16T20:44:04Z

tedana/decomposition/_utils.py

+    -------
+    mmix_wt : :obj:`numpy.ndarray`
+        Wavelet-transformed data.
+    cAlen : :obj:`int`


Is there a reason to call it cAlen here and cAl in idwtmat() ?

I just copied the functions from the old version of tedana and removed the correct_size argument in idwtmat, which pywt must have dropped at some point.

emdupre · 2018-07-16T20:49:01Z

tedana/decomposition/_utils.py

+    llt = len(np.hstack(pywt.dwt(mmix[0], 'db2')))
+    mmix_wt = np.zeros([mmix.shape[0], llt])
+    for ii in range(mmix_wt.shape[0]):
+        wtx = pywt.dwt(mmix[ii], 'db2')


From the pywavelets documentation, it seems that wtx would be the approximation coefficients for the specified row of mmix. Then cAlen would be the length of the first entry of wtx in the last for-loop.

If you agree with that assessment, could we improve these variable names ?

I think it's the approximation and detail coefficients, not that I'm sure what those are. cAlen is the number of coefficients (i.e., len of cA). Oh, what about n_coefs_approx instead of cAlen, coefs_tup instead of wtx, n_coefs_total instead of llt, mmix_coefs instead of mmix_wt?

What does mmix stand for again?

emdupre · 2018-07-16T20:49:24Z

tedana/decomposition/_utils.py

+    ----------
+    mmix_wt : :obj:`numpy.ndarray`
+        Wavelet-transformed data.
+    cAl : :obj:`int`


See previous comments ^

emdupre · 2018-07-16T20:51:52Z

tedana/workflows/t2smap.py


    LGR.info('Computing optimal combination')
-    tsoc = np.array(model.make_optcom(catd, t2s, tes, mask, combmode),


Can you explain these changes for me ?

Most of these changes are just taken from the current version of the tedana workflow, except for changing the output variable names from the fit_decay functions, which I did because I think the new ones are more informative. I also dropped the part where it writes out the t2ss and s0vs maps because fit_decay_ts doesn't return those and I honestly have no clue what they are.

Are there any changes in particular that don't make sense?

Yes, sorry, that wasn't a helpful question on my part ! Your explanation solved it anyway, though: I was confused why we were renaming tsoc to OCcatd.

emdupre · 2018-07-16T20:55:32Z

Also, it looks like this may partially address #25 ?

Also use more informative variable names in wavelet transformation functions.

tsalo · 2018-07-17T02:17:00Z

The tests are now passing, and there are several working (though basic) tests for the t2smap workflow. Yay increased coverage!

I guess it does address #25 a bit, though my changes are mostly to the workflow rather than the fitting functions.

I'm not sure whether the wvpca works or not, but other than that I think it looks good. WDYT?

tsalo · 2018-07-17T02:23:37Z

Oh and I did notice two situations where model.optcom fails:

If the limited T2* map is used (how it was written, but not how tedana is) instead of the full one (combmode='t2s' and fitmode='all', both of which are the defaults)
If the full T2* timeseries (how it's hardcoded) is used with combmode='ste' and fitmode='ts'

emdupre · 2018-07-17T14:05:14Z

This is great ! 🎉

Do you think we should have some testing for the wvpca functions (i.e., dwtmat and idwtmat) ? I know that they're pulled almost directly from pywt – so I'm not sure how much sense a unit test makes – but now is the time I'd prefer to add them, if we're going to !

Regarding failing options: I had no idea the @pytest.mark.xfail decorator existed ! But more importantly: can you explain the first failure a little more for me ? What do you mean "not how tedana is" ?

tsalo · 2018-07-17T14:34:37Z

I can definitely add some simple tests. Something to make sure that the data are returned in the right shape at least. I don't know enough about wavelets to come up with anything better atm.

Honestly, I had to hunt down a way to allow a test to fail. I didn't want to hide this weird condition that was failing but didn't want to have to deal with a continually breaking test, which would throw off the PRs once merged.

By "not how tedana is" I meant that t2smap was using the limited t2s map for optcom before, but the tedana workflow was using the full t2smap. So I guess how t2smap was set up originally, it would fail (at least on the test data in the repository). These failures could all be entirely dependent on the data for all I know though.

emdupre · 2018-07-17T15:39:47Z

Do you think the failure has anything to do with this ? I can't find if we ever resolved that.

tsalo · 2018-07-17T16:32:42Z

We continued talking about the differences somewhere (maybe on Gitter?), but I can't find the discussion. Anyway, I don't think it could be the source of the failure, because my current version uses start_echo=1 in t2smap to match the tedana workflow.

I just realized that ste doesn't use the t2s estimates (of course), so dealing with the alpha weights for the averaging as if they're 3D (SxExT) (as one would have when using fit_decay_ts instead of fit_decay) instead of 2D (SxE) (the default) would cause an error.

I am changing how make_optcom decides to handle the alpha weights, but I also think we should change how make_optcom is organized overall, because t2s and combmode are kind of redundant. At the very least, t2s should be an optional argument, since it isn't used when combmode = 'ste'. I'm going to commit my proposed changes so you can see what I was thinking, and then we can discuss them in the code comments.

tsalo · 2018-07-17T16:35:58Z

tedana/model/combine.py

@@ -8,26 +8,27 @@
 LGR = logging.getLogger(__name__)


-def make_optcom(data, t2s, tes, mask, combmode, verbose=True):
+def make_optcom(data, tes, mask, t2s=None, combmode='t2s', verbose=True):


@emdupre I think we could probably consolidate these into one argument, but in the interest of clarity two is fine. If we move t2s so it's a keyword argument, then it can default to None (which works when combmode is ste but not when it's t2s). Then below we check the inputs to make sure they make sense. WDYT?

Just to clarify: what would happen when combmode is ste if we make t2s the one (keyword) argument ? Or am I misunderstanding the suggestion ?

Maybe combmode = [t2s array] and combmode = 'ste' could be the two options? I admit, that's not the best suggestion... Or I guess when t2s = None, ste is used by default? So it would be def make_optcom(data, tes, mask, t2s=None, verbose=True): and if t2s is an array, then t2s combination is used, but if it's None, then ste combination is used?

That makes sense to me, unless we foresee adding other combmode options.. @handwerkerd might have thoughts on this !

emdupre · 2018-07-17T17:31:07Z

I had the same thought, but I can't find it anywhere ! I was having this same issue with the discussion we had on #88...

I see your point about having t2s as a keyword argument -- I think your idea makes good sense ! I am, though, a little confused about the suggestion to reduce it down to one argument; see q above.

For now, I'm happy to get circle passing again, move forward with the two argument idea, and reorganize optcom in another PR. I just wasn't sure if we could nail down the pytest failure !

tsalo · 2018-07-17T17:35:19Z

Whoops, I missed one of the make_optcom calls. Hopefully circle will pass now.

emdupre · 2018-07-17T18:26:38Z

This LGTM ! I'll open a new issue re: a t2s keyword argument for make_optcom so we can continue the discussion ( / get @handwerkerd's thoughts) there.

Thanks again, @tsalo !!

tsalo added 9 commits May 25, 2018 13:12

Merge remote-tracking branch 'ME-ICA/master'

f8811e3

Merge remote-tracking branch 'ME-ICA/master'

d3c13e7

Merge remote-tracking branch 'ME-ICA/master'

51adcd3

Merge remote-tracking branch 'ME-ICA/master'

58ee2de

Fix typo.

271262d

Update t2smap workflow.

fea2911

Hopefully it works.

Update t2smap CLI.

f448954

Update requirements.

ee541e7

emdupre reviewed Jul 16, 2018

View reviewed changes

tsalo added 8 commits July 16, 2018 17:12

Merge remote-tracking branch 'ME-ICA/master' into wavelet

9856609

Add generated files to workflow docstrings

5e9fa0c

Also use more informative variable names in wavelet transformation functions.

Fix typos in t2smap workflow.

c8c271a

Add t2smap workflow test.

7b96912

Fix t2smap tests.

4469617

Run t2smap workflow in travis.

b774ccf

Get t2smap tests working.

ed1668c

Fix tests maybe?

87ea146

tsalo added 2 commits July 17, 2018 10:42

Add wavelet transform test and squish bug!

eee78cc

Of course pywt is "Python Webtrends connector"

fe7fe1d

Resolve issues with 'ste' and volume-wise T2* in make_optcom.

8bcd585

tsalo commented Jul 17, 2018

View reviewed changes

Fix style issue.

641d1e6

Fix make_optcom call in fitmodels_direct.

5e51b26

emdupre merged commit 712e0bb into ME-ICA:master Jul 17, 2018

emdupre mentioned this pull request Jul 17, 2018

Revisit make_optcom #93

Closed

tsalo mentioned this pull request Jul 24, 2018

Refactor t2smap #25

Closed

tsalo mentioned this pull request Nov 20, 2020

Examine differences in results between original MEICA & current tedana #622

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[ENH, FIX] Add wavelet denoising and fix t2smap workflow #90

[ENH, FIX] Add wavelet denoising and fix t2smap workflow #90

tsalo commented Jul 16, 2018

codecov bot commented Jul 16, 2018 •

edited

Loading

emdupre left a comment

emdupre Jul 16, 2018

tsalo Jul 16, 2018

emdupre Jul 16, 2018

tsalo Jul 16, 2018

emdupre Jul 16, 2018

emdupre Jul 16, 2018

tsalo Jul 16, 2018

emdupre Jul 17, 2018

emdupre commented Jul 16, 2018

tsalo commented Jul 17, 2018

tsalo commented Jul 17, 2018

emdupre commented Jul 17, 2018

tsalo commented Jul 17, 2018

emdupre commented Jul 17, 2018

tsalo commented Jul 17, 2018

tsalo Jul 17, 2018

emdupre Jul 17, 2018

tsalo Jul 17, 2018

emdupre Jul 17, 2018

emdupre commented Jul 17, 2018

tsalo commented Jul 17, 2018

emdupre commented Jul 17, 2018


		LGR.info('Computing optimal combination')
		tsoc = np.array(model.make_optcom(catd, t2s, tes, mask, combmode),

[ENH, FIX] Add wavelet denoising and fix t2smap workflow #90

[ENH, FIX] Add wavelet denoising and fix t2smap workflow #90

Conversation

tsalo commented Jul 16, 2018

codecov bot commented Jul 16, 2018 • edited Loading

Codecov Report

emdupre left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

emdupre commented Jul 16, 2018

tsalo commented Jul 17, 2018

tsalo commented Jul 17, 2018

emdupre commented Jul 17, 2018

tsalo commented Jul 17, 2018

emdupre commented Jul 17, 2018

tsalo commented Jul 17, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

emdupre commented Jul 17, 2018

tsalo commented Jul 17, 2018

emdupre commented Jul 17, 2018

codecov bot commented Jul 16, 2018 •

edited

Loading