%pylab inline

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'

import nipy.modalities.fmri.hrf as hrfs

t = np.arange(0, 25, 0.1)
hrf = hrfs.glovert(t)
plt.plot(t, hrf)
plt.xlabel('time (seconds)')

dt = 1/8. # time interval in seconds
t = np.arange(0, 40, dt) # times
experiment = np.zeros_like(t)
experiment[t == 5] = 1 # Event at 5 seconds
plt.plot(experiment)

hrf_t = np.arange(0, 24, dt)
hrf = hrfs.glovert(hrf_t)

hrf_experiment = np.convolve(experiment, hrf)
plt.plot(hrf_experiment)
hrf_experiment.shape

def dumb_convolve(data, kernel):
    data_len = len(data)
    kernel_len = len(kernel)
    out_len = data_len + kernel_len - 1
    out = np.zeros((out_len,))
    for i, value in enumerate(data):
        out[i:i+kernel_len] += value * kernel
    return out

plt.plot(dumb_convolve(experiment, hrf))

experiment[t == 15] = 1
plt.plot(dumb_convolve(experiment, hrf))

block = np.zeros((400,))
block[50:150] = 1
plt.plot(np.convolve(block, hrf))

dt = 1. # time interval in seconds
t_low = np.arange(0, 40, 1) # dt == 1
experiment_low = np.zeros_like(t_low)
experiment_low[5] = 1
experiment_low[15] = 1
hrf_low = hrfs.glovert(np.arange(24))
plt.plot(hrf_low)

data_len = len(experiment_low)
kernel_len = len(hrf_low)
out_len = data_len + kernel_len - 1
# Let's stack the summed time courses
components = np.zeros((data_len, out_len))
for i, value in enumerate(experiment_low):
   components[i, i:i+kernel_len] = value * hrf_low
# Then sum to get the result
out = np.sum(components, axis=0)

plt.plot(out)

plt.imshow(components)
plt.plot((0, data_len), (0, data_len), 'r')
plt.axis('tight')

def ugly_convolve(data, kernel):
    # Starts the same as the dumb convolve
    data_len = len(data)
    kernel_len = len(kernel)
    out_len = data_len + kernel_len - 1
    out = np.zeros((out_len,))
    # Different here though
    for i in range(out_len):
        # Sum up the value of the data * kernel going backwards through the data
        val = 0
        # Go back over the kernel and data
        for k in range(kernel_len):
            # Go back from current point in data
            data_ind = i - k
            # If we've run out of data, assume 0
            if data_ind < 0 or data_ind > data_len - 1:
                d_val = 0
            else:
                d_val = data[data_ind]
            # Multiply the corresponding kernel value with data
            val = val + d_val * kernel[k]
        out[i] = val
    return out

ugly_out = ugly_convolve(experiment_low, hrf_low)
plt.plot(ugly_out)
assert np.all(ugly_out == out)

up_hrfs_rows = kernel_len - 1 + data_len + kernel_len - 1
up_hrfs = np.zeros((up_hrfs_rows, out_len))
backwards_hrf = hrf_low[::-1]
for i in range(0, kernel_len -1 + data_len):
    up_hrfs[i:i+kernel_len, i] = backwards_hrf
plt.imshow(up_hrfs)

up_hrfs_trimmed = up_hrfs[kernel_len-1:kernel_len-1+data_len, :]
plt.imshow(up_hrfs_trimmed)
plt.plot((0, data_len), (0, data_len), 'r')
plt.axis('tight')
up_hrfs_trimmed.shape

fig, axes = plt.subplots(4, 1)
for i in range(4):
    axes[i].plot(up_hrfs_trimmed[:, 28+i])

data = np.tile(experiment_low[:, None], (1, out_len))
plt.imshow(data)

plt.imshow(data * up_hrfs_trimmed)
plt.plot((0, data_len), (0, data_len), 'r')
plt.axis('tight')
assert np.all(data * up_hrfs_trimmed == components)

upwards_out = np.sum(data * up_hrfs_trimmed, 0)
assert np.all(out == upwards_out)

filter_matrix = up_hrfs_trimmed.T
plt.imshow(filter_matrix)
filtered_out = filter_matrix.dot(experiment_low)
assert np.all(filtered_out == out)

plt.plot(filter_matrix[30, :])