%% make up a brain

brain =  zeros([10 10]);

brain(4:6,2:3) = 1;
brain(4:6,8:9) = 1;
brain(4:6,5:6) = 2;
imagesc(brain)

%% decide on experiment (this becomes design matrix later

% Stimulus 1 (flanker A)

TR = 1.5;
basicBlock = mod(round((4:163)/8),2);

%% prediction of neural activity

%% Basic model for activity during a block

tau_n = 2; % decay constant
t = 0:TR:TR*7; % time during 1 block
nAmp = 1; % relative strength of neural response
nSS = 0.5; % sustained neural response
neural = nSS + (nAmp - nSS)*exp(-t/tau_n);

%% String those together for full experiment
noFlankers = zeros([1 168]);
TR = 1.5;
t = 0:TR:TR*7; % time during 1 block
nAmp = 1; % relative strength of neural response
nSS = 0.5; % sustained neural response
blockLength = 16; % in units of TR
nBlocks = 10;
for n = 1:nBlocks
    noFlankers(n*blockLength:(n+1)*blockLength-9) = ...
        nSS + (nAmp - nSS)*exp(-t/tau_n);
end

%% String those together for full experiment
orthFlankers = zeros([1 168]);
TR = 1.5;
t = 0:TR:TR*7; % time during 1 block
nAmp = 1.5; % relative strength of neural response
nSS = 0.5; % sustained neural response
blockLength = 16; % in units of TR
nBlocks = 10;
for n = 1:nBlocks
    orthFlankers(n*blockLength:(n+1)*blockLength-9) = ...
        nSS + (nAmp - nSS)*exp(-t/tau_n);
end

%% String those together for full experiment
parFlankers = zeros([1 168]);
TR = 1.5;
t = 0:TR:TR*7; % time during 1 block
nAmp = 1; % relative strength of neural response
nSS = 0.33; % sustained neural response
blockLength = 16; % in units of TR
nBlocks = 10;
for n = 1:nBlocks
    parFlankers(n*blockLength:(n+1)*blockLength-9) = ...
        nSS + (nAmp - nSS)*exp(-t/tau_n);
end

t = 0:TR:TR*(length(noFlankers)-1);

subplot(3,1,1)

plot(t,noFlankers,'k','linewidth',2); hold on
plot(t,orthFlankers,'r','linewidth',2)
plot(t,parFlankers,'g','linewidth',2); hold off
set(gca,'FontSize',20)
legend('no flankers','orthog. flank','parallel flank')
title('predicted neural activity')

hirf = spm_hrf(1.5);

BOLD_no = conv(noFlankers,hirf);
BOLD_orth = conv(orthFlankers,hirf);
BOLD_par = conv(parFlankers,hirf);

subplot(3,1,2)
plot(t,BOLD_no(1:length(t)),'k','linewidth',2); hold on
plot(t,BOLD_orth(1:length(t)),'r','linewidth',2)
plot(t,BOLD_par(1:length(t)),'g','linewidth',2); hold off
set(gca,'FontSize',20)
legend('no flankers','orthog. flank','parallel flank')
title('predicted BOLD response')

CNR = 8;
measBOLD_no = BOLD_no(1:length(t)) + 1/CNR*randn(size(t));
measBOLD_orth = BOLD_orth(1:length(t)) + 1/CNR*randn(size(t));
measBOLD_par = BOLD_par(1:length(t)) + 1/CNR*randn(size(t));

subplot(3,1,3)
plot(t,measBOLD_no,'k','linewidth',2); hold on
plot(t,measBOLD_orth,'r','linewidth',2)
plot(t,measBOLD_par,'g','linewidth',2); hold off
legend('no flankers','orthog. flank','parallel flank')
set(gca,'FontSize',20)
title('predicted BOLD response')
set(gcf,'Color',[1 1 1])
xlabel('time (seconds)')

data1 = zeros([size(brain) length(measBOLD_no)]);
data2 = zeros([size(brain) length(measBOLD_no)]);
data3 = zeros([size(brain) length(measBOLD_no)]);
for r = 1:size(brain,1)
    for c = 1:size(brain,2)
            measBOLD_no = BOLD_no(1:length(t)) + 1/CNR*randn(size(t));
            measBOLD_orth = BOLD_orth(1:length(t)) + 1/CNR*randn(size(t));
            measBOLD_par = BOLD_par(1:length(t)) + 1/CNR*randn(size(t));
        if brain(r,c) == 2
            data1(r,c,:) = measBOLD_no;
            data2(r,c,:) = measBOLD_orth;
            data3(r,c,:) = measBOLD_par;
        elseif brain(r,c) == 1;
            data2(r,c,:) = measBOLD_orth;
            data3(r,c,:) = measBOLD_par;
        end
    end
end

subplot(2,1,1); imagesc(brain,[0 1]); axis image off
subplot(2,3,4); imagesc(data1(:,:,100),[0 1]); axis image off
subplot(2,3,5); imagesc(data2(:,:,100),[0 1]); axis image off
subplot(2,3,6); imagesc(data3(:,:,100),[0 1]); axis image off
colormap(gray); set(gcf,'Color',[1 1 1])
