function ICE13b(exampleType)
% Usage:  ICE13b(exampleType)
%   Demonstration of GLM to estimate HRF in single voxel
%   for block design experiments
%    exampleType = 1:  slow ER; "same" HIRFs; no trends.
%    exampleType = 2:  slow ER; "same" HIRFs; trends.
%    exampleType = 3:  rapid ER; "same" HIRFs; no trends.
%    exampleType = 4:  rapid ER; diff HIRFs; no trends.
%    exampleType = 5:  rapid ER; diff HIRFs; trends.
% Uses functions:
%   caoMakeBlock.m
%   caoMakeStimConvMtx.m
%   caoComputeTTA.m
%   spm_hrf.m (which calls spm_Gpdf)

blockLength = 10;
nStimTypes = 1;
nBlocksPerScan = 10; % my code adds a "on" half-cycle up front
nStim = 10; % of each type
TR = 1.5;
nPtsEst = blockLength;

switch exampleType
    case 1; disp(' block exp''t; cannonical HIRF; no trends.')
        hrf = 1.5*spm_hrf(TR);
        x = zeros(nPtsEst + 3,1); x(11:13) = 0;
    case 2; disp(' block exp''t; cannonical HIRF; trends.')
        hrf = 1.5*spm_hrf(TR);
        x = zeros(nPtsEst + 3,1); x(11:13) = rand([3 1]) - .5;
    case 3; disp(' block exp''t; slow HIRF; no trends.')
        hrf = spm_hrf(TR,[4 16 1 1 3 0 32]);
        x = zeros(nPtsEst + 3,1); x(11:13) = rand([3 1]) - .5;
end

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

% Expand columns of designMatrix to estimate full HRF
designHRF = caoMakeStimConvMtx(designMatrix,nPtsEst);

% 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'];

figure(1); subplot(1,2,1)
imagesc(A); axis off; title('design mtx, full model')
colormap(gray); set(gcf,'Color',[1 1 1])

    figure(2); set(gcf,'Color',[1 1 1])
    plot(hrf,'b','linewidth',2); 
    set(gca,'FontSize',14,'FontWeight','bold')
    title('input HIRF'); xlabel('Time from stim. onset (TRs)'); ylabel('BOLD (a.u.)')
    x(1:10) = hrf(1:nPtsEst);


%  Use full model to generate some data with noise
BOLD = A*x;
noise = .2*randn(size(BOLD));
y = BOLD + noise;

figure(3); set(gcf,'Color',[1 1 1])
subplot(2,1,1)
plot(designMatrix); axis tight; 
set(gca,'FontSize',14,'FontWeight','bold'); title('design')
subplot(2,1,2)
plot(BOLD,'k'); hold on; plot(y,'r'); hold off; axis tight
set(gca,'FontSize',14,'FontWeight','bold'); title('simulated data')
xlabel('Time (TRs)'); ylabel('BOLD (a.u.)')

figure(4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "Analyze" with full model
ATAinvAT = inv(A'*A)*A';
x_est = ATAinvAT*y;
subplot(1,3,1)
plot(hrf(1:nPtsEst),'b:'); hold on
plot(x_est(1:10),'b','linewidth',2); hold off
 set(gca,'FontSize',12,'FontWeight','bold')
    axis tight;  title('full model'); set(gcf,'Color',[1 1 1])
 ylabel('BOLD (a.u.)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "Analyze" with TTA
subplot(1,3,2)
[TTA,TTAsem,count] = caoComputeTTA(design,y',nPtsEst);
plot(hrf(1:nPtsEst),'b:'); 
plot(TTA,'linewidth',2); hold off; axis tight; 
 set(gca,'FontSize',12,'FontWeight','bold')
 title('TTA'); xlabel('time from stim onset (TRs)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "Analyze" assuming shape of HIRF
subplot(1,3,3)
cannonicalHRF = spm_hrf(TR); clear A_conv
A_conv = conv(designMatrix(:,1),cannonicalHRF);
A_conv = [A_conv(1:size(designMatrix,1)) A(:,nPtsEst+(1:3))];

ATAinvAT = inv(A_conv'*A_conv)*A_conv';
x_est = ATAinvAT*y;
plot(hrf(1:nPtsEst),'b:'); hold on
plot(x_est(1)*cannonicalHRF(1:10),'b','linewidth',2); hold off
  set(gca,'FontSize',12,'FontWeight','bold')
    legend('input','est',4); axis tight; title('cannonical model')
 set(gcf,'Color',[1 1 1])
 
figure(1); subplot(1,2,2)
imagesc(A_conv); axis off; title('design mtx, assumed HIRF')
colormap(gray); set(gcf,'Color',[1 1 1])