% ICE12 - Demonstration of GLM in multiple voxels, assuming shape of HIRF

% Decide on "brain" and "anatomy"
brain = zeros([10 10]);
brain(2:4,2:4) = 1;
brain(6:8,6:8) = 2;
brain(5:9,[5 9]) = 3;
brain([5 9],6:8) = 3;
imagesc(brain); axis image

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ER design
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
waitUpFront = 5; % TRs so everything reaches steady-state and I have baseline
nStimTypes = 3;
ISIrange = [4:8]; % from 4 to 6 TRs between stim presentation
nStim = 10; % of each type
design = caoMakeER(waitUpFront,nStimTypes,ISIrange,nStim);

% Simulate different types of neural data
% Different amplitude of response to each stim type
neural = zeros(prod(size(brain)),length(design));
for b = 1:length(brain(:))
    if brain(b) == 1
        relResp = [1 .5 1.5]; % weaker response to stim2; stronger for stim3
		for n = 1:length(design)
            if design(n) % if there's a stimulus
                neural(b,n) = relResp(design(n)); % pick the appropriate response
            end
		end
    elseif brain(b) == 2
        relResp = [.5 1 1.5]; % weaker response to stim1; stronger for stim3
		for n = 1:length(design)
            if design(n) % if there's a stimulus
                neural(b,n) = relResp(design(n)); % pick the appropriate response
            end
		end
    elseif brain(b) == 3
        relResp = [1 1 1]; % same response to all stim
		for n = 1:length(design)
            if design(n) % if there's a stimulus
                neural(b,n) = relResp(design(n)); % pick the appropriate response
            end
		end
    end
end
neural = reshape(neural,[10 10 length(design)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic block design
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
blockLength = 10; % in units of TR; cycle = [blockLength(off) blockLength(on)];
nBlocksPerScan = 10; % my code adds a "on" half-cycle up front
design = caoMakeBlock(blockLength,nBlocksPerScan);

designMatrix = zeros(length(design),nStimTypes);
for n = 1:nStimTypes
    designMatrix(find(design == n),n) = 1;
end

% look at designMatrix and some sample voxels
clf; subplot(4,1,1); plot(designMatrix); axis tight
% a voxel of type 1
subplot(4,1,2); plot(squeeze(neural(2,2,:)),'c'); axis([0 185 0 2])
% type 2
subplot(4,1,3); plot(squeeze(neural(6,6,:)),'y'); axis([0 185 0 2])
% type 3
subplot(4,1,4); plot(squeeze(neural(5,5,:)),'r'); axis([0 185 0 2])


% Convolve design matrix with assumed HIRF
for n = 1:nStimTypes
    designHRF(:,n) = conv(designMatrix(:,n),hrf);
end
designHRF = designHRF(1:size(designMatrix,1),:);
designHRF = designHRF - repmat(mean(designHRF,1),[size(designHRF,1) 1]);
designHRF = designHRF./max(designHRF(:));

% Add some trends ...
drift = -1:2/(length(design)-1):1;
halfCos = cos(pi*(1:length(design))/length(design) - pi/2); 
  halfCos = halfCos - mean(halfCos);
  halfCos = halfCos/max(halfCos(:));
threeHalvesCos = cos(3*pi*(1:length(design))/length(design) - pi/2); 
  threeHalvesCos = threeHalvesCos - mean(threeHalvesCos(:));
  threeHalvesCos = threeHalvesCos./max(threeHalvesCos(:));
  
A = [designHRF drift' halfCos' threeHalvesCos'];
imagesc(A); axis off
colormap(gray); set(gcf,'Color',[1 1 1])

% Use the model to generate some (noiseless) data
for b = 1:prod(size(brain))
    if brain(b) == 1
        x(1:3,b) = [1 ; .5 ; 1.5];
    elseif brain(b) == 2
        x(1:3,b) = [.5 ; 1 ; 1.5]; 
    elseif brain(b) == 3
        x(1:3,b) = [1 ; 1 ; 1];
    end
    x(4:6,b) = rand([3 1]) - .5;
end

BOLD = A*x;

% now add some noise
noise = .2*randn(size(BOLD));

y = BOLD + noise;

%%  make a little tool that lets me click on a spot in the "brain" and see
%%  the underlying simulated BOLD
%% Could be done much more gracefully ...
done = 0;
figure(1); imagesc(brain); axis image
figure(2); clf
brainData = reshape(y',[10 10 size(A,1)]);
while ~done
    figure(1)
    [plot_y plot_x p] = impixel; % better right click!
    figure(2); plot(squeeze(brainData(plot_x,plot_y,:)))
    done = input('Done?  ');
end

%% Actually estimate regression weights
ATAinvAT = inv(A'*A)*A';
x_est = ATAinvAT*y;

% Repeat above exercise, showing fits ...

y_est = A*x_est;
done = 0;
figure(1); imagesc(brain); axis image
figure(2); clf
brainEst = reshape(y_est',[10 10 size(A,1)]);
while ~done
    figure(1)
    [plot_y plot_x p] = impixel; % better right click!
    figure(2); plot(squeeze(brainData(plot_x,plot_y,:)))
    hold on; plot(squeeze(brainEst(plot_x,plot_y,:)),'k'); hold off
    done = input('Done?  ');
end

%% Look at maps of responses to each stimulus


subplot(2,3,1); imagesc(reshape(x_est(1,:),[10 10]),[0 2]); axis image off
subplot(2,3,2); imagesc(reshape(x_est(2,:),[10 10]),[0 2]); axis image off
subplot(2,3,3); imagesc(reshape(x_est(3,:),[10 10]),[0 2]); axis image off

subplot(2,3,4); imagesc(reshape(x(1,:),[10 10]),[0 2]); axis image off
subplot(2,3,5); imagesc(reshape(x(2,:),[10 10]),[0 2]); axis image off
subplot(2,3,6); imagesc(reshape(x(3,:),[10 10]),[0 2]); axis image off