# First set up some plotting things as usual
%pylab inline
import matplotlib as mpl
mpl.rcParams["figure.figsize"] = (9, 6)
mpl.rcParams["axes.grid"] = True
mpl.rcParams["lines.linewidth"] = 1.5

"""
Explanation: Linear Systems Tutorial
Class: Psych 204B
Tutorial: Linear Systems
Author: Wandell
Date: 3.31.04
Duration: 30 minutes
Copyright: Stanford University, Brian Wandell
Checked: 10/13/09 JW
Checked: 10/13/10 BW
Translated to Python by Michael Waskom, October 2012
End of explanation
"""


from scipy.special import gammaln
def spm_hrf(RT, fMRI_T=16):
    """Python implementation of spm_hrf"""
    _spm_Gpdf = lambda x, h, l: exp(h * log(l) + (h - 1) * log(x) - (l * x) - gammaln(h))
    dt = RT / float(fMRI_T)
    u = np.arange(1, int(32 / dt) + 1)
    hrf = np.concatenate(([0.0], _spm_Gpdf(u, 6 , dt) - _spm_Gpdf(u, 16, dt) / 6))
    idx = arange(int((32 / RT)) + 1) * fMRI_T
    hrf = hrf[idx]
    hrf = hrf / sum(hrf)
    return hrf

"""
Explanation: This tutorial introduces the student to the basic methods of linear
systems theory, focusing on the tools of shift-invariant linear systems.
The tutorial introduces the idea of superposition,  shift-invariance,
and impulse response functions. These ideas are fundamental to many
methods and analyses used in functional neuroimaging and through
engineering and science.
Linear systems theory is a method of characterizing certain types of
common systems.  A system is something that has an input and an output,
and thus we can think of it as a function output = L(input).
For example, the system might be an optical element, like a lens, that
takes an image as an input and produces another image as an output. Or,
 this system may be a biological system that we are measuring, for
instance the BOLD response that is the consequence of neural activity.
Characterizing the complete input-output properties of a system by 
exhaustive  measurement is usually impossible. When a system qualifies
as a linear system, it is possible to use the responses to a small set
of inputs to predict the response to any possible input. This can save
the scientist enormous amounts of work, and makes it possible to
characterize the system completely.
Not all systems can be described using linear systems theory. To use
linear methods to describe a system, it must satisfy the basic principle
of superposition.
The principle of superposition
The principle of superposition in words is simply this:
If the response to input A is Ra, and the response to input B is
Rb, then the response to the sum of A and B (A+B) is the sum of the
responses: Ra + Rb.
In equation form this is

If L(A) = Ra and L(B) = Rb, then L(A + B) = Ra + Rb

There is one simple consequence of superposition that is often taken as
a special case and named the principle of homogeneity.
Suppose that the inputs are equal, namely, A = B. 
It is easy to see that L(A + A) = Ra + Ra and thus L(2A) = 2Ra. 
In general, if s is any scalar, a linear system will satisfy the rule:

If L(A) = Ra, then L(sA) = sRa

An enormous number of simple systems satisfy superposition over a fairly
large range of inputs.  Linear systems theory is a basic tool used
throughout engineering, physics, and statistical modeling.  
In this tutorial, we will build up some of the tools of linear systems
theory.  These include the use of matrices to describe a linear system,
the very important special case of a shift-invariant linear system, and
the associated ideas of the impulse response function
All of these ideas are used in functional neuroimaging.  In fact, nearly
every field that uses linear systems theory develops its own terminology
(jargon) to describe the general ideas of linear systems theory in the
specific context of that field.   
In neuroimaging, the hemodynamic response function (HRF) is the impulse
response function. The analysis of the modulations in block-design
experiments are closely connected to the system transfer function.
Impulses and shift-invariance
time varying stimuli can be thought of as a series of impulses at
different levels. We can approximate any complex stimulus as if it were
simply the sum of a number of shifted and scaled impulses. 
A shift-invariant linear system is a special case when the system
responds in the same way to every impulse, no matter when the impulse
occurs in time.  That is, it doesn't matter what happened just before the
impulse or what happens just after the impulse.  The response to the
impulse will be the same.
The only difference in the response to different impulses is that the
responses are shifted in time to correspond to the impulse and multiplied
by a scalar that equals the size of the impulse.  
The significance of the principle of superposition is this:
If we know the impulse response function, we can predict the response to
any stimulus (combinations of impulses).  We compute the response to each
impulse and then add up all the results.
Practically, to characterize shift-invariant linear systems, then we need
to measure only one thing: the way the system responds to an impulse of a
particular intensity. This response is called the impulse response
function of the system. 
The analysis of BOLD data presupposes that the responses are
shift-invariant linear systems.  In this field, the impulse response
function is called the hemodynamic response function, or HRF.
In neuroimaging, the principles of  superposition and shift invariance
are essential to the ideas used in rapid event-related designs: The
underlying assumption is that if you have two neural (or cognitive)
events A and B occurring one after the other, then the BOLD response of A
and the BOLD response of B  will be identical (but B will be shifted in
time with respect to A) and thus total BOLD signal after these two events
is a sum of two independent responses shifted in time:
BOLD(A+B) = BOLD(A) + BOLD(B).
Now, we do some simple calculations. But first, we define a function to plot
a model hemodynamic response function.
End of explanation
"""


RT = 1  # Repetition time in seconds (sample every second)
hrf = spm_hrf(RT)

nT = len(hrf)   # Number of time steps
t = arange(nT)  # Individual time samples of the HRF

plot(t, hrf)
xlabel('Time (sec)')
ylabel('HRF level');

"""
Explanation: Let's take a look at this model (or "canonical" function)
End of explanation
"""


print(sum(hrf))

"""
Explanation: Notice that the hrf values sum to 1.
This means that convolution with the HRF will preserve the mean of the input stimulus level.
End of explanation
"""


def plot_imp_hrf(times):
    """Make a plot of impulse functions and estimated hemodynamic response."""
    
    # Get the predicted model
    hrf = spm_hrf(1)
    nT = len(hrf)
    t = arange(nT)
    
    # Make impulse vectors for each input stim
    stims = []
    for time in times:
        stim = zeros(nT)
        stim[time] = 1
        stims.append(stim)
    
    # Plot the stimulus
    figure(figsize=(9, 4))
    subplot(1, 2, 1)
    plot(zeros(nT))
    for time in times:
        plot([time, time], [0, .5])
    axis([0, 30, -0.5, 1])
    xlabel("Time (s)")
    title("Stimulus")
    
    # Convolve the stimulus with the hrf model
    sum_stim = sum(stims, axis=0)
    bold = convolve(sum_stim, hrf)
    
    # Plot the expected BOLD response
    subplot(1, 2, 2)
    plot(bold)
    axis([0, 30, -0.5, 1])
    xlabel("Time (s)")
    title("BOLD Response")

"""
Explanation: Also, notice that the default SPM HRF has a post-stimulus undershoot.
Next we will plot several instances of impulse functions with the corresponding estimated hemodynamic response using this model. To create the estimated response trace, we will use the convolve function from the numpy package (loaded automatically as we are using IPython in pylab mode; see above).
According to the principle of superposition the output should be the sum of the individual events. Thus, the output can be written as the sum of the outputs of stimuli 1 and 2.
Convolution is a mathematical function that applies the superposition principle: basically it sums over all impulse responses shifted in time
output = sum(hrf * input[t])

Since we'll be doing this a lot, let's write a function that takes a list of impulse times and makes the two plots.
End of explanation
"""


plot_imp_hrf([1])

"""
Explanation: Now let's see what this plot looks like for a simple single stimulus.
End of explanation
"""


plot_imp_hrf([1])
plot_imp_hrf([8])

"""
Explanation: We see that if we move around the time of the stimulus impulse, it also changes the timing of the expected BOLD response
End of explanation
"""


plot_imp_hrf([1, 8])

"""
Explanation: Now, plot the stimulus and output of a linear system that responds to both stimulus impulses.
End of explanation
"""


plot_imp_hrf([1, 5])

"""
Explanation: Block Design Experiments
We will start with blocks of 2 events and then continue to longer blocks.
We will examine how the number of events and the spacing between events
affects the predicted bold signal.
Example 1: Two events that are spaced 4 seconds apart
End of explanation
"""


plot_imp_hrf([1, 3])

"""
Explanation: Example 2: Two events that are spaced 2 seconds apart
Note that here you get one peak and not 2 peaks.  Why?
End of explanation
"""


plot_imp_hrf([1, 3, 5])

"""
Explanation: Example 3:  3 events that are spaced  2s apart
End of explanation
"""


plot_imp_hrf([1, 3, 5, 7, 9])

"""
Explanation: Example 4: 5 stimuli given 2 seconds -> block which is 10 seconds long
What changed as you increased the number of stimuli? Why?
End of explanation
"""


