clear; clc;
warning('off');

sesslist = {'LAMP3_012_190617_161856';
            'LAMP3_016_190619_180636';
            'LAMP3_018_190620_140743';
            'LAMP3_019_190620_153139';
            'LAMP3_020_190629_152350';
            'LAMP3_022_190630_151526';
            'LAMP3_007_190613_210018';
            'LAMP2_042_190614_120612'};
        
iplist = [1,7,8];
datapath = 'C:\Users\Liang\Dropbox\AchProject\Fig3 data\Non-stimulation\';

%% load data
plotdata = zeros(0); alphamap = false(81);
for ip = iplist
    session = sesslist{ip,1};
    load([datapath,session,filesep,'Achmapc.mat'],'AChmapc');
    AChmapc_z = (AChmapc-nanmean(AChmapc,'all'))./nanstd(AChmapc,[],'all');
    AChmapc_z(isnan(AChmapc_z)) = -3;
    AChmapc_z_sm = imgaussfilt(AChmapc_z,1,'FilterDomain','spatial');
    alphamap = alphamap | ~isnan(AChmapc);
    plotdata(:,:,ip) = AChmapc_z;
    scadata(:,:,ip) = (AChmapc-nanmean(AChmapc,'all'))./nanstd(AChmapc,[],'all');
    clear scatdata;
end

plotdata_m = mean(plotdata,3);
plotdata_m(alphamap) = (plotdata_m(alphamap)-nanmean(plotdata_m(alphamap)))./nanstd(plotdata_m(alphamap));

figure
im = imagesc(plotdata_m); set(gca,'YDir','normal');
xticks([]); yticks([]); clim([-1,1]*3); 
set(im,'AlphaData',alphamap);
cb = colorbar(); cb.Label.String = 'dFoverF (z-scored)'; 

%% postion bined
figure
im = imagesc(plotdata_m); set(gca,'YDir','normal');
xticks([]); yticks([]); clim([-1,1]*3); 
set(im,'AlphaData',alphamap);
cb = colorbar(); cb.Label.String = 'dFoverF (z-scored)'; 

delay = round([([-0.18,0.18]+0.02+0.4)/0.01,([-0.4,-0.1]+0.4)/0.01])+1;
pc = patch(([-0.18,0.18,0.18,-0.18]+0.02+0.4)/0.01,([-0.4,-0.4,-0.1,-0.1]+0.4)/0.01,'r');
pc.FaceAlpha = 0.1; pc.EdgeColor = 'none';

center = round([([-0.1,0.1]+0.02+0.4)/0.01,([-0.1,0.4]+0.4)/0.01])+1;
pc = patch(([-0.1,0.1,0.1,-0.1]+0.02+0.4)/0.01,([-0.1,-0.1,0.4,0.4]+0.4)/0.01,'b');
pc.FaceAlpha = 0.1; pc.EdgeColor = 'none';

achdelay = squeeze(nanmean(scadata(delay(1):delay(2),delay(3):delay(4),:),[1,2]));
achcenter = squeeze(nanmean(scadata(center(1):center(2),center(3):center(4),:),[1,2]));

left1 = round([([-0.4,-0.18]+0.02+0.4)/0.01,([-0.4,-0.1]+0.4)/0.01])+1;
pc = patch(([-0.4,-0.18,-0.18,-0.4]+0.02+0.4)/0.01,([-0.4,-0.4,-0.1,-0.1]+0.4)/0.01,'c');
pc.FaceAlpha = 0.2; pc.EdgeColor = 'none';
left2 = round([([-0.4,-0.1]+0.02+0.4)/0.01,([-0.1,0.4]+0.4)/0.01])+1;
pc = patch(([-0.4,-0.1,-0.1,-0.4]+0.02+0.4)/0.01,([-0.1,-0.1,0.4,0.4]+0.4)/0.01,'c');
pc.FaceAlpha = 0.2; pc.EdgeColor = 'none';

right1 = round([([0.18,0.38]+0.02+0.4)/0.01,([-0.4,-0.1]+0.4)/0.01])+1;
pc = patch(([0.18,0.38,0.38,0.18]+0.02+0.4)/0.01,([-0.4,-0.4,-0.1,-0.1]+0.4)/0.01,'y');
pc.FaceAlpha = 0.2; pc.EdgeColor = 'none';
right2 = round([([0.1,0.38]+0.02+0.4)/0.01,([-0.1,0.4]+0.4)/0.01])+1;
pc = patch(([0.1,0.38,0.38,0.1]+0.02+0.4)/0.01,([-0.1,-0.1,0.4,0.4]+0.4)/0.01,'y');
pc.FaceAlpha = 0.2; pc.EdgeColor = 'none';

rightsz1 = numel(scadata(right1(1):right1(2),right1(3):right1(4),1));
rightach1 = squeeze(nanmean(scadata(right1(1):right1(2),right1(3):right1(4),:),[1,2]));
rightsz2 = numel(scadata(right2(1):right2(2),right2(3):right2(4),1));
rightach2 = squeeze(nanmean(scadata(right2(1):right2(2),right2(3):right2(4),:),[1,2]));
achright = (rightach1*rightsz1 + rightach2*rightsz2)./(rightsz1+rightsz2);

leftsz1 = numel(scadata(left1(1):left1(2),left1(3):left1(4),1));
leftach1 = squeeze(nanmean(scadata(left1(1):left1(2),left1(3):left1(4),:),[1,2]));
leftsz2 = numel(scadata(left2(1):left2(2),left2(3):left2(4),1));
leftach2 = squeeze(nanmean(scadata(left2(1):left2(2),left2(3):left2(4),:),[1,2]));
achleft = (leftach1*leftsz1 + leftach2*leftsz2)./(leftsz1+leftsz2);

non_data = [achdelay,achcenter,achleft,achright];
save('ach_non','non_data');