RF: Optimize ICC_rep_anova with a memoized helper function

[effigies]

@bbfrederick Post #3453, I got motivated to look at your PR again. I grokked the code a little better and realized we could do this more Pythonically.

WDYT? I set the maxsize=1 to match your effective cache size. I don't really know how big these matrices can get, so I'm not sure whether that was a memory optimization or just an assumption that the number of switches back and forth between shapes is likely to be small.

@kristoferm94 Would be interested in your thoughts. This does break up the improved matrix ordering you provided, but I suspect the caching benefit is going to dominate here.

As some basic tests:

In [1]: import numpy as np

In [2]: rng = np.random.default_rng()

In [3]: from nipype.algorithms.icc import *

In [4]: %time ICC_rep_anova(rng.normal(size=(50, 10)))
CPU times: user 26.4 ms, sys: 67 ms, total: 93.4 ms
Wall time: 15.4 ms
Out[4]: 
(-0.0010955283641284907,
 -0.0012372157263673023,
 1.130569661036011,
 0.01640374953027528,
 9,
 441)

In [5]: ICC_projection_matrix.cache_clear()

In [6]: %time ICC_rep_anova(rng.normal(size=(50, 10)))
CPU times: user 431 µs, sys: 18.8 ms, total: 19.3 ms
Wall time: 3.14 ms
Out[6]: 
(0.052737606398966394,
 0.049204080066708754,
 0.8837938966422438,
 0.028379415383302194,
 9,
 441)

In [7]: %timeit ICC_rep_anova(rng.normal(size=(50, 10)))
160 µs ± 7.95 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Cells 4-6 is just to distinguish compilation effects. Cell 7 shows the effect of caching.

[codecov]

Codecov Report

Merging #3454 (ad2249a) into master (3dc858e) will increase coverage by 0.00%.
The diff coverage is 94.11%.

@@           Coverage Diff           @@
##           master    #3454   +/-   ##
=======================================
  Coverage   65.24%   65.24%           
=======================================
  Files         309      309           
  Lines       40833    40838    +5     
  Branches     5376     5377    +1     
=======================================
+ Hits        26641    26645    +4     
  Misses      13118    13118           
- Partials     1074     1075    +1     

Flag	Coverage Δ
unittests	65.02% <94.11%> (+<0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
nipype/algorithms/icc.py	58.97% <94.11%> (+1.44%)	⬆️

---

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 3dc858e...ad2249a. Read the comment docs.

[bbfrederick]

Sorry it took so long to get back to you - kept finding myself on mobile with no screen big enough to look at code.

Yes, that looks like a good switch - I wasn't aware of lru_cache (I need to get out more)! Maxsize=1 is fine - my use case is calculating it on ~260k voxels in a set of 3D images from 1000 subjects with 4 conditions - exactly the same matrix size each time, so the repeated pseudoinverse was killing me.... I don't imagine you'd need to switch the matrix up very frequently, although I'm sure you could come up with some pathological application that would.

[effigies]

Sounds good. I agree that there's probably not a huge advantage in storing more than one precomputed pseudoinverse. I'll go ahead and merge. If someone can come up with a relevant use case for increasing the cache, we can revisit.