파일 업로드 'Matlab_function'

Matlab_function/GetAvailablePCAUnit.m

@@ -0,0 +1,19 @@
+% function Data=GetAvailablePCAUnit(Data1, Data2)
+% FR_Thre=0.25; 
+% k=1; Data=zeros(size(Data1,1),1);
+% for clstindex=1:size(Data1,2)
+%     if mean(Data2(:,clstindex)) > FR_Thre% && ~(length(find(Data2(:,clstindex)==0)) == 30)
+%         Data(:,k)=Data1(:,clstindex); k=k+1;
+%        
+%     end
+% end
+
+function [Data, Data_Name] = GetAvailablePCAUnit(Data2, FR_Thre,Name)
+% FR_Thre=0.5; 
+k=1; Data=zeros(size(Data2,1),1);
+for clstindex=1:size(Data2,2)
+    if nanmean(Data2(:,clstindex)) > FR_Thre && (length(find(Data2(:,clstindex)==0)) <= size(Data2,1)*0.5)
+        Data(:,k)=Data2(:,clstindex); 
+        Data_Name(k,1)=Name(clstindex); k=k+1;
+    end
+end

Matlab_function/GetExcludeNan.m

@@ -0,0 +1,8 @@
+function data = GetExcludeNan(y)
+k=1;
+for i=1:length(y)
+    if ~isnan(y(i))
+        data(k)=y(i);
+        k=k+1;
+    end
+end

Matlab_function/GetFiringRate.m

@@ -0,0 +1,3 @@
+function AvgFiringRate = GetFiringRate(SkaggsRateMap)
+AvgFiringRate=nansum(nansum(SkaggsRateMap))/length(find(~isnan(SkaggsRateMap)==1));
+

Matlab_function/GetGaussianSmoothing.m

@@ -0,0 +1,23 @@
+function Data_Smoothing = GetGaussianSmoothing(Data,Timewindow)
+%% Smoothing processing
+Data_Smoothing=zeros(1,size(Data,2));
+gaussFilter = gausswin(Timewindow,2.5);
+gaussFilter=gaussFilter';
+gaussFilter = gaussFilter / sum(gaussFilter); % Normalize.
+
+% % Data transformation for smoothing
+if size(Data,2)==1 && size(Data,1) > 1
+    Data=Data';
+end
+
+if size(Data,1) == 1 %% 1 * n matrix
+    Data_Smoothing=conv(Data, gaussFilter);    
+    %     Eliminate first and last bins contained unsufficient information.
+    %     that created by smoothing process
+    Data_Smoothing=Data_Smoothing((floor(Timewindow/2)+1):end-floor(Timewindow/2));
+elseif size(Data,1) > 1 %% m * n matrix
+    for i=1:size(Data,1)
+        temp(i,:)=conv(Data(i,:), gaussFilter);    
+        Data_Smoothing(i,:)=temp(i,(floor(Timewindow/2)+1):end-floor(Timewindow/2));
+    end
+end

Matlab_function/GetGroupingVar.m

@@ -0,0 +1,15 @@
+function Data = GetGroupingVar(data,varargin)
+if nargin == 1
+    for i=1:length(data)
+        Data(i,1)=i;
+    end
+elseif nargin == 2
+    if ischar(varargin{1})
+        for i=1:length(data)
+            Data{i,1} = varargin{1};
+        end
+    else
+        n=varargin{1};
+        Data = n*ones(length(data),1);
+    end
+end

Matlab_function/GetPC_GaussianFiltering.m

@@ -0,0 +1,15 @@
+function z = GetPC_GaussianFiltering(x, window, gaussFilter)
+
+for i=1:(window-1)/2
+    temp = x((1:i+(window-1)/2),:);
+    z(i,:) = (temp'*gaussFilter((window+1)/2+1-i:window,i))';
+end
+for i=(window+1)/2:size(x,1)-(window-1)/2
+    temp = x((i-(window-1)/2:i+(window-(window+1)/2)),:);
+    z(i,:) = (temp'*gaussFilter(:,(window+1)/2))';
+end
+j=1;
+for i=size(x,1)-(window-1)/2+1:size(x,1)
+    temp=x((end-window+1+j:end),:);
+    z(i,:) = (temp'*gaussFilter(1:window-j,(window+1)/2-j))'; j=j+1;
+end

Matlab_function/GetPopulationCorrelation.m

@@ -0,0 +1,10 @@
+function [Data1, Data2] = GetPopulationCorrelation(data1, data2)
+k=1;
+for i=1:size(data1,1)
+    if ~(sum(double(isnan(data1(i,:)))) || sum(double(isnan(data2(i,:)))))
+        Data1(k,:)=data1(i,:);
+        Data2(k,:)=data2(i,:);
+        k=k+1;
+    end
+end
+

Matlab_function/GetSpaInfo.m

@@ -0,0 +1,29 @@
+function SpatialInformation = GetSpaInfo(OCCMat, SkaggsrateMat)
+%Calculate spatial information [Skaggs']
+%SpatialInformation = GetSpaInfo(OCCMat, SkaggsrateMat)
+%
+%Related function(s)
+%jkABM.m
+%outputs from jkABM.m are needed to calculate spatial information scores.
+%
+%Input format
+% OCCMat	[nROW x nCOL matrix]									raw occupancy amp
+% SkaggsrateMat [nROW x nCOL matrix]								Skaggs' ABM firing rate map
+%
+%Output format
+% SpatialInformation [1 x 1]										bit/spike spatial information scores
+%
+%Originally from VB codes which Inah Lee has.
+%Translate VB codes to matlab was done by [Jangjin Kim, July-13-2008]
+%Verification was done
+
+criteria = .0001;
+
+SumRate = nansum(nansum(SkaggsrateMat));
+SumOcc = sum(sum(~isnan(OCCMat)));
+MeanRate = SumRate / SumOcc;
+
+H1 = sum(SkaggsrateMat(SkaggsrateMat > criteria) * -1 .* log2(SkaggsrateMat(SkaggsrateMat > criteria)) / SumOcc);
+H0 = MeanRate * -1 * log2(MeanRate);
+
+SpatialInformation = (H0 - H1) / MeanRate;

Matlab_function/GetSpatialCorrelation.m

@@ -0,0 +1,56 @@
+function SpaCorr12 = GetSpatialCorrelation(SkaggsRatemap1, SkaggasRatemap2)
+
+SkaggsRatemap1_1D=0;
+SkaggasRatemap2_1D=0;
+
+for i=1:size(SkaggsRatemap1,1)
+    for j=1:size(SkaggsRatemap1,2)
+        SkaggsRatemap1_1D(end+1)=SkaggsRatemap1(i,j);
+        SkaggasRatemap2_1D(end+1)=SkaggasRatemap2(i,j);
+    end
+end
+SkaggsRatemap1_1D=SkaggsRatemap1_1D(2:end);
+SkaggasRatemap2_1D=SkaggasRatemap2_1D(2:end);
+
+SkaggsRatemap1_1D_ne=0;
+SkaggasRatemap2_1D_ne=0;
+
+for i=1:length(SkaggsRatemap1_1D)
+    if ~isnan(SkaggsRatemap1_1D(i)) && ~isnan(SkaggasRatemap2_1D(i))
+        SkaggsRatemap1_1D_ne(end+1)=SkaggsRatemap1_1D(i);
+        SkaggasRatemap2_1D_ne(end+1)=SkaggasRatemap2_1D(i);
+    end
+end
+SkaggsRatemap1_1D_ne=SkaggsRatemap1_1D_ne(2:end);
+SkaggasRatemap2_1D_ne=SkaggasRatemap2_1D_ne(2:end);
+SpaCorr12=corr2(SkaggsRatemap1_1D_ne,SkaggasRatemap2_1D_ne);
+
+
+
+% SkaggsRatemap1_1D=0;
+% SkaggasRatemap2_1D=0;
+% 
+% for i=1:48
+%     for j=1:72
+%         SkaggsRatemap1_1D(end+1)=SkaggsRatemap1(i,j);
+%         SkaggasRatemap2_1D(end+1)=SkaggasRatemap2(i,j);
+%     end
+% end
+% SkaggsRatemap1_1D=SkaggsRatemap1_1D(2:end);
+% SkaggasRatemap2_1D=SkaggasRatemap2_1D(2:end);
+% 
+% SkaggsRatemap1_1D_ne=0;
+% SkaggasRatemap2_1D_ne=0;
+% 
+% for i=1:length(SkaggsRatemap1_1D)
+%     if ~isnan(SkaggsRatemap1_1D(i)) && ~isnan(SkaggasRatemap2_1D(i))
+%         SkaggsRatemap1_1D_ne(end+1)=SkaggsRatemap1_1D(i);
+%         SkaggasRatemap2_1D_ne(end+1)=SkaggasRatemap2_1D(i);
+%     end
+% end
+% SkaggsRatemap1_1D_ne=SkaggsRatemap1_1D_ne(2:end);
+% SkaggasRatemap2_1D_ne=SkaggasRatemap2_1D_ne(2:end);
+% SpaCorr12=corr2(SkaggsRatemap1_1D_ne,SkaggasRatemap2_1D_ne);
+% 
+
+

Matlab_function/GetSpatialInfo_pvalue_ALT.m

@@ -0,0 +1,85 @@
+function Spainfo_pvalue = GetSpatialInfo_pvalue_ALT(thisCLST, Pos, SpaInfo)
+randN = 100;%
+infoArray = nan(randN, 1);
+randArray = nan(randN, 1);
+SpaInfo_Random_Array=0;
+
+XsizeOfVideo = 720;
+YsizeOfVideo = 480;
+samplingRate = 30;
+scaleForRateMap = 10;
+binXForRateMap = XsizeOfVideo / scaleForRateMap;
+binYForRateMap = YsizeOfVideo / scaleForRateMap;
+% parfor randI = 1:randN
+spkN = size(thisCLST.PositionLocked(thisCLST.Flag_Position_pvalue), 1);
+if spkN ~= 0
+    for randI = 1:randN
+        %input x, y, t: position data. it must include only behavior session data.
+        %t_spk: spk data time stamps. it must include only behavior session data.
+        %SpaInfo : Real spatial information score
+        
+        x_spkR = nan(spkN, 1);
+        y_spkR = nan(spkN, 1);
+        t_spkR = nan(spkN, 1);
+        
+        occN = size(Pos.t(Pos.Flag_Position_pvalue),1);
+        T=Pos.t(Pos.Flag_Position_pvalue);
+        startTime =T(1);
+        endTime = T(end);
+        Endtime=T(end);
+        if occN < 1500 %minimum value
+            error('occ data must be greater than or equare to 1500');
+        end
+        
+        for i=1:spkN
+            while true
+                cRand = rand();
+                shift = occN * cRand; %shift amount (unit: occ)
+                
+                if shift > 300 && (occN-shift) > 300 % sufficiantly far from original position
+                    Shift(i)=shift; % randomize the every spike's positions.
+                    break;
+                end
+            end
+        end
+        t_spk=thisCLST.PositionLocked(thisCLST.Flag_Position_pvalue);
+        x_spk = thisCLST.x(thisCLST.Flag_Position_pvalue);
+        y_spk = thisCLST.y(thisCLST.Flag_Position_pvalue);
+        t_pos=Pos.t(Pos.Flag_Position_pvalue);
+        x_pos=Pos.x(Pos.Flag_Position_pvalue);
+        y_pos=Pos.y(Pos.Flag_Position_pvalue);
+        
+        spk_index=zeros(1,length(t_spk));
+        for i=1:length(t_pos)
+            for j=1:length(t_spk)
+                if(t_spk(j)==t_pos(i))
+                    spk_index(j)=i;
+                end
+            end
+        end
+        
+        spkR_index=spk_index+floor(Shift);
+        outRange = spkR_index > occN;
+        spkR_index(outRange) = spkR_index(outRange) - occN;
+        
+        for i=1:length(x_spk)
+            x_spkR(i) = x_pos(spkR_index(i));
+            y_spkR(i) = y_pos(spkR_index(i));
+            t_spkR(i) = t_pos(spkR_index(i));
+        end
+        
+        %%
+        [occMatR, spkMatR, rawMatR, skaggsrateMatR] = abmFiringRateMap( ...
+            [t_spkR, x_spkR, y_spkR], ...
+            [Pos.t(Pos.Flag_Position_pvalue), Pos.x(Pos.Flag_Position_pvalue), Pos.y(Pos.Flag_Position_pvalue)], ...
+            binYForRateMap, binXForRateMap, scaleForRateMap, samplingRate);
+        
+        SpaInfo_Random= GetSpaInfo(occMatR, skaggsrateMatR);
+        SpaInfo_Random_Array(randI) = SpaInfo_Random;
+        rand_Array(randI) = cRand;
+    end
+    
+    Spainfo_pvalue = sum(SpaInfo_Random_Array > SpaInfo) / randN;
+else
+    Spainfo_pvalue = NaN;
+end

Matlab_function/GetSpikeOutboundFlag.m

@@ -0,0 +1,177 @@
+function [thisCLST, Time, TrialNumber] = GetSpikeOutboundFlag(thisCLST, TrialNumber, Time)
+rowSpike=length(thisCLST.t);
+thisCLST.Flag_OutboundSpk_Track=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk=zeros(rowSpike,1);
+Time.latency=zeros(1,TrialNumber.Trial);
+for j=1:TrialNumber.Trial
+    if j==1
+        Time.latency(j)=Time.Choice(j,1)-Time.Entry(j,1);
+        for i=1:rowSpike
+            if (thisCLST.t(i)>=(Time.Entry(j,1)))&&(thisCLST.t(i)<Time.Choice(j,1))&&(thisCLST.SmoothingVelocity(i) > 5)
+                thisCLST.Flag_OutboundSpk_Track(i)=j;  
+            end
+        end
+
+    elseif (j~=1) %
+        for i=1:rowSpike
+            if (thisCLST.t(i)>=(Time.Choice(j-1,1)+Time.LeaveRewardZone(j-1,1)))&&(thisCLST.t(i)<Time.Choice(j,1))&&(thisCLST.SmoothingVelocity(i) > 5)
+                thisCLST.Flag_OutboundSpk_Track(i)=j;                    
+            end
+        end
+    end
+end
+%
+for j=1:TrialNumber.Trial
+    if j==1
+        for i=1:rowSpike
+            if (thisCLST.t(i)>=Time.Entry(j,1))&&(thisCLST.t(i)<Time.Choice(j,1))&&(thisCLST.SmoothingVelocity(i) > 5)
+                thisCLST.Flag_OutboundSpk(i)=j;  
+            end
+        end
+
+    elseif (j~=1) %
+        for i=1:rowSpike
+            if (thisCLST.t(i)>=Time.Choice(j-1,1))&&(thisCLST.t(i)<Time.Choice(j,1))&&(thisCLST.SmoothingVelocity(i) > 5)
+                thisCLST.Flag_OutboundSpk(i)=j;                    
+            end
+        end
+    end
+end
+
+%% Individual
+thisCLST.Flag_OutboundSpk_Individual_Track=cell(1,TrialNumber.Trial);
+thisCLST.Flag_OutboundSpk_Individual=cell(1,TrialNumber.Trial);
+for j=1:TrialNumber.Trial
+    thisCLST.Flag_OutboundSpk_Individual_Track{1,j}=logical(thisCLST.Flag_OutboundSpk_Track==j);
+    thisCLST.Flag_OutboundSpk_Individual{1,j}=logical(thisCLST.Flag_OutboundSpk==j);
+end
+
+%% Total
+thisCLST.Flag_OutboundSpk_Odd_Total_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Odd_Total{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Odd_Total{1,1}=zeros(row,1);
+for i=1:2:TrialNumber.Trial
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Odd_Total_Track{1,1}=thisCLST.Flag_OutboundSpk_Odd_Total_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Odd_Total{1,1}=thisCLST.Flag_OutboundSpk_Odd_Total{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Odd_Total{1,1}=thisCLST.Flag_OutboundSpk_Odd_Total{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+thisCLST.Flag_OutboundSpk_Even_Total_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Even_Total{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Even_Total{1,1}=zeros(row,1);
+for i=2:2:TrialNumber.Trial
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Even_Total_Track{1,1}=thisCLST.Flag_OutboundSpk_Even_Total_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Even_Total{1,1}=thisCLST.Flag_OutboundSpk_Even_Total{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Even_Total{1,1}=thisCLST.Flag_OutboundSpk_Even_Total{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+
+%% B1
+thisCLST.Flag_OutboundSpk_Odd_B1_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Odd_B1{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Odd_B1{1,1}=zeros(row,1);
+for i=TrialNumber.Block(1)+1:2:TrialNumber.Block(2)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Odd_B1_Track{1,1}=thisCLST.Flag_OutboundSpk_Odd_B1_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Odd_B1{1,1}=thisCLST.Flag_OutboundSpk_Odd_B1{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Odd_B1{1,1}=thisCLST.Flag_OutboundSpk_Odd_B1{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+thisCLST.Flag_OutboundSpk_Even_B1_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Even_B1{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Even_B1{1,1}=zeros(row,1);
+for i=TrialNumber.Block(1)+2:2:TrialNumber.Block(2)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Even_B1_Track{1,1}=thisCLST.Flag_OutboundSpk_Even_B1_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Even_B1{1,1}=thisCLST.Flag_OutboundSpk_Even_B1{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Even_B1{1,1}=thisCLST.Flag_OutboundSpk_Even_B1{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+
+%% B2
+thisCLST.Flag_OutboundSpk_Odd_B2_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Odd_B2{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Odd_B2{1,1}=zeros(row,1);
+for i=TrialNumber.Block(2)+1:2:TrialNumber.Block(3)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Odd_B2_Track{1,1}=thisCLST.Flag_OutboundSpk_Odd_B2_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Odd_B2{1,1}=thisCLST.Flag_OutboundSpk_Odd_B2{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Odd_B2{1,1}=thisCLST.Flag_OutboundSpk_Odd_B2{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+thisCLST.Flag_OutboundSpk_Even_B2_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Even_B2{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Even_B2{1,1}=zeros(row,1);
+for i=TrialNumber.Block(2)+2:2:TrialNumber.Block(3)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Even_B2_Track{1,1}=thisCLST.Flag_OutboundSpk_Even_B2_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Even_B2{1,1}=thisCLST.Flag_OutboundSpk_Even_B2{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Even_B2{1,1}=thisCLST.Flag_OutboundSpk_Even_B2{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+
+%% B3
+thisCLST.Flag_OutboundSpk_Odd_B3_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Odd_B3{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Odd_B3{1,1}=zeros(row,1);
+for i=TrialNumber.Block(3)+1:2:TrialNumber.Block(4)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Odd_B3_Track{1,1}=thisCLST.Flag_OutboundSpk_Odd_B3_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Odd_B3{1,1}=thisCLST.Flag_OutboundSpk_Odd_B3{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Odd_B3{1,1}=thisCLST.Flag_OutboundSpk_Odd_B3{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+thisCLST.Flag_OutboundSpk_Even_B3_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Even_B3{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Even_B3{1,1}=zeros(row,1);
+for i=TrialNumber.Block(3)+2:2:TrialNumber.Block(4)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Even_B3_Track{1,1}=thisCLST.Flag_OutboundSpk_Even_B3_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Even_B3{1,1}=thisCLST.Flag_OutboundSpk_Even_B3{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Even_B3{1,1}=thisCLST.Flag_OutboundSpk_Even_B3{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+
+%% B4
+thisCLST.Flag_OutboundSpk_Odd_B4_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Odd_B4{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Odd_B4{1,1}=zeros(row,1);
+for i=TrialNumber.Block(4)+1:2:TrialNumber.Block(5)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Odd_B4_Track{1,1}=thisCLST.Flag_OutboundSpk_Odd_B4_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Odd_B4{1,1}=thisCLST.Flag_OutboundSpk_Odd_B4{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Odd_B4{1,1}=thisCLST.Flag_OutboundSpk_Odd_B4{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end
+thisCLST.Flag_OutboundSpk_Even_B4_Track{1,1}=zeros(rowSpike,1);
+thisCLST.Flag_OutboundSpk_Even_B4{1,1}=zeros(rowSpike,1);
+% thisCLST.Flag_VM_OutboundSpk_Even_B4{1,1}=zeros(row,1);
+for i=TrialNumber.Block(4)+2:2:TrialNumber.Block(5)
+    if ~sum(i==TrialNumber.VOID)
+        thisCLST.Flag_OutboundSpk_Even_B4_Track{1,1}=thisCLST.Flag_OutboundSpk_Even_B4_Track{1,1}|thisCLST.Flag_OutboundSpk_Individual_Track{1,i};
+        thisCLST.Flag_OutboundSpk_Even_B4{1,1}=thisCLST.Flag_OutboundSpk_Even_B4{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         if ~sum(i==TrialNumber.VL)
+%             thisCLST.Flag_VM_OutboundSpk_Even_B4{1,1}=thisCLST.Flag_OutboundSpk_Even_B4{1,1}|thisCLST.Flag_OutboundSpk_Individual{1,i};
+%         end
+    end
+end

Matlab_function/GetSpikePosition.m

@@ -0,0 +1,32 @@
+function thisCLST = GetSpikePosition(thisCLST, Pos)
+rowSpike=length(thisCLST.Timestamp);
+row=length(Pos.t);
+thisCLST.x=zeros(rowSpike,1);
+thisCLST.y=zeros(rowSpike,1);
+thisCLST.t=zeros(rowSpike,1);
+thisCLST.SmoothingVelocity=zeros(rowSpike,1);
+
+k=1;
+for i=1:rowSpike
+    i
+    for j=k:row-3
+        if (thisCLST.Timestamp(i,1)>=Pos.t(j))&&(thisCLST.Timestamp(i,1)<Pos.t(j+1))
+            if(thisCLST.Timestamp(i,1)-Pos.t(j))>(Pos.t(j+1)-thisCLST.Timestamp(i,1))
+                thisCLST.x(i)=Pos.x(j+1);
+                thisCLST.y(i)=Pos.y(j+1);
+                thisCLST.t(i)=Pos.t(j+1);
+                thisCLST.SmoothingVelocity(i)=Pos.SmoothingVelocity(j+1);
+                k=j;
+                break;
+            elseif (thisCLST.Timestamp(i,1)-Pos.t(j))<=(Pos.t(j+1)-thisCLST.Timestamp(i,1))
+                thisCLST.x(i)=Pos.x(j);
+                thisCLST.y(i)=Pos.y(j);
+                thisCLST.t(i)=Pos.t(j);
+                thisCLST.SmoothingVelocity(i)=Pos.SmoothingVelocity(j);
+                k=j;
+                break;
+            end
+            
+        end
+    end
+end

Matlab_function/GetWeightedEuclideanDistance.m

@@ -0,0 +1,6 @@
+function WeightedEUdist=GetWeightedEuclideanDistance(x,y,z)
+SumOfSquare=0;
+for j=1:length(x)
+    SumOfSquare=SumOfSquare+z(j)*0.01*(x(j)-y(j))^2;
+end
+WeightedEUdist=sqrt(SumOfSquare);

Matlab_function/Jin_MeanSTE_Line.m

@@ -0,0 +1,54 @@
+function [Mean, STD, STE] = Jin_MeanSTE_Line(x,y,varargin)
+y = GetExcludeNan(y);
+
+% markersize=25;
+markersize=5;
+linewidth=1;
+% Example)
+% Width=0.3;
+% Color.color=53; Color.alpha=0.4;
+% h = JMeanSTE_Line(1,a.data1,Width,Color)
+
+if nargin == 2
+    Width = 0.25;
+    Color.color=4; Color.alpha=1;
+end
+
+if nargin == 3
+    Color = varargin{1};
+    Width =  0.25;
+end
+
+if nargin == 4
+    Width = varargin{1};
+    Color = varargin{2};
+end
+
+
+c=mapcolor(100,0);
+Mean=nanmean(y);
+STD=nanstd(y);
+STE=nanstd(y)/sqrt(length(y));
+
+x1=x-Width;
+x2=x+Width;
+y1=0;
+y2=Mean;
+
+% p1=plot(x,Mean,'.'); p1.MarkerSize=markersize; p1.Color=c(Color.color,:);
+% Std
+if Mean >= 0
+    l1=line([x x],[(Mean-STE) (Mean+STE)],'LineWidth',1,'LineStyle','-'); l1.Color=c(Color.color,:); l1.LineWidth=linewidth;
+    l1=line([(x-Width/3) (x+Width/3)],[(Mean+STE) (Mean+STE)],'LineWidth',1,'LineStyle','-'); l1.Color=c(Color.color,:); l1.LineWidth=linewidth;
+    l1=line([(x-Width/3) (x+Width/3)],[(Mean-STE) (Mean-STE)],'LineWidth',1,'LineStyle','-'); l1.Color=c(Color.color,:); l1.LineWidth=linewidth;
+elseif Mean < 0
+    l1=line([x x],[Mean+STE (Mean-STE)],'LineWidth',1,'LineStyle','-'); l1.Color=c(Color.color,:);
+    l1=line([(x-Width/3) (x+Width/3)],[(Mean+STE) (Mean+STE)],'LineWidth',1,'LineStyle','-'); l1.Color=c(Color.color,:); l1.LineWidth=linewidth;
+    l1=line([(x-Width/3) (x+Width/3)],[(Mean-STE) (Mean-STE)],'LineWidth',1,'LineStyle','-'); l1.Color=c(Color.color,:);
+end
+
+g=gca;
+g.YColor=[0 0 0];
+g.XColor=[0 0 0];
+g.ZColor=[0 0 0];
+g.LineWidth=0.6;

Matlab_function/PeriEventTimeHistogram_Spike.m

@@ -0,0 +1,104 @@
+function [PETHspike]=PeriEventTimeHistogram_Spike(TrialNumber, thisCLST, Time, RasterInfo)
+    rowSpike=length(thisCLST.Timestamp);
+
+    PETHspike.Choice=zeros(TrialNumber.Trial,1);
+    for j=1:TrialNumber.Trial
+        h=1;
+        for i=1:rowSpike
+            if(thisCLST.Timestamp(i)>=Time.Choice(j,1)-RasterInfo.Left)&&(thisCLST.Timestamp(i)<Time.Choice(j,1)+RasterInfo.Right)
+                PETHspike.Choice(j,h)=thisCLST.Timestamp(i);   
+                PETHspike.Choice_tagging(i)=1; % 180402 tagging spikes in the raster plot.
+                h=h+1;
+            end
+        end
+    end
+
+    %%
+    [Row_Spk_First,Col_Spk_First]=size(PETHspike.Choice);
+    for j=1:Row_Spk_First
+        for i=1:Col_Spk_First
+            if (PETHspike.Choice(j,i)==0)
+                PETHspike.Choice(j,i)=NaN;
+            end
+        end
+    end
+    for j=1:Row_Spk_First
+        for i=1:Col_Spk_First
+            PETHspike.Choice(j,i)=PETHspike.Choice(j,i)-Time.Choice(j,1);
+        end
+    end
+    
+    PETHspike.Return=zeros(TrialNumber.Trial,1);
+    for j=1:TrialNumber.Trial
+        h=1;
+        for i=1:rowSpike
+            if(thisCLST.Timestamp(i)>=Time.ReturnEntry(j,1)-RasterInfo.Left)&&(thisCLST.Timestamp(i)<Time.ReturnEntry(j,1)+RasterInfo.Right)
+                PETHspike.Return(j,h)=thisCLST.Timestamp(i);   
+                PETHspike.Return_tagging(i)=1; % 180402 tagging spikes in the raster plot.
+                h=h+1;
+            end
+        end
+    end
+
+    %%
+    [Row_Spk_First,Col_Spk_First]=size(PETHspike.Return);
+    for j=1:Row_Spk_First
+        for i=1:Col_Spk_First
+            if (PETHspike.Return(j,i)==0)
+                PETHspike.Return(j,i)=NaN;
+            end
+        end
+    end
+    for j=1:Row_Spk_First
+        for i=1:Col_Spk_First
+            PETHspike.Return(j,i)=PETHspike.Return(j,i)-Time.ReturnEntry(j,1);
+        end
+    end
+    
+end
+
+
+% Block 1 ~ 4
+% for i = 1:Row_Spk_First  
+%     if i <= TrialNumber.Odd(1)
+%         for j=1:Col_Spk_First
+%             if ~isnan(Alt_Raster_Spk_FirstRoute_10(i,j)) && Alt_Raster_Spk_FirstRoute_10(i,j)>-RasterInfo.RasterLeft*RasterResolution && Alt_Raster_Spk_FirstRoute_10(i,j)<RasterRight*RasterResolution
+%                 if Alt_Raster_Spk_FirstRoute_10(i,j)<0
+%                     Raster_First_Binning(1,(RasterInfo.RasterLeft*RasterResolution-fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))=Raster_First_Binning(1,RasterInfo.RasterLeft*RasterResolution-(fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))+1;
+%                 elseif Alt_Raster_Spk_FirstRoute_10(i,j)>=0
+%                     Raster_First_Binning(1,((RasterInfo.RasterLeft*RasterResolution+1)+fix(Alt_Raster_Spk_FirstRoute_10(i,j))))=Raster_First_Binning(1,(RasterInfo.RasterLeft*RasterResolution+1)+(fix(Alt_Raster_Spk_FirstRoute_10(i,j))))+1;
+%                 end
+%             end
+%         end
+%     elseif ~sum(i==RewardOmissionSecond)&& i>block(2)/2 && i<=block(3)/2
+%         for j=1:Col_Spk_First
+%             if ~isnan(Alt_Raster_Spk_FirstRoute_10(i,j)) && Alt_Raster_Spk_FirstRoute_10(i,j)>-RasterInfo.RasterLeft*RasterResolution && Alt_Raster_Spk_FirstRoute_10(i,j)<RasterRight*RasterResolution
+%                 if Alt_Raster_Spk_FirstRoute_10(i,j)<0
+%                     Raster_First_Binning(2,(RasterInfo.RasterLeft*RasterResolution-fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))=Raster_First_Binning(2,RasterInfo.RasterLeft*RasterResolution-(fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))+1;
+%                 elseif Alt_Raster_Spk_FirstRoute_10(i,j)>=0
+%                     Raster_First_Binning(2,((RasterInfo.RasterLeft*RasterResolution+1)+fix(Alt_Raster_Spk_FirstRoute_10(i,j))))=Raster_First_Binning(2,(RasterInfo.RasterLeft*RasterResolution+1)+(fix(Alt_Raster_Spk_FirstRoute_10(i,j))))+1;
+%                 end
+%             end
+%         end 
+%     elseif ~sum(i==RewardOmissionSecond)&& i>block(3)/2 && i<=block(4)/2
+%         for j=1:Col_Spk_First
+%             if ~isnan(Alt_Raster_Spk_FirstRoute_10(i,j)) && Alt_Raster_Spk_FirstRoute_10(i,j)>-RasterInfo.RasterLeft*RasterResolution && Alt_Raster_Spk_FirstRoute_10(i,j)<RasterRight*RasterResolution
+%                 if Alt_Raster_Spk_FirstRoute_10(i,j)<0
+%                     Raster_First_Binning(3,(RasterInfo.RasterLeft*RasterResolution-fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))=Raster_First_Binning(3,RasterInfo.RasterLeft*RasterResolution-(fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))+1;
+%                 elseif Alt_Raster_Spk_FirstRoute_10(i,j)>=0
+%                     Raster_First_Binning(3,((RasterInfo.RasterLeft*RasterResolution+1)+fix(Alt_Raster_Spk_FirstRoute_10(i,j))))=Raster_First_Binning(3,(RasterInfo.RasterLeft*RasterResolution+1)+(fix(Alt_Raster_Spk_FirstRoute_10(i,j))))+1;
+%                 end
+%             end
+%         end 
+%     elseif ~sum(i==RewardOmissionSecond)&& i>block(4)/2 && i<=block(5)/2
+%         for j=1:Col_Spk_First
+%             if ~isnan(Alt_Raster_Spk_FirstRoute_10(i,j)) && Alt_Raster_Spk_FirstRoute_10(i,j)>-RasterInfo.RasterLeft*RasterResolution && Alt_Raster_Spk_FirstRoute_10(i,j)<RasterRight*RasterResolution
+%                 if Alt_Raster_Spk_FirstRoute_10(i,j)<0
+%                     Raster_First_Binning(4,(RasterInfo.RasterLeft*RasterResolution-fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))=Raster_First_Binning(4,RasterInfo.RasterLeft*RasterResolution-(fix(Alt_Raster_Spk_FirstRoute_10(i,j))*(-1)))+1;
+%                 elseif Alt_Raster_Spk_FirstRoute_10(i,j)>=0
+%                     Raster_First_Binning(4,((RasterInfo.RasterLeft*RasterResolution+1)+fix(Alt_Raster_Spk_FirstRoute_10(i,j))))=Raster_First_Binning(4,(RasterInfo.RasterLeft*RasterResolution+1)+(fix(Alt_Raster_Spk_FirstRoute_10(i,j))))+1;
+%                 end
+%             end
+%         end 
+%     end
+

Matlab_function/Stat_ANOVA1.m

@@ -0,0 +1,5 @@
+function [Pvalue, result, gnames] = Stat_ANOVA1(Y,GROUP) 
+    [Pvalue.Main_OneWay,result,stats] = anova1(Y,GROUP,'off');
+    [Pvalue.MultipleComparison_OneWay,~,~,gnames] = multcompare(stats, 'display', 'off');
+end
+

Matlab_function/Stat_ANOVA2.m

@@ -0,0 +1,33 @@
+function [Pvalue, result, gnames] = Stat_ANOVA2(Y,GROUP,NAME) 
+        
+% Example 1
+% a=[1;2;3;4;5]
+% b=[2;2;3;4;4]
+% c=[1;3;3;2;2]
+% d=[2;3;4;4;1]
+% A=[1;2;3;4;5]
+% B=[2;2;3;4;4]
+% C=[1;3;3;2;2]
+% D=[2;3;4;4;1]
+% Y = [a; b; c; d;...
+%      A; B; C; D];
+% Species = [GetGroupingVar(a,1); GetGroupingVar(b,1); GetGroupingVar(c,1); GetGroupingVar(d,1); ...
+%     GetGroupingVar(A,2); GetGroupingVar(B,2); GetGroupingVar(C,2); GetGroupingVar(D,2)];
+% Block = [GetGroupingVar(a,1); GetGroupingVar(b,2); GetGroupingVar(c,3); GetGroupingVar(d,4); ...
+%     GetGroupingVar(A,1); GetGroupingVar(B,2); GetGroupingVar(C,3); GetGroupingVar(D,4)];
+% NAME = {'Species', 'Block'};
+% [Pvalue, result, gnames] = Stat_ANOVA2(Y,[Species Block],NAME)
+
+% Example 2
+% Y = [dHP_Odd; dHP_Even; iHP_Odd; iHP_Even];
+% Region = [GetGroupingVar(dHP_Odd,1); GetGroupingVar(dHP_Even,1); ...
+%     GetGroupingVar(iHP_Odd,2); GetGroupingVar(iHP_Even,2)];
+% Trajectory = [GetGroupingVar(dHP_Odd,1); GetGroupingVar(dHP_Even,2);...
+%     GetGroupingVar(iHP_Odd,1); GetGroupingVar(iHP_Even,2)];
+% NAME = {'Region', 'Trajectory'};
+% [Pvalue, result, gnames] = Stat_ANOVA2(Y,[Region Trajectory],NAME)
+
+
+[Pvalue.Main_TwoWay,result,stats] = anovan(Y,GROUP,'model','interaction','varnames',NAME,'display','off');
+[Pvalue.MultipleComparison, ~, ~, gnames] = multcompare(stats,'Dimension',[1,2],'display','off');
+end

Matlab_function/Stat_ANOVA2_Mixed.m

@@ -0,0 +1,196 @@
+function [SSQs, DFs, MSQs, Fs, Ps]=Stat_ANOVA2_Mixed(X,suppress_output)
+% simple function for mixed (between- and within-subjects) ANOVA
+% 
+% Based loosely on BWAOV2 (http://www.mathworks.com/matlabcentral/fileexchange/5579-bwaov2) by Antonio Trujillo-Ortiz
+%  (insofar as I used that function to figure out the basic equations, as it is apparently very hard to find documentation
+%  on mixed-model ANOVAs on the Internet). However, the code is all original.
+% 
+% The major advantage of this function over the existing BWAOV2 is that it corrects a bug that occurs when the groups
+%  have unequal numbers of subjects, as pointed out in the Matlab Central File Exchange by Jaewon Hwang. The code is also,
+%  at least in my opinion, much cleaner.
+% 
+% At present this function only supports mixed models with a single between-subjects factor and a single within-subjects
+%  (repeated measures) factor, each with as many levels as you like. I would be happy to add more bells and whistles in
+%  future editions, such as the ability to define multiple factors, apply Mauchly's test and add in non-sphericity 
+%  corrections when necessary, etc. I'm a better programmer than I am a statistician, though, so if anyone out there would 
+%  like to lend a hand with the math (e.g., feed me the equations for these features), I'd be more than happy to implement 
+%  them. Email matthew DOT r DOT johnson AT aya DOT yale DOT edu
+% 
+% Also feel free to modify this file for your own purposes and upload your changes to the Matlab Central File Exchange if 
+%  you like.
+% 
+% I have checked this function against the example data in David M. Lane's HyperStat online textbook, which is the same 
+%  data that breaks BWAOV2 (http://davidmlane.com/hyperstat/questions/Chapter_14.html, question 6). I have also checked it
+%  against SPSS and gotten identical results. However, I haven't tested every possible case so bugs may remain. Use at 
+%  your own risk. If you find bugs and let me know, I'm happy to try to fix them.
+% 
+% ===============
+%      USAGE
+% ===============
+% 
+% Inputs:
+% 
+% X: design matrix with four columns (future versions may allow different input configurations)
+%     - first column  (i.e., X(:,1)) : all dependent variable values
+%     - second column (i.e., X(:,2)) : between-subjects factor (e.g., subject group) level codes (ranging from 1:L where 
+%         L is the # of levels for the between-subjects factor)
+%     - third column  (i.e., X(:,3)) : within-subjects factor (e.g., condition/task) level codes (ranging from 1:L where 
+%         L is the # of levels for the within-subjects factor)
+%     - fourth column (i.e., X(:,4)) : subject codes (ranging from 1:N where N is the total number of subjects)
+% 
+% suppress_output: defaults to 0 (meaning it displays the ANOVA table as output). If you don't want to display the table,
+%  just pass in a non-zero value
+% 
+% Outputs:
+% 
+% SSQs, DFs, MSQs, Fs, Ps : Sum of squares, degrees of freedom, mean squares, F-values, P-values. All the same values 
+%  that are shown in the table if you choose to display it. All will be cell matrices. Values within will be in the same 
+%  order that they are shown in the output table.
+% 
+% Enjoy! -MJ
+
+% Example by JSW
+% dHP_Pre=[PopulationAvgFr.Pre_dHP_PYR_1]';
+% dHP_Main=[PopulationAvgFr.Main_dHP_PYR_1]';
+% iHP_Pre=[PopulationAvgFr.Pre_iHP_PYR_1]';
+% iHP_Main=[PopulationAvgFr.Main_iHP_PYR_1]';
+% 
+% Y = [dHP_Pre; dHP_Main; iHP_Pre; iHP_Main];
+% BS = [GetGroupingVar(dHP_Pre,1); GetGroupingVar(dHP_Main,1); GetGroupingVar(iHP_Pre,2); GetGroupingVar(iHP_Main,2)];
+% WS = [GetGroupingVar(dHP_Pre,1); GetGroupingVar(dHP_Main,2); GetGroupingVar(iHP_Pre,1); GetGroupingVar(iHP_Main,2)];
+% S = [GetGroupingVar(dHP_Pre); GetGroupingVar(dHP_Main); GetGroupingVar(iHP_Pre)+length(dHP_Pre); GetGroupingVar(iHP_Main)+length(dHP_Pre)];
+% X = [Y BS WS S];    
+% [SSQs, DFs, MSQs, Fs, Ps] = Stat_ANOVA2_Mixed(X);
+if nargin < 1, 
+   error('No input');
+end;
+
+if nargin < 2 || isempty(suppress_output)
+    suppress_output=0;
+end
+
+all_dvs=X(:,1);
+all_bs_labels=X(:,2);
+all_ws_labels=X(:,3);
+all_subj_labels=X(:,4);
+
+bs_levels=sort(unique(all_bs_labels));
+ws_levels=sort(unique(all_ws_labels));
+subj_levels=sort(unique(all_subj_labels));
+n_bs_levels=length(bs_levels);
+n_ws_levels=length(ws_levels);
+n_subjects=length(subj_levels);
+
+if any(bs_levels(:)~=(1:n_bs_levels)') || any(ws_levels(:)~=(1:n_ws_levels)') || any(subj_levels(:)~=(1:n_subjects)')
+    error('Levels of factors (including subject labels) must be numbered 1:L (where L is the # of levels for that factor');
+end
+for i=1:n_bs_levels
+    for j=1:n_ws_levels
+        this_cell_inds=find(all_bs_labels==i & all_ws_labels==j);
+        if isempty(this_cell_inds)
+            error('At least one empty cell found');
+        end
+        n_subs_per_cell(i,j)=length(this_cell_inds); %#ok<AGROW>
+        cell_totals(i,j)=sum(all_dvs(this_cell_inds)); %#ok<AGROW>
+    end
+    
+    if any(n_subs_per_cell(i,:)~=n_subs_per_cell(i,1))
+        error('At least one subject missing at least one repeated measure (or is possibly entered more than once)');
+        %technically this is not a failsafe check, as it could be that subject ! has only conditions A & B,
+        %  whereas subject 2 has only condition C, which still gives equal values for all the conditions.
+        %  We'll double-check for that circumstance below
+    end
+end
+
+for k=1:n_subjects
+    this_subj_inds=find(all_subj_labels==k);
+    if length(this_subj_inds)~=n_ws_levels
+        %our second check for this issue
+        error('At least one subject missing at least one repeated measure (or is possibly entered more than once)');
+    end
+    subj_totals(k)=sum(all_dvs(this_subj_inds)); %#ok<AGROW>
+end
+
+correction_term = sum(all_dvs)^2 / length(all_dvs);
+SStot = sum(all_dvs.^2) - correction_term;
+%don't really need this for calculations, but can uncomment if we want to print
+% DFtot = length(all_dvs) - 1; %total degrees of freedom
+
+%subject "factor" (i.e. differences in subject means)
+SSsub = sum(subj_totals .^ 2)/n_ws_levels - correction_term;
+
+%between-subjects factor
+SStmp=[];
+for i=1:n_bs_levels
+    SStmp(i)=(sum(cell_totals(i,:))^2) / sum(n_subs_per_cell(i,:)); %#ok<AGROW>
+end
+SSbs    = sum(SStmp) - correction_term;
+DFbs    = n_bs_levels - 1;
+MSbs    = SSbs / DFbs;
+%error terms for between-subjects factor
+ERRbs   = SSsub - SSbs;
+DFERRbs = n_subjects - n_bs_levels;
+MSERRbs = ERRbs / DFERRbs;
+
+%correction with harmonic mean of cell sizes if cell sizes are not all equal
+n_subs_hm=harmmean(n_subs_per_cell(:));
+cell_totals_hm = (cell_totals ./ n_subs_per_cell) * n_subs_hm;
+correction_term_hm = sum(cell_totals_hm(:))^2 / (n_subs_hm * n_bs_levels * n_ws_levels);
+n_subs_per_cell_hm = ones(n_bs_levels,n_ws_levels) * n_subs_hm;
+
+%within-subjects factor
+SStmp=[];
+for j=1:n_ws_levels
+    SStmp(j)=(sum(cell_totals_hm(:,j))^2) ./ sum(n_subs_per_cell_hm(:,j)); %#ok<AGROW>
+end
+SSws  = sum(SStmp) - correction_term_hm;
+DFws  = n_ws_levels - 1;
+MSws  = SSws / DFws;
+
+%uncorrected version of within-subjects factor for calculating interaction
+SStmp=[];
+for j=1:n_ws_levels
+    SStmp(j)=(sum(cell_totals(:,j))^2) ./ sum(n_subs_per_cell(:,j)); %#ok<AGROW>
+end
+SSws_unc = sum(SStmp) - correction_term;
+
+%interaction of between-subjects and within-subjects factor
+SStmp = sum((cell_totals(:) .^ 2) ./ n_subs_per_cell(:));
+SSint = SStmp - SSbs - SSws_unc - correction_term;
+DFint = DFbs * DFws;
+MSint = SSint / DFint;
+
+%error terms (for both within-subjects factor and interaction)
+ERRws   = SStot - SSbs - ERRbs - SSws_unc - SSint;
+DFERRws = DFERRbs * DFws;
+MSERRws = ERRws / DFERRws;
+
+%F-values
+Fbs  = MSbs  / MSERRbs;
+Fws  = MSws  / MSERRws;
+Fint = MSint / MSERRws;
+
+%P-values
+Pbs  = 1-fcdf(Fbs, DFbs, DFERRbs);
+Pws  = 1-fcdf(Fws, DFws, DFERRws);
+Pint = 1-fcdf(Fint,DFint,DFERRws);
+
+if ~suppress_output
+    disp(' ');
+    disp(' ');
+    fprintf(1,'     Source                          SSq        df       MSq         F         p        \n');
+    fprintf(1,'-------------------------------------------------------------------------------------\n');
+    fprintf(1,'     Between-subjects factor   %9.3f %9i %9.3f %9.3f %9.7f\n',SSbs,DFbs,MSbs,Fbs,Pbs);
+    fprintf(1,'     Between-subjects error    %9.3f %9i %9.3f\n',ERRbs,DFERRbs,MSERRbs);
+    fprintf(1,'     Within-subjects factor    %9.3f %9i %9.3f %9.3f %9.7f\n',SSws,DFws,MSws,Fws,Pws);
+    fprintf(1,'     Within x between int.     %9.3f %9i %9.3f %9.3f %9.7f\n',SSint,DFint,MSint,Fint,Pint);
+    fprintf(1,'     Within-subjects error     %9.3f %9i %9.3f\n',ERRws,DFERRws,MSERRws);
+    disp(' ');
+    disp(' ');
+end
+
+SSQs = { SSbs; ERRbs;   SSws; SSint; ERRws   };
+DFs  = { DFbs; DFERRbs; DFws; DFint; DFERRws };
+MSQs = { MSbs; MSERRbs; MSws; MSint; MSERRws };
+Fs   = { Fbs;  [];      Fws;  Fint;  []      };
+Ps   = { Pbs;  [];      Pws;  Pint;  []      };

Matlab_function/abmFiringRateMap.m

@@ -0,0 +1,145 @@
+function [occMat spkMat rawMat skaggsrateMat] = abmFiringRateMap(tmazeMat, finAlzPosMat, imROW, imCOL, thisSCALE, samplingRate)
+%Generate firing rate map using Skaggs' adaptive binning methods
+%[occMat spkMat rawMat skaggsrateMat] = abmFiringRateMap(tmazeMat, finAlzPosMat, imROW, imCOL, thisSCALE, samplingRate)
+%
+%Input format
+%  tmazeMat	 [n x 3 matrix]									timestamp x-coord y-coord
+%  finAlzPosMat [n x 3 matrix]									timestamp x-coord y-coord
+%  imROW [1 x 1]												number of rows of ratemap
+%  imCOL [1 x 1]												number of columns of ratemap
+%  thisSCALE [1 x 1]											if you want to use original resolution, it should be 1. If you want to scale-down the resolusion by 10, it should be 10.
+%  samplingRate [1 x 1]											sampling rate for behavioral traces. The current settings use 30 Hz.
+%
+%Output format
+% occMat [imROW x imCOL matrix]									raw occupancy map
+% spkMat [imROW x imCOL matrix]									raw spike position map
+% rawMat [imROW x imCOL matrix; spkMat ./ (occMat .* samplingRate)]		raw firing rate map
+% skaggsrateMat [imROW x imCOL matrix]								Skaggs' ABM firing rate map
+%
+%Originally from VB codes which Inah Lee has.
+%Translate VB codes to matlab was done by [Jangjin Kim, July-13-2008]
+%Verification was done
+%reviced by SB 4/29/2013 : add 'to avoid error'
+%revised by JSW 4/4/2019 : add empty tmazeMat error
+
+    %to avoid error
+    tmazeMatSize = size(tmazeMat);
+    if tmazeMatSize(1) == 0 && tmazeMatSize(2) == 0
+        tmazeMat = zeros(0,3);
+    end
+
+    %pINDEX
+    pXCOORD = 2;																						%x coordinates [from Neuralynx]
+    pYCOORD = 3;																						%y coordinates [from Neuralynx]
+    alphaSET = .0001;
+
+    occMat = nan(imROW, imCOL);
+    spkMat = zeros(imROW, imCOL);
+    rawMat = zeros(imROW, imCOL);
+    skaggsrateMat = zeros(imROW, imCOL);
+
+    for rowRUN = 1:1:imROW
+        for colRUN = 1:1:imCOL
+            if ~isempty(tmazeMat) %revised by JSW 4/4/2019 : add empty tmazeMat error
+            numspk =  size(tmazeMat(tmazeMat(:, pXCOORD) ./ thisSCALE > (colRUN - 1) & tmazeMat(:, pXCOORD) ./ thisSCALE <= colRUN & ... 
+                                    tmazeMat(:, pYCOORD) ./ thisSCALE > (rowRUN - 1) & tmazeMat(:, pYCOORD) ./ thisSCALE <= rowRUN, :));
+                if numspk(1, 1) == 0
+                    spkMat(rowRUN, colRUN) = 0;
+                else
+                    spkMat(rowRUN, colRUN) = numspk(1, 1);
+                end	%numspk(1, 1) == 0
+            else
+                spkMat(rowRUN, colRUN) = 0;
+            end
+            if ~isempty(finAlzPosMat)
+                numocc = size(finAlzPosMat(finAlzPosMat(:, pXCOORD) ./ thisSCALE > (colRUN - 1) & finAlzPosMat(:, pXCOORD) ./ thisSCALE <= colRUN & ... 
+                                    finAlzPosMat(:, pYCOORD) ./ thisSCALE > (rowRUN - 1) & finAlzPosMat(:, pYCOORD) ./ thisSCALE <= rowRUN, :));
+                if numocc(1, 1) == 0
+                    occMat(rowRUN, colRUN) = nan;
+                else
+                    occMat(rowRUN, colRUN) = numocc(1, 1);
+                end	%numocc(1, 1) == 0
+            else
+                occMat(rowRUN, colRUN) = nan;
+            end
+        end	%colRUN = 1:1:imCOL
+    end	%rowRUN = 1:1:imROW
+
+    rawMat = spkMat ./ occMat .* samplingRate;
+
+    if sum(sum(spkMat)) > 0
+        nanIND = find(isnan(occMat));
+        occMat(isnan(occMat)) = 0;
+        i = -20:20; j = -20:19;
+        Ti = zeros(28, 400);
+        Tj = zeros(28, 400);
+
+        for colRUN = 1:1:imCOL
+            for rowRUN = 1:1:imROW
+                if occMat(rowRUN, colRUN) > 0
+                    %Equivalent to BuildRTable
+                    Nr = zeros(1, 400);
+
+                    for ii = 1:1:41
+                        for jj = 1:1:40
+                            %r(1, (ii - 1) * 41 + jj) = i(1, ii)^2 + j(1, jj)^2; 
+                            %r(ii, jj) = i(1, ii)^2 + j(1, jj)^2; 
+                            r = i(1, ii)^2 + j(1, jj)^2 + 1;
+                            if r <= 400
+                                N = Nr(r);
+                                Ti(N + 1, r) = i(1, ii);
+                                Tj(N + 1, r) = j(1, jj);
+                                Nr(r) = Nr(r) + 1;
+                            end	%r <= 400
+                        end	%jj = 1:1:40
+                    end	%ii = 1:1:41
+
+                    nSPK = spkMat(rowRUN, colRUN);
+                    nOCC = occMat(rowRUN, colRUN);
+
+                    lambda = nSPK / nOCC * samplingRate;
+
+                    for colRUN2 = 1:1:imCOL								%i
+                        for rowRUN2 = 1:1:imROW							%j
+                            if occMat(rowRUN2, colRUN2) == 0
+                                skaggsrateMat(rowRUN2, colRUN2) = nan;
+                            else
+                                EnoughPoints = false;
+                                occCount = 1;						%Assume there's at least 1 spk and 1 occ
+                                spkCount = 1;
+                                r = 0;
+                                while (r <= 200 & EnoughPoints == false) % r <= 200 original
+                                    for h = 1:1:Nr(r + 1)
+                                        ii = (colRUN2 - 1) + Ti(h, r + 1);
+                                        jj = (rowRUN2 - 1) + Tj(h, r + 1);	%-1 to compensate VB code
+                                        if (ii >= 0 & ii < 64) & (jj >= 0 & jj < 48)
+                                            occCount = occCount + occMat(jj + 1, ii + 1);
+                                            spkCount = spkCount + spkMat(jj + 1, ii + 1);
+                                        end	%(ii >= 0 & ii < 64) & (jj >= 0 & jj < 48)
+                                    end	%h = 1:1:Nr(r + 1)
+                                    if alphaSET^2 * occCount^2 * r * spkCount > 1
+                                        EnoughPoints = true;
+                                    end	%alphaSET^2 * occCount^2 * r * spkCount > 1
+                                    r = r + 1;
+                                end	%(r <= 200 & EnoughPoints == false)
+                                if occCount < samplingRate / 2
+                                    skaggsrateMat(rowRUN2, colRUN2) = 0;
+                                else
+                                    skaggsrateMat(rowRUN2, colRUN2) = spkCount / occCount * samplingRate;
+                                end	%occCount < samplingRate / 2
+                            end	%occMat(colRUN2, rowRUN2) == 0
+                            if occMat(rowRUN2, colRUN2) ~= 0 & isnan(skaggsrateMat(rowRUN2, colRUN2))
+                                rowRUN2 = rowRUN2;
+                            end	%occMat(rowRUN2, colRUN2) ~= 0 & isnan(skaggsrateMat(rowRUN2, colRUN2))
+                            if rowRUN2 == 16 & colRUN2 == 39				%Add 1 to compensate VB indices
+                                rowRUN2 = rowRUN2;
+                            end	%rowRUN2 == 16 & colRUN2 == 39
+                        end	%rowRUN2 = 1:1:imROW
+                    end	%colRUN2 = 1:1:imCOL
+                end	%occMat(rowRUN, colRUN) > 0
+            end	%rowRUN = 1:1:imROW
+        end	%colRUN = 1:1:imCOL
+        occMat(nanIND) = nan;
+    else
+        skaggsrateMat = rawMat;
+    end	%sum(sum(spkMat)) > 0

Matlab_function/distinguishable_colors.m

@@ -0,0 +1,152 @@
+function colors = distinguishable_colors(n_colors,bg,func)
+% DISTINGUISHABLE_COLORS: pick colors that are maximally perceptually distinct
+%
+% When plotting a set of lines, you may want to distinguish them by color.
+% By default, Matlab chooses a small set of colors and cycles among them,
+% and so if you have more than a few lines there will be confusion about
+% which line is which. To fix this problem, one would want to be able to
+% pick a much larger set of distinct colors, where the number of colors
+% equals or exceeds the number of lines you want to plot. Because our
+% ability to distinguish among colors has limits, one should choose these
+% colors to be "maximally perceptually distinguishable."
+%
+% This function generates a set of colors which are distinguishable
+% by reference to the "Lab" color space, which more closely matches
+% human color perception than RGB. Given an initial large list of possible
+% colors, it iteratively chooses the entry in the list that is farthest (in
+% Lab space) from all previously-chosen entries. While this "greedy"
+% algorithm does not yield a global maximum, it is simple and efficient.
+% Moreover, the sequence of colors is consistent no matter how many you
+% request, which facilitates the users' ability to learn the color order
+% and avoids major changes in the appearance of plots when adding or
+% removing lines.
+%
+% Syntax:
+%   colors = distinguishable_colors(n_colors)
+% Specify the number of colors you want as a scalar, n_colors. This will
+% generate an n_colors-by-3 matrix, each row representing an RGB
+% color triple. If you don't precisely know how many you will need in
+% advance, there is no harm (other than execution time) in specifying
+% slightly more than you think you will need.
+%
+%   colors = distinguishable_colors(n_colors,bg)
+% This syntax allows you to specify the background color, to make sure that
+% your colors are also distinguishable from the background. Default value
+% is white. bg may be specified as an RGB triple or as one of the standard
+% "ColorSpec" strings. You can even specify multiple colors:
+%     bg = {'w','k'}
+% or
+%     bg = [1 1 1; 0 0 0]
+% will only produce colors that are distinguishable from both white and
+% black.
+%
+%   colors = distinguishable_colors(n_colors,bg,rgb2labfunc)
+% By default, distinguishable_colors uses the image processing toolbox's
+% color conversion functions makecform and applycform. Alternatively, you
+% can supply your own color conversion function.
+%
+% Example:
+%   c = distinguishable_colors(25);
+%   figure
+%   image(reshape(c,[1 size(c)]))
+%
+% Example using the file exchange's 'colorspace':
+%   func = @(x) colorspace('RGB->Lab',x);
+%   c = distinguishable_colors(25,'w',func);
+
+% Copyright 2010-2011 by Timothy E. Holy
+
+  % Parse the inputs
+  if (nargin < 2)
+    bg = [1 1 1];  % default white background
+  else
+    if iscell(bg)
+      % User specified a list of colors as a cell aray
+      bgc = bg;
+      for i = 1:length(bgc)
+	bgc{i} = parsecolor(bgc{i});
+      end
+      bg = cat(1,bgc{:});
+    else
+      % User specified a numeric array of colors (n-by-3)
+      bg = parsecolor(bg);
+    end
+  end
+  
+  % Generate a sizable number of RGB triples. This represents our space of
+  % possible choices. By starting in RGB space, we ensure that all of the
+  % colors can be generated by the monitor.
+  n_grid = 30;  % number of grid divisions along each axis in RGB space
+  x = linspace(0,1,n_grid);
+  [R,G,B] = ndgrid(x,x,x);
+  rgb = [R(:) G(:) B(:)];
+  if (n_colors > size(rgb,1)/3)
+    error('You can''t readily distinguish that many colors');
+  end
+  
+  % Convert to Lab color space, which more closely represents human
+  % perception
+  if (nargin > 2)
+    lab = func(rgb);
+    bglab = func(bg);
+  else
+    C = makecform('srgb2lab');
+    lab = applycform(rgb,C);
+    bglab = applycform(bg,C);
+  end
+
+  % If the user specified multiple background colors, compute distances
+  % from the candidate colors to the background colors
+  mindist2 = inf(size(rgb,1),1);
+  for i = 1:size(bglab,1)-1
+    dX = bsxfun(@minus,lab,bglab(i,:)); % displacement all colors from bg
+    dist2 = sum(dX.^2,2);  % square distance
+    mindist2 = min(dist2,mindist2);  % dist2 to closest previously-chosen color
+  end
+  
+  % Iteratively pick the color that maximizes the distance to the nearest
+  % already-picked color
+  colors = zeros(n_colors,3);
+  lastlab = bglab(end,:);   % initialize by making the "previous" color equal to background
+  for i = 1:n_colors
+    dX = bsxfun(@minus,lab,lastlab); % displacement of last from all colors on list
+    dist2 = sum(dX.^2,2);  % square distance
+    mindist2 = min(dist2,mindist2);  % dist2 to closest previously-chosen color
+    [~,index] = max(mindist2);  % find the entry farthest from all previously-chosen colors
+    colors(i,:) = rgb(index,:);  % save for output
+    lastlab = lab(index,:);  % prepare for next iteration
+  end
+end
+
+function c = parsecolor(s)
+  if ischar(s)
+    c = colorstr2rgb(s);
+  elseif isnumeric(s) && size(s,2) == 3
+    c = s;
+  else
+    error('MATLAB:InvalidColorSpec','Color specification cannot be parsed.');
+  end
+end
+
+function c = colorstr2rgb(c)
+  % Convert a color string to an RGB value.
+  % This is cribbed from Matlab's whitebg function.
+  % Why don't they make this a stand-alone function?
+  rgbspec = [1 0 0;0 1 0;0 0 1;1 1 1;0 1 1;1 0 1;1 1 0;0 0 0];
+  cspec = 'rgbwcmyk';
+  k = find(cspec==c(1));
+  if isempty(k)
+    error('MATLAB:InvalidColorString','Unknown color string.');
+  end
+  if k~=3 || length(c)==1,
+    c = rgbspec(k,:);
+  elseif length(c)>2,
+    if strcmpi(c(1:3),'bla')
+      c = [0 0 0];
+    elseif strcmpi(c(1:3),'blu')
+      c = [0 0 1];
+    else
+      error('MATLAB:UnknownColorString', 'Unknown color string.');
+    end
+  end
+end

Matlab_function/mapcolor.m

@@ -0,0 +1,19 @@
+function c=mapcolor(i,flag)
+% flag : 1 
+  c = distinguishable_colors(i);
+  if flag==1
+      fig=figure;
+      fig.Position=[0 0 2300 1000];
+      image(reshape(c,[10 size(c,1)/10 3]));
+      
+%       EPS file save    
+%         filename_ai=['E:\EphysAnalysis\Result\ColorMap.eps'];
+%         print( gcf, '-painters', '-r300', filename_ai, '-depsc');
+% 
+        SaveRoot=['E:\EphysAnalysis\Result'];
+        FileName=['ColorMap.png'];
+        SaveFigure(gcf,SaveRoot, FileName);
+
+
+  end
+end