11# -*- coding: utf-8 -*-
22import os
3+ from functools import lru_cache
34import numpy as np
45from numpy import ones , kron , mean , eye , hstack , tile
56from numpy .linalg import pinv
@@ -86,67 +87,61 @@ def _list_outputs(self):
8687 return outputs
8788
8889
89- def ICC_rep_anova (Y , nocache = False ):
90+ @lru_cache (maxsize = 1 )
91+ def ICC_projection_matrix (shape ):
92+ nb_subjects , nb_conditions = shape
93+ 94+ x = kron (eye (nb_conditions ), ones ((nb_subjects , 1 ))) # sessions
95+ x0 = tile (eye (nb_subjects ), (nb_conditions , 1 )) # subjects
96+ X = hstack ([x , x0 ])
97+ return X @ pinv (X .T @ X , hermitian = True ) @ X .T
98+ 99+ 100+ def ICC_rep_anova (Y , projection_matrix = None ):
90101 """
91102 the data Y are entered as a 'table' ie subjects are in rows and repeated
92103 measures in columns
104+
93105 One Sample Repeated measure ANOVA
106+
94107 Y = XB + E with X = [FaTor / Subjects]
95108
96- This is a hacked up (but fully compatible) version of ICC_rep_anova
97- from nipype that caches some very expensive operations that depend
98- only on the input array shape - if you're going to run the routine
99- multiple times (like, on every voxel of an image), this gives you a
100- HUGE speed boost for large input arrays. If you change the dimensions
101- of Y, it will reinitialize automatially. Set nocache to True to get
102- the original, much slower behavior. No, actually, don't do that.
109+ ``ICC_rep_anova`` involves an expensive operation to compute a projection
110+ matrix, which depends only on the shape of ``Y``, which is computed by
111+ calling ``ICC_projection_matrix(Y.shape)``. If arrays of multiple shapes are
112+ expected, it may be worth pre-computing and passing directly as an
113+ argument to ``ICC_rep_anova``.
114+
115+ If only one ``Y.shape`` will occur, you do not need to explicitly handle
116+ these, as the most recently calculated matrix is cached automatically.
117+ For example, if you are running the same computation on every voxel of
118+ an image, you will see significant speedups.
119+
120+ If a ``Y`` is passed with a new shape, a new matrix will be calculated
121+ automatically.
103122 """
104- global icc_inited
105- global current_Y_shape
106- global dfc , dfe , dfr
107- global nb_subjects , nb_conditions
108- global x , x0 , X
109- global centerbit
110- 111- try :
112- current_Y_shape
113- if nocache or (current_Y_shape != Y .shape ):
114- icc_inited = False
115- except NameError :
116- icc_inited = False
117- 118- if not icc_inited :
119- [nb_subjects , nb_conditions ] = Y .shape
120- current_Y_shape = Y .shape
121- dfc = nb_conditions - 1
122- dfe = (nb_subjects - 1 ) * dfc
123- dfr = nb_subjects - 1
123+ [nb_subjects , nb_conditions ] = Y .shape
124+ dfc = nb_conditions - 1
125+ dfr = nb_subjects - 1
126+ dfe = dfr * dfc
124127
125128 # Compute the repeated measure effect
126129 # ------------------------------------
127130
128131 # Sum Square Total
129- mean_Y = mean (Y )
130- SST = ((Y - mean_Y ) ** 2 ).sum ()
131- 132- # create the design matrix for the different levels
133- if not icc_inited :
134- x = kron (eye (nb_conditions ), ones ((nb_subjects , 1 ))) # sessions
135- x0 = tile (eye (nb_subjects ), (nb_conditions , 1 )) # subjects
136- X = hstack ([x , x0 ])
137- centerbit = X @ pinv (X .T @ X , hermitian = True ) @ X .T
132+ demeaned_Y = Y - mean (Y )
133+ SST = np .sum (demeaned_Y ** 2 )
138134
139135 # Sum Square Error
140- predicted_Y = centerbit @ Y .flatten ("F" )
141- residuals = Y .flatten ("F" ) - predicted_Y
142- SSE = (residuals ** 2 ).sum ()
143- 144- residuals .shape = Y .shape
136+ if projection_matrix is None :
137+ projection_matrix = ICC_projection_matrix (Y .shape )
138+ residuals = Y .flatten ("F" ) - (projection_matrix @ Y .flatten ("F" ))
139+ SSE = np .sum (residuals ** 2 )
145140
146141 MSE = SSE / dfe
147142
148143 # Sum square session effect - between columns/sessions
149- SSC = (( mean (Y , 0 ) - mean_Y ) ** 2 ). sum ( ) * nb_subjects
144+ SSC = np . sum ( mean (demeaned_Y , 0 ) ** 2 ) * nb_subjects
150145 MSC = SSC / dfc / nb_subjects
151146
152147 session_effect_F = MSC / MSE
@@ -162,6 +157,4 @@ def ICC_rep_anova(Y, nocache=False):
162157 e_var = MSE # variance of error
163158 r_var = (MSR - MSE ) / nb_conditions # variance between subjects
164159
165- icc_inited = True
166- 167160 return ICC , r_var , e_var , session_effect_F , dfc , dfe
0 commit comments