%Lecture_043012.m

%% Estimating 'beta weights'

% Standard single-event forward convolution:
maxt = 240;  %seconds
dt = 1;
t = [0:dt:(maxt-dt)]';
m = length(t);

% template hdr
maxth = 15;
th = [0:dt:(maxth-1)];
n = length(th);

nCascades = 2;
k = 1;
h = Gamma(nCascades,k,th)';

%stimulus
period = 40;
s=    mod(t,period)<= period/2;

figure(1)
clf
subplot(2,1,1)
stairs(t,s);
set(gca,'XLim',[0,maxt]);
set(gca,'YLim',[-.1,1.1]);
set(gca,'XTick',[0:period/2:maxt]);

%design matrix:
X = zeros(m,n);
temp = s;
for i=1:n
    X(:,i) = temp;
    temp = [0;temp(1:end-1)];
end

figure(2)
clf
imagesc(X);
axis equal
axis off
colormap(gray);

%% template fMRI response:
y = X*h;

figure(1)
subplot(2,1,2)
plot(t,y);
set(gca,'XLim',[0,maxt]);
set(gca,'XTick',[0:period/2:maxt]);

%% Add noise dc and trend to make 'fake' data:

k1= 3;  %scale parameter
k2= 100*rand(1);  %dc
k3=  10*randn(1);  %linear trend factor

resp = k1*y + .5*randn(size(t)) + ...
    k2+ ...
    k3*linspace(0,1,m)';

figure(3)
clf
plot(t,resp);
set(gca,'XLim',[0,maxt]);
set(gca,'XTick',[0:period/2:maxt]);


%% reconstruct the scale value, dc and trend.


%% The old way: reconstruct the hdr 
XX =  [X,ones(m,1),linspace(0,1,m)'];

hest = pinv(XX)*resp;

%% Estimating beta weights or how to estimate k1,k2 and k3:
Xh = X*h;  %unscaled guess at the fMRI response

XXh =  [Xh,ones(m,1),linspace(0,1,m)'];

kEst = pinv(XXh)*resp;

%kEst should look like [k1,k2,k3]'

%% Beta weights for multiple event types

nTypes = 4;
order = shuffle(1:4);

% blocked design stimulus
id = 1:15;
s = zeros(m,1);
for i = 1:3
    for j = 1:nTypes
        s(id) = order(j);
        id = id+20;
    end
end

figure(1)
clf
subplot(2,1,1)
stairs(t,s);
set(gca,'XLim',[0,maxt]);
set(gca,'XTick',[0:20:maxt]);

%design matrix
X = [];
for j=1:nTypes
    Xj = zeros(m,n);
    temp = s==j;

    for i=1:n
        Xj(:,i) = temp;
        temp = [0;temp(1:end-1)];
    end
    X = [X,Xj];
end

figure(2)
clf
imagesc(X);
colormap(gray);
axis equal
axis off


%% fake data again (amplitude, noise, dc trend)

amps = [1,2,3,4]; 

hh = h*amps;
hh = hh(:);

y = X*hh;
dc = 100*rand(1)
trend = 10*randn(1);

resp = y+.5*randn(size(t))+dc+trend*linspace(0,1,m)';

figure(1)
subplot(2,1,2)
plot(t,resp);
set(gca,'XLim',[0,maxt]);
set(gca,'XTick',[0:period/2:maxt]);


%% generate matrix H, an m x 4 matrix with columns containing h

H = zeros(n*nTypes,nTypes);

id = (1:n)';
for i=1:nTypes
    H(id,i) = h;
    id = id+n;
end

%predicted response for each event type

XX = X*H;

figure(1)
subplot(2,1,1);
hold on
plot(t,XX,'-');

%% estimate the amplitudes, dc and trend values

XXH = [XX,ones(m,1),linspace(0,1,m)'];

kEst = pinv(XXH)*resp;

% kEst(1:4) contains the four amplitudes

respHat = XXH*kEst;

subplot(2,1,2)
hold on
plot(t,respHat,'r-');

