function ICE13(exampleType)
% Usage:  ICE13(exampleType)
%   Demonstration of GLM to estimate HRF in single voxel
%    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:
%   caoMakeER.m
%   caoMakeStimConvMtx.m
%   caoComputeTTA.m
%   spm_hrf.m (which calls spm_Gpdf)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ER design
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters for all examples
waitUpFront = 5; % TRs so everything reaches steady-state and I have baseline
nStimTypes = 3; % this is fixed for now because of the HIRF models
nStim = 50; % of each type
TR = 1.5;
nPtsEst = 15;

switch exampleType
    case 1; disp(' slow ER; same HIRFs; no trends.')
        ISIrange = [10]; % from 4 to 6 TRs between stim presentation
        hrf1 = 1.5*spm_hrf(TR);
        hrf2 = spm_hrf(TR);
        hrf3 = 0.8*spm_hrf(TR);
        x = zeros(nStimTypes*nPtsEst + 3,1); x(31:33) = 0;
    case 2; disp(' slow ER; same HIRFs; trends.')
        ISIrange = [10]; % from 4 to 6 TRs between stim presentation
        hrf1 = 1.5*spm_hrf(TR);
        hrf2 = spm_hrf(TR);
        hrf3 = 0.8*spm_hrf(TR);
        x = zeros(nStimTypes*nPtsEst + 3,1); x(31:33) = rand([3 1]) - 0.5;
    case 3; disp(' rapid ER; same HIRFs; no trends.')
        ISIrange = [4:6]; % from 4 to 6 TRs between stim presentation
        hrf1 = 1.5*spm_hrf(TR);
        hrf2 = spm_hrf(TR);
        hrf3 = 0.8*spm_hrf(TR);;
        x = zeros(nStimTypes*nPtsEst + 3,1); x(31:33) = 0;
    case 4; disp(' rapid ER; diff HIRFs; no trends.')
        ISIrange = [4:6]; % from 4 to 6 TRs between stim presentation
        hrf1 = spm_hrf(TR,[6 16 1 1 6 0 32]);
        hrf2 = spm_hrf(TR,[4 16 1 1 3 0 32]);
        hrf3 = spm_hrf(TR,[6 24 1 1 10 0 32]);
        x = zeros(nStimTypes*nPtsEst + 3,1); x(31:33) = 0;
    case 5; disp(' rapid ER; diff HIRFs; trends.')
        ISIrange = [4:6]; % from 4 to 6 TRs between stim presentation
        hrf1 = spm_hrf(TR,[6 16 1 1 6 0 32]);
        hrf2 = spm_hrf(TR,[4 16 1 1 3 0 32]);
        hrf3 = spm_hrf(TR,[6 24 1 1 10 0 32]);
        x = zeros(nStimTypes*nPtsEst + 3,1); x(31:33) = rand([3 1]) - 0.5;
end


design = caoMakeER(waitUpFront,nStimTypes,ISIrange,nStim);
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(hrf1,'b','linewidth',2); hold on
plot(hrf2,'g','linewidth',2); plot(hrf3,'r','linewidth',2); hold off
set(gca,'FontSize',14,'FontWeight','bold')
legend('response to stim1','response to stim2','response to stim3')
xlabel('Time from stim. onset (TRs)'); ylabel('BOLD (a.u.)')

% Use the model to generate some (noiseless) data
x(1:nPtsEst) = hrf1(1:nPtsEst);
x(nPtsEst+(1:nPtsEst)) = hrf2(1:nPtsEst);
x(2*nPtsEst+(1:nPtsEst)) = hrf3(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(hrf1(1:nPtsEst),'b:'); hold on
plot(hrf2(1:nPtsEst),'g:'); plot(hrf3(1:nPtsEst),'r:');
plot(x_est(1:nPtsEst),'b','linewidth',2); 
plot(x_est(nPtsEst + (1:nPtsEst)),'g','linewidth',2);
plot(x_est(2*nPtsEst + (1:nPtsEst)),'r','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(hrf1(1:nPtsEst),'b:'); hold on
plot(hrf2(1:nPtsEst),'g:'); plot(hrf3(1:nPtsEst),'r:');
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
for n = 1:3
    A_conv(:,n) = conv(designMatrix(:,n),cannonicalHRF);
end
A_conv = [A_conv(1:size(designMatrix,1),:) A(:,nStimTypes*nPtsEst + (1:3))];

ATAinvAT = inv(A_conv'*A_conv)*A_conv';
x_est = ATAinvAT*y;
plot(hrf1(1:nPtsEst),'b:'); hold on
plot(hrf2(1:nPtsEst),'g:'); plot(hrf3(1:nPtsEst),'r:');
plot(x_est(1)*cannonicalHRF(1:nPtsEst),'b','linewidth',2); 
plot(x_est(2)*cannonicalHRF(1:nPtsEst),'g','linewidth',2);
plot(x_est(3)*cannonicalHRF(1:nPtsEst),'r','linewidth',2); hold off
  set(gca,'FontSize',12,'FontWeight','bold')
    legend('s1','s2','s3',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])