%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\oxytocin\OXYT010\OXYT010_004_220126_112539',
            'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
            'F:\oxytocin\OXYT010\OXYT010_005_220126_150151',
            'F:\oxytocin\OXYT010\OXYT010_007_220203_092824',
            'F:\oxytocin\OXYT014\OXYT014_001_220429_155806',
            'F:\oxytocin\OXYT012\OXYT012_001_220323_130702',
            ]
dF_tranAll = [];dF_RNtranAll = [];dF_RWtranAll = [];dF_NWtranAll =[];
for ii = 1:6 %size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    
    info = read_Intan_RHD2000_file(basepath);
    analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
    lfp_fbtime = find(analogin(1,:)>0.6,1);lfp_fbtime(1,2) = find(analogin(1,:)>0.6,1,'last');
    delete analogin
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
        load([name,ext],[name,'_AChSti']);
        ACh_sti = eval([name,'_AChSti']);
        fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1,'last'));
        fbIDX = round(fbtime*20000);
        load([name,ext],[name,'_AChSig1']);
        ACh_signal = eval([name,'_AChSig1']);
        ACh_signal.data = downsample(ACh_signal.values((fbIDX(1,1)):(fbIDX(1,2))),200);
        ACh_signal.time = downsample(ACh_signal.times((fbIDX(1,1)):(fbIDX(1,2))),200)- fbtime(1,1);
        p = [];
        p = polyfit(ACh_signal.time,ACh_signal.data,5);
        Ybaseline = polyval(p,ACh_signal.time);
        figure
        plot(ACh_signal.time, ACh_signal.data); hold on;
        plot(ACh_signal.time,Ybaseline);
        dFoverF = smoothdata((ACh_signal.data - Ybaseline)./Ybaseline,'movmean',10);
        hold on;
        plot(ACh_signal.time, dFoverF);xlim([0 inf]);
    else load('dFoverF.mat')
    end
    
    index = find((SleepState.ints.REMstate(:,2)-SleepState.ints.REMstate(:,1))>20)
    SWS_REM = SleepState.ints.REMstate(index,1)-lfp_fbtime(1)/20000;
    dF_transit = [];
    for jj = 1:length(SWS_REM)
        idx = round((SWS_REM(jj) + [-20 20])*100);
        dF_transit(:,jj) = dFoverF(idx(1):idx(2));
    end
    dF_transit = dF_transit';
    dF_tranAll = [dF_tranAll;dF_transit];
    
   %% find REM-NREM and 
   RN_idx = []; RW_idx= [];
   for kk = 1:length(index)
     dx = find((SleepState.ints.NREMstate(:,1)-SleepState.ints.REMstate(index(kk),2)>0 & SleepState.ints.NREMstate(:,1)-SleepState.ints.REMstate(index(kk),2)<2 ))
     dxl = find((SleepState.ints.WAKEstate(:,1)-SleepState.ints.REMstate(index(kk),2)>0 & SleepState.ints.WAKEstate(:,1)-SleepState.ints.REMstate(index(kk),2)<2 ))
     if isempty(dx)==0 
         RN_idx = [RN_idx,index(kk)]
     end
     if isempty(dxl)==0
         RW_idx = [RW_idx,index(kk)]
     end
   end
   REM_RN = SleepState.ints.REMstate(RN_idx,2)-lfp_fbtime(1)/20000;
   dF_RNtransit = [];
   for jj = 1:length(REM_RN)
       idx = round((REM_RN(jj) + [-20 20])*100);
       dF_RNtransit(:,jj) = dFoverF(idx(1):idx(2));
   end
   dF_RNtransit = dF_RNtransit';
   dF_RNtranAll = [dF_RNtranAll;dF_RNtransit];
   
   REM_RW = SleepState.ints.REMstate(RW_idx,2)-lfp_fbtime(1)/20000;
   dF_RWtransit = [];
   for jj = 1:length(REM_RW)
       idx = round((REM_RW(jj) + [-20 20])*100);
       dF_RWtransit(:,jj) = dFoverF(idx(1):idx(2));
   end
   dF_RWtransit = dF_RWtransit';
   dF_RWtranAll = [dF_RWtranAll;dF_RWtransit];
   %% find WAKE-NREM
   NREM = SleepState.ints.NREMstate;
   index = 1:length(NREM);
   NW_idx = [];
   for kk = 1:length(index)
     dx = find((SleepState.ints.NREMstate(:,1)-SleepState.ints.WAKEstate(index(kk),2)>0 & SleepState.ints.NREMstate(:,1)-SleepState.ints.WAKEstate(index(kk),2)<2 ))
     if isempty(dx)==0 
         NW_idx = [NW_idx,index(kk)]
     end
   end
   NREM_WA = SleepState.ints.NREMstate(NW_idx,2)-lfp_fbtime(1)/20000;
   dF_NWtransit = [];
   for jj = 2:length(NREM_WA)-2
       idx = round((NREM_WA(jj) + [-50 50])*100);
       dF_NWtransit(:,jj) = dFoverF(idx(1):idx(2));
   end
   dF_NWtransit = dF_NWtransit';
   dF_NWtranAll = [dF_NWtranAll;dF_NWtransit];
end
figure;
outsr = 100;
imagesc([1:size(dF_tranAll,2)./outsr]'-20,[1:size(dF_tranAll,1)]',...
    zscore(dF_tranAll,[],2));
xlabel('Times from transit (s)'); ylabel('REM rumber #'); title('NREM-REM')

figure;
outsr = 100;
imagesc([1:size(dF_RNtranAll,2)./outsr]'-20,[1:size(dF_RNtranAll,1)]',...
    zscore(dF_RNtranAll,[],2));
xlabel('Times from transit (s)'); ylabel('REM rumber #');title('REM-NREM')

figure;
outsr = 100;
imagesc([1:size(dF_RWtranAll,2)./outsr]'-20,[1:size(dF_RWtranAll,1)]',...
    zscore(dF_RWtranAll,[],2));
xlabel('Times from transit (s)'); ylabel('REM rumber #');title('REM-WAKE')

figure;
outsr = 100;
imagesc([1:size(dF_NWtranAll,2)./outsr]'-50,[1:size(dF_NWtranAll,1)]',...
    zscore(dF_NWtranAll,[],2));
xlabel('Times from transit (s)'); ylabel('REM rumber #');title('NREM-WAKE')
