Retrieval of yeo7 information has off-by-1 error in python code #347
Description
This came to my attention via the following error:
File /local/scratch/BrainStat/brainstat/stats/SLM.py:325, in SLM._surfstat_to_brainstat_rft(self)
323 yeo_names, _ = fetch_yeo_networks_metadata(7)
324 yeo_names.insert(0, "Undefined")
--> 325 yeo7_index = yeo7[self.P["peak"]["vertid"][i]]
326 if "yeo7" not in self.P["peak"]:
327 self.P["peak"]["yeo7"] = []
IndexError: index 64984 is out of bounds for axis 0 with size 64984
Full Error
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
Cell In[14], line 36
22 interact = {
23 "model": term_ses * term_panss + term_sub,
24 "contrast": (ses_post * df["panss-ratio"]) - (ses_pre * df["panss-ratio"]),
25 "label": "intr",
26 }
27 slm = SLM(
28 ses["model"],
29 ses["contrast"] * -1,
(...)
34 two_tailed=False,
35 )
---> 36 slm.fit(np.asanyarray(surface.sel(desc="thickness", smoothing=8)))
File /local/scratch/BrainStat/brainstat/stats/SLM.py:158, in SLM.fit(self, Y)
156 self._unmask()
157 if self.correction is not None:
--> 158 self.multiple_comparison_corrections(student_t_test)
File /local/scratch/BrainStat/brainstat/stats/SLM.py:254, in SLM.multiple_comparison_corrections(self, student_t_test)
251 self.P[field][field2] = [self.P[field][field2]]
252 self.Q = Q1
--> 254 self._surfstat_to_brainstat_rft()
File /local/scratch/BrainStat/brainstat/stats/SLM.py:325, in SLM._surfstat_to_brainstat_rft(self)
323 yeo_names, _ = fetch_yeo_networks_metadata(7)
324 yeo_names.insert(0, "Undefined")
--> 325 yeo7_index = yeo7[self.P["peak"]["vertid"][i]]
326 if "yeo7" not in self.P["peak"]:
327 self.P["peak"]["yeo7"] = []
IndexError: index 64984 is out of bounds for axis 0 with size 64984
From what I can tell, "clustid" and "vertid" are being saved as 1-based indices in brainstat.stats._multiple_comparisons::peak_clus, but being used at brainstat/stats/SLM.py:325 to directly index a numpy array, which uses 0-based indexing. In my example, I just happened to have a peak on the last voxel of the mesh, triggering the IndexError. But of course, this would mean potentially the wrong networks are being indexed even when the code runs without error.
I'm not familiar enough with the code to know if this problem creeps up elsewhere. For what it's worth, it's definitely surprising to me as a python user that the vertids are 1-based, although it's probably impossible to change this now. But I think at least this should be prominently documented.
Hopefully this is sufficient information, it's a bit tricky to make a minimal example, but I can try to pull one together if necessary.