non_tr_onsets.Rmd

---
jupyter:
  jupytext:
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.11.5
  kernelspec:
    display_name: Python 3 (ipykernel)
    language: python
    name: python3
---

# Using onsets that do not start on a TR

You’ll need the file `new_cond.txt` to type along with this
demonstration.

```{python}
import numpy as np
import matplotlib.pyplot as plt
import nibabel as nib
```

Imagine we are analyzing our `example image`.  It has a TR of 2.5, and 173 TRs.

```{python}
TR = 2.5
n_trs = 173
```

The actual condition file for this dataset is `ds114_sub009_t2r1_cond.txt`.

```{python}
# Load the function to fetch the data file we need.
import nipraxis
# Fetch the data file.
cond_fname = nipraxis.fetch_file('ds114_sub009_t2r1_cond.txt')
# Show the file name of the fetched data.
cond_fname
```

You may remember it has a *block design* with blocks of length 12 TRs while
the subject is doing the task.

What if we had a different *event related* condition file like the one at
`new_cond.txt`?

```{python}
# Fetch the data file.
new_cond_fname = nipraxis.fetch_file('new_cond.txt')
# Show the file name of the fetched data.
new_cond_fname
```

```{python}
cond_data = np.loadtxt(new_cond_fname)
cond_data
onsets_seconds = cond_data[:, 0]
durations_seconds = cond_data[:, 1]
amplitudes = cond_data[:, 2]
```

Notice that the onsets of the events can happen in the middle of the volumes
(well after the volumes have started).

```{python}
onsets_in_scans = onsets_seconds / TR
onsets_in_scans
```

Notice also that the events have *amplitudes* between 1 and 3.  The events of
amplitude 3 we expect to have an evoked brain response three times higher than
events with amplitude 1.

What to do about the events with onsets that don’t exactly align with the
start of the TRs (volumes)?

One option would be to round the event onsets to the nearest TR.  This will
mean that the event model will be different from our expected response by TR
seconds / 2 == 1.25 seconds in this case.

Can we do better than that?

We could make a neural and hemodynamic regressor at a finer time resolution
than the TRs, and later sample this regressor at the TR onset times.

This is what we do next.

```{python}
tr_divs = 100.0  # finer resolution has 100 steps per TR
```

With each TR divided into 100 intervals, one element corresponds to time
intervals of 1/100 of a TR:

```{python}
high_res_times = np.arange(0, n_trs, 1 / tr_divs) * TR
```

We will soon create a new neural prediction time-course where one element
corresponds to 1 / 100 of a TR:

```{python}
high_res_neural = np.zeros(high_res_times.shape)
```

We have the onset indices in terms of TRs, but now we want the onset indices
in terms of the new vector with 100 elements per TR:

```{python}
high_res_onset_indices = onsets_in_scans * tr_divs
high_res_onset_indices
```

In the same way, the durations were in seconds.  We divide by the TR to get
duration in terms of scans, then multiply by 100 to get the number in terms of
elements in the new `high_res_neural` time-course.

```{python}
high_res_durations = durations_seconds / TR * tr_divs
high_res_durations
```

Now we fill in the `high_res_neural` time course by setting values between
the start and the end of each event with the matching amplitudes:

```{python}
for hr_onset, hr_duration, amplitude in zip(
           high_res_onset_indices, high_res_durations, amplitudes):
    hr_onset = int(round(hr_onset))  # index - must be int
    hr_duration = int(round(hr_duration))  # makes index - must be int
    high_res_neural[hr_onset:hr_onset + hr_duration] = amplitude
plt.plot(high_res_times, high_res_neural)
plt.xlabel('Time (seconds)')
plt.ylabel('High resolution neural prediction')
```

We can use the hemodynamic response function we created earlier:

```{python}
from scipy.stats import gamma

def hrf(times):
    """ Return values for HRF at given times """
    # Gamma pdf for the peak
    peak_values = gamma.pdf(times, 6)
    # Gamma pdf for the undershoot
    undershoot_values = gamma.pdf(times, 12)
    # Combine them
    values = peak_values - 0.35 * undershoot_values
    # Scale max to 0.6
    return values / np.max(values) * 0.6
```

We are going to convolve at this higher time resolution.  First we need to
sample the HRF at this finer time resolution, to match the neural prediction:

```{python}
hrf_times = np.arange(0, 24, 1 / tr_divs)
hrf_at_hr = hrf(hrf_times)
```

Next we convolve the sampled HRF with the high resolution neural time course:

```{python}
high_res_hemo = np.convolve(high_res_neural, hrf_at_hr)
# Drop tail from convolution
high_res_hemo = high_res_hemo[:len(high_res_neural)]
plt.plot(high_res_times, high_res_hemo)
plt.xlabel('Time (seconds)')
plt.ylabel('High resolution convolved values')
len(high_res_times)
```

We can see that this is sampled at high resolution on the x axis by looking at
the first 20 TRs-worth of data:

```{python}
top_index = int(20 * tr_divs)
plt.plot(high_res_times[:top_index], high_res_hemo[:top_index], 'x:')
```

We can then subsample this high-resolution time course to get the values
corresponding to the start of each original TR (volume):

```{python}
tr_indices = np.arange(n_trs)
hr_tr_indices = np.round(tr_indices * tr_divs).astype(int)
tr_hemo = high_res_hemo[hr_tr_indices]
tr_times = tr_indices * TR  # times of TR onsets in seconds
plt.plot(tr_times, tr_hemo)
plt.xlabel('Time (seconds)')
plt.ylabel('Convolved values at TR onsets')
```

The first 20 TRs-worth of data shows these values are sampled every TR rather
than every 1/100th of a TR:

```{python}
plt.plot(tr_times[:20], tr_hemo[:20], 'x:')
plt.xlabel('Time (seconds)')
plt.ylabel('Convolved values at TR onsets')
```