adjusted numerical issues in time indexing

ibs-lab · Apr 4, 2024 · 7e11171 · 7e11171
1 parent 983e73e
commit 7e11171
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 9 deletions.
diff --git a/src/cedalion/glm/design_matrix.py b/src/cedalion/glm/design_matrix.py
@@ -130,7 +130,10 @@ def make_hrf_regressors(
     dt = 1 / t.cd.sampling_rate
     n_pre = round(trange[0] / dt)
     n_post = round(trange[1] / dt)
-    t_hrf = np.arange(n_pre * dt, (n_post + 1) * dt, dt)
+    # np.arange results can be non-consistent when using non-integer steps
+    #t_hrf = np.arange(n_pre * dt, (n_post + 1) * dt, dt)   
+    # using linspace instead
+    t_hrf = np.linspace(trange[0], trange[1], abs(n_post) + abs(n_pre) + 1)
     nt = len(t)
 
     # prune good stim
@@ -140,6 +143,11 @@ def make_hrf_regressors(
     onset = np.zeros((nt, n_cond))
     n_trials = np.zeros(n_cond)
 
+    # handle case of single condition
+    if n_cond == 1:
+        cond_names = np.array([cond_names])
+        s = np.expand_dims(s, axis=1)
+
     for i_cond in range(n_cond):
         lst_t = np.where(s[:, lst_cond[i_cond]] == 1)[0]
         lstp = np.where((lst_t + n_pre >= 0) & (lst_t + n_post < nt))[0]

diff --git a/src/cedalion/glm/solve.py b/src/cedalion/glm/solve.py
@@ -224,6 +224,8 @@ def get_HRFs(
         hrfs: xarray.DataArray containing HRFs for every condition and every channel (time x chromo x channels x conditions)
     """
 
+    dt = 1 / predicted_hrf.cd.sampling_rate
+
     unit = predicted_hrf.pint.units
     predicted_hrf = predicted_hrf.pint.dequantify()
 
@@ -235,12 +237,20 @@ def get_HRFs(
     )
 
     # get time axis for HRFs:
-    time_hrf = predicted_hrf.sel(
-        time=slice(
-            stim_onsets.sel(condition=conds[0]) + HRFmin,
-            stim_onsets.sel(condition=conds[0]) + HRFmax,
-        ),
-    ).time - stim_onsets.sel(condition=conds[0])
+
+    # time_hrf = predicted_hrf.sel(
+    #    time=slice(
+    #        stim_onsets.sel(condition=conds[0]) + HRFmin,
+    #        stim_onsets.sel(condition=conds[0]) + HRFmax,
+    #    ),
+    # ).time - stim_onsets.sel(condition=conds[0])
+
+    n_pre = round(HRFmin / dt)
+    n_post = round(HRFmax / dt)
+    # np.arange results can be non-consistent when using non-integer steps
+    # t_hrf = np.arange(n_pre * dt, (n_post + 1) * dt, dt)
+    # using linspace instead
+    time_hrf = np.linspace(HRFmin, HRFmax, abs(n_post) + abs(n_pre) + 1)
 
     hrfs = xr.DataArray(
         np.zeros(
@@ -266,8 +276,12 @@ def get_HRFs(
             hrf = predicted_hrf.sel(
                 chromo=chromo,
                 time=slice(
-                    stim_onsets.sel(condition=cond) + HRFmin,
-                    stim_onsets.sel(condition=cond) + HRFmax,
+                    predicted_hrf.time.sel(
+                        time=stim_onsets.sel(condition=cond) + HRFmin, method="nearest"
+                    ),
+                    predicted_hrf.time.sel(
+                        time=stim_onsets.sel(condition=cond) + HRFmax, method="nearest"
+                    ),
                 ),
             )
             # change time axis to relative time