---
title: "Fontanier et al - Timescale extraction and descriptive statistics"
author: "Vincent Fontanier, Matthieu Sarazin, Frederic M. Stoll, Bruno Delord, and Emmanuel Procyk"
date: "16/07/2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

This document provides access to timescale extraction and their descriptive statistics as described in the manuscript "Local determinants of network dynamics underly the temporal signature of cognitive states in the frontal cortex".

Data availablity :
Spike timeseries and units information are accessible in the Zenodo folder

Tested the 16/07/2020

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252    LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C                   LC_TIME=French_France.1252    

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lme4_1.1-23         Matrix_1.2-17       minpack.lm_1.2-1    car_3.0-8           carData_3.0-4      
 [6] knitr_1.28          broom_0.5.6         ggsci_2.9           dendextend_1.13.4   NbClust_3.0        
[11] timeSeries_3062.100 timeDate_3043.102   ggExtra_0.9         scales_1.1.1        pracma_2.2.9       
[16] RColorBrewer_1.1-2  reshape2_1.4.4      factoextra_1.0.7    FactoMineR_2.3      lattice_0.20-38    
[21] zoo_1.8-8           MASS_7.3-51.4       gridExtra_2.3       dplyr_1.0.0         tidyr_1.1.0        
[26] ggplot2_3.3.1      

loaded via a namespace (and not attached):
 [1] nlme_3.1-140         tools_3.6.1          backports_1.1.7      R6_2.4.1             mgcv_1.8-28         
 [6] colorspace_1.4-1     withr_2.2.0          tidyselect_1.1.0     curl_4.3             compiler_3.6.1      
[11] cli_2.0.2            flashClust_1.01-2    labeling_0.3         stringr_1.4.0        digest_0.6.25       
[16] foreign_0.8-71       minqa_1.2.4          rio_0.5.16           pkgconfig_2.0.3      htmltools_0.4.0     
[21] highr_0.8            fastmap_1.0.1        rlang_0.4.6          readxl_1.3.1         rstudioapi_0.11     
[26] shiny_1.4.0.2        generics_0.0.2       farver_2.0.3         zip_2.0.4            magrittr_1.5        
[31] leaps_3.1            Rcpp_1.0.4.6         munsell_0.5.0        fansi_0.4.1          abind_1.4-5         
[36] viridis_0.5.1        lifecycle_0.2.0      scatterplot3d_0.3-41 stringi_1.4.6        yaml_2.2.1          
[41] plyr_1.8.6           promises_1.1.1       ggrepel_0.8.2        forcats_0.5.0        crayon_1.3.4        
[46] miniUI_0.1.1.1       haven_2.3.1          splines_3.6.1        hms_0.5.3            pillar_1.4.4        
[51] boot_1.3-22          glue_1.4.1           data.table_1.12.8    vctrs_0.3.1          nloptr_1.2.2.1      
[56] httpuv_1.5.4         cellranger_1.1.0     gtable_0.3.0         purrr_0.3.4          assertthat_0.2.1    
[61] xfun_0.14            openxlsx_4.1.5       mime_0.9             xtable_1.8-4         later_1.1.0.1       
[66] viridisLite_0.3.0    tibble_3.0.1         cluster_2.1.0        statmod_1.4.34       ellipsis_0.3.1   


\newpage


Libraries loading
```{r Libraries loading, message=FALSE, warning=FALSE, include=T, echo=FALSE}
library(ggplot2)
theme_set(theme_bw())
library(tidyr)
library(dplyr)
library(gridExtra)
# See https://cran.r-project.org/web/packages/egg/vignettes/Ecosystem.html for grobs arrange
library(MASS) #caution this will mask dplyr select so we will need to dplyr::select()
library(zoo)
library(lattice)
library(FactoMineR)
library(factoextra)
library(reshape2)
library(RColorBrewer)
library(pracma)
library(scales)
library(grid)
library(gridExtra)
library(ggExtra)
library(timeSeries)
library(NbClust)
library(dendextend)
library("ggsci")
library(broom)
library(knitr)
library(car)
library(minpack.lm)
library(lme4)
```

Path and useful variables

For non-parametric statistical testing with interaction we will use the BLOM transformation method 
Briefly it is a subcase of Rank-Based Inverse Normal Transformations (refered as INT)
It is described by the transformation:
y <- qnorm((r-c)/(N-2c+1))
where qnorm is the standard normal quantile function
r is the rank of your variable
c = 3/8, the Blom constant
N is the sample size

```{r}
blom <- function(r,N){
  blomfactor <- 3/8;
  out <- qnorm((r-blomfactor)/(N-2*blomfactor+1))
  return(out)
}
```

```{r}
# Compute sem and confidence intervals
sem <- function(x,na.rm){
  if(missing(na.rm)){
    na.rm <- F;
  }
  if(na.rm==T){x<-x[!is.na(x)]};
  out <- sd(x)/sqrt(length(x));
  return(out)
}
cilow <- function(x,na.rm) {
  if(missing(na.rm)){
    na.rm <- F;
  }
  out <- mean(x,na.rm=na.rm)-2*sem(x,na.rm=na.rm);
  return(out)
}
ciup <- function(x,na.rm) {
  if(missing(na.rm)){
    na.rm <- F;
  }
  out <- mean(x,na.rm=na.rm)+2*sem(x,na.rm=na.rm);
  return(out)
}

#Useful function to bin continuous numerical data
rm(x,my.binsize,x.min,x.max,x.breaks,x.labs,xbin)

fun.bin <- function(x,my.binsize) {
  
  x.min <- floor(min(x));
  x.max <- ceiling(max(x));
  
  x.breaks <- seq(from=x.min,to=x.max+my.binsize,by=my.binsize);
  x.labs <- colMeans(rbind(x.breaks[-length(x.breaks)],c(x.breaks[-1])));
  
  xbin <- cut(x,breaks = x.breaks,labels = x.labs,right = F);
  xbin <- as.numeric(as.character(xbin));

  return(xbin)
}

if (Sys.info()["nodename"] == "PC_MANU_ATIV9"){
  path <- "C://Users/eprocyk/IBox/TimeScale/"
  
}else if(Sys.info()["nodename"] =="DESKTOP-MANU"){
   path <- "C://Users/E Procyk/IBOX/TimeScale/"
   
}else if(Sys.info()["nodename"] =="VINCENT-VAIO"){
   path <- "C://Users/Vincent/Nextcloud/TimeScale//"

}else if(Sys.info()["nodename"] =="R31STAGIAIRE"){
   path <- "C://Users/Julien/Nextcloud/TimeScale/"
}else if(Sys.info()["nodename"] =="DESKTOP-T35MINL"){
   path <- "C://Users/vince/Nextcloud/TimeScale/"
}else{
  path  <- "D://My Dropbox/My Dropbox/Main_drop/R/Check/TimeScales/"
}

theme_fig <- function (rat) { 
    theme_bw() %+replace% 
        theme(aspect.ratio=rat,panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),legend.position = "none")
}

colMCC = "#666666"
colLPFC = "#5FC6F1"

colArea <- c(colMCC,colLPFC);
```

Load information about the units
```{r}
##  Information about units
load(file=paste0(path,"neurons_info.Rdata"))
```


# Timescales extraction from spike timeseries (spike autocorrelogram method)

Briefly for each unit the steps are:

\begin{enumerate}

\item Compute the autoc

\begin{enumerate}
\item Compute the lagged differences of the spiketimeseries from 1st to 100th lag order
\item Sum up lagged differences from 0 to 1000ms in 3.33ms wide bins
\item Remove the 10 first ms of the Autoc
\item Compute the probability density function of the autocorrelogram defined as :
Pi = Xi/S(X) * 1/deltaT

with Pi the PDF value for bin i, Xi the count of bin i, S(X) the sum of counts and deltaT the temporal bin width

(minmax normalize bincounts in the previous version)
\item smooth the counts (loess with span=0.1)
\end{enumerate}

\item Detect peak and putative dip/additional peak
\begin{itemize}
\item Detect the peak of the autoc which is its max. If the max is the first bin then the peak is defined as the first value for which n-1<n<n+1
\item Putative dip detection
\begin{itemize}
\item Dip detection from the first peak of the autoc
\item Putative dip if the min of the first 100ms is a local minimum and below 75% (50% in V1) of the fullrange of the autoc
\end{itemize}

\item If dip detected, detection of a 2nd peak:
\begin{itemize}
\item Latency from dip to 2nd dip must be of at least 12ms
\item Max of this
\end{itemize}
\end{itemize}


\item Fit from the 1st peak till the end of the autoc:
Use of the nlsLM function from minpack.lm package (formula= counts ~ A * exp(-mids/TAU) +B)
Starting values are sampled in the plausible solution space:
\begin{itemize}
\item A: [0,2*range(AC)]
\item B: [0,2*min(AC)]
\item TAU: [0 1000]
\end{itemize}
Note that previously starting values were set to the linearized solutions of the exponential equation, this resulted in some singular gradient matrix errors, probably due to limitations in the solution space due to the starting values.
\item If a putative dip was detected:
\begin{itemize}
\item Fit from 1st peak to dip
\item Fit from 2nd peak to dip
\end{itemize}

\item Post-hoc treatment of fits
A fit is validated if :
\begin{itemize}
\item a, B and TAU >0
\end{itemize}
If all the fits are validated, the units is marked as dip if the cumulative RMSE of dips fits is inferior to the mono RMSE. Otherwise it is marked as mono.
To be marked as dip, the two dip fits should be validated.
If the mono fit is validated and the unit is not marked as dip, then it is marked as mono.
If fits are not validated the unit is marked as non-fitted.

\end{enumerate}

We extracted TAU from units fitted with a single exponential. For this study and regarding their small number, we put dip units aside.


## Compute the autoc
We define a function that returns the ISI and the Autoc for a given timeserie
```{r,message=FALSE, warning=FALSE}

rm(spiketimeserie,curr_diff,autocorr,lg,dum,curr_autoc,unit,my_list,fr)

compute_autoc <- function(spiketimeserie,unit){

## Extract the ISIs
  curr_ISI <- data.frame();
  for (lg in 1:100){
    curr_diff <- diff(spiketimeserie, lag=lg);
    
    curr_diff <- curr_diff[curr_diff>0];
    curr_diff <- curr_diff[curr_diff<1];
    curr_diff <- curr_diff*1000;
    
    dum <- hist(curr_diff,breaks=seq(from=0,to=1000,length.out =300), plot=FALSE)
    dum <- data.frame(Unit=unit,"mids"=dum$mids, "counts"=dum$counts,"lag"=lg)
    curr_ISI <- rbind(dum,curr_ISI)
  }

  # Sum the ISIs to obtain the autoc
  curr_autoc <- curr_ISI %>% group_by(mids) %>%
    dplyr::summarise(Unit=unique(Unit),counts=sum(counts)) %>% ungroup();
  
  
  # Compute the probability density function of the autocorrelogram defined as :
# Pi = Xi/S(X) * 1/deltaT
# with Pi the PDF value for bin i, Xi the count of bin i, S(X) the sum of counts and deltaT the temporal bin width
  
  curr_autoc_tot <- sum(curr_autoc$counts);
  curr_binwidth <- mean(diff(curr_autoc$mids));
  
  curr_autoc <- curr_autoc %>% dplyr::mutate(countsn=counts/!!curr_autoc_tot * 1/!!curr_binwidth);
  
# smooth the autoc
 
  curr_autoc$countsnsm <- fitted(loess(countsn~mids, data=curr_autoc, span = 0.1));
  
    # Remove the first 10ms of the autoc to avoid noise contamination
  curr_autoc <- dplyr::filter(curr_autoc,mids >=10);
  
  # Pdf again
  curr_autoc_tot <- sum(curr_autoc$countsnsm);
  curr_autoc <- curr_autoc %>% dplyr::mutate(countsnsmn=countsnsm/!!curr_autoc_tot * 1/!!curr_binwidth);
  
  
 # Only keep ISIs up to the 5th lagged difference
 curr_ISI <- curr_ISI %>% dplyr::filter(lag<=5);
  
  # pdf of the ISIs
  curr_isi_tot <- sum(curr_ISI$counts);
  curr_ISI <- curr_ISI %>% dplyr::mutate(countsn=counts/!!curr_isi_tot * 1/!!curr_binwidth);
 
 
  # Store the autocorrelograms
  my_list <- list();
  my_list$Autoc<- curr_autoc;
  my_list$ISI<- curr_ISI;
  
  #Firing rate during the spike timeseries in 1ms bins 
  
  fr <- hist(spiketimeserie,breaks=seq(from=min(spiketimeserie,na.rm = T),to=max(spiketimeserie,na.rm = T)+1,by=1), plot=FALSE,include.lowest = T)
  
  fr <- data.frame("mids"=fr$mids, "counts"=fr$counts)
  #Remove last bin cause it was ceiled
  fr <- fr %>% dplyr::filter(mids != max(mids));
  
  #put bin from 0 to n sec 
  fr$mids <- fr$mids-min(fr$mids);
  
  fr$unit <- unit;

  my_list$fr<- fr;

  return(my_list)   
}
  
```

## Detect peak and putative dip/additional peak

 Wrap this procedure in a local function to compute it on units autocs

```{r,message=FALSE, warning=FALSE}
rm(autoc,peek.first,df.peak,unit,curr_min,curr_max,curr_dip,curr_local_min,curr_begin,curr_end,curr_post_peak)

dip_fun <- function(autoc,unit){
  
  # # Remove bins smaller than 10ms
  # autoc <- dplyr::filter(autoc,mids >=10);
  
  # Peak detection
  peek.first <- which(autoc$counts==max(autoc$counts));
  if(peek.first==1){
    peek.first <- 2;
    
    while(!(autoc$counts[peek.first-1] < autoc$counts[peek.first] & autoc$counts[peek.first+1] < autoc$counts[peek.first])){
      peek.first<-peek.first+1;
    }
    
  }
  df.peak <- data.frame(peak.first.lat=autoc$mids[peek.first],peak.first=autoc$counts[peek.first],unit);
  
  # Putative dip detection from first peak
  curr_post_peak <- dplyr::filter(autoc,mids >= df.peak$peak.first.lat);
  
  ## Global (post peak) min and max
  curr_min <- min(curr_post_peak$counts);
  curr_max <- max(curr_post_peak$counts);
  
  ## Minimum of the AC from first peak until 100ms after it
  curr_begin <- dplyr::filter(curr_post_peak,between(mids,df.peak$peak.first.lat,df.peak$peak.first.lat+100));
  
  curr_local_min <- min(curr_begin$counts);
  curr_local_min_lat <- curr_begin$mids[curr_begin$counts == curr_local_min];
  
  #This minimum should be a local minimum (that is AC(tmin-1)>AC(tmin)<AC(tmin+1))
  curr_local_min_idx <- which(curr_post_peak$mids == curr_local_min_lat);
  
  if((curr_post_peak$counts[curr_local_min_idx-1] > curr_post_peak$counts[curr_local_min_idx] & curr_post_peak$counts[curr_local_min_idx+1] > curr_post_peak$counts[curr_local_min_idx]) ==T){
    ## Is the local min inferior to 75% of the global range (50% in V1)?
    curr_dip <- curr_max-curr_local_min >=  3/4 * (curr_max -curr_min);
  }else{
    curr_dip <- F;
  }
  
  # Second peak detection if putative dip
  
  if(curr_dip==TRUE){
    df.peak$dip.lat <- curr_begin$mids[curr_begin$counts==curr_local_min];
    df.peak$dip <- curr_local_min;

    curr_end <- dplyr::filter(autoc,mids>=df.peak$dip.lat+12);
      
    df.peak$peak.second.lat <- curr_end$mids[curr_end$counts==max(curr_end$counts)];
    df.peak$peak.second  <- max(curr_end$counts);
    
  }else{
    df.peak$dip.lat <- NA;
    df.peak$dip <- NA;
    df.peak$peak.second.lat <- NA;
    df.peak$peak.second <- NA;
  }
  
  return(df.peak)
}
```

## Fit function

We wrap in a local function the fit. This functions requires the smoothed autocorrelogram and returns the Autoc metrics.


### Exponential function
```{r,message=FALSE, warning=FALSE}
rm(x,a,b,TAU)
expfit <- function(x,a,b,TAU) a*(exp(-x/TAU))+b;

rm(A,B,TAU,TAU.se,TAU.p,mf,st,autocsm,unit,fit.err)


exp_fit_fun <- function(autocsm,unit,st){

  #Fit exponential and extract Tau
  tryCatch({
    
    
    mf <- nlsLM(formula= counts ~ A * exp(-mids/TAU) +B , data=autocsm, start=st);
    
    
    A <- summary(mf)$coefficients[2,1];
    B <- summary(mf)$coefficients[3,1];
    TAU <- summary(mf)$coefficients[1,1];
    TAU.se <- summary(mf)$coefficients[1,2];
    TAU.p <- summary(mf)$coefficients[1,4];
    
    
    autocsm <- autocsm %>% dplyr::mutate(fitval=expfit(mids,A,B,TAU));
    
    fit.err <- sqrt(sum((autocsm$fitval-autocsm$counts)^2));
    
    
    # Return the fit 
  
  df.fit <- data.frame("unit"=unit,"TAU" = TAU,"TAU.se"=TAU.se,"TAU.p"=TAU.p,"A" = A,"B"=B,"fit.err"=fit.err);
    
  }, error=function(e){
    
    df.fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA);
  })
  
  if(!exists("df.fit")){
    df.fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA);
  }
  
  
  return(df.fit)
}

```


### Define starting values for the fits
50 Random selection in selected range, we keep the same starting values for all the units
A [0 ; 2*range(Autoc)]
B [0 ; 2*min(Autoc)]
TAU [0 ; 1000]
```{r}
rm(autoc,st)

f.start_val <- function(autoc){
  st <- data.frame(TAU=runif(50, min = 0, max = 1000),
                        A=runif(50, min = 0, max = 2*diff(range(autoc$counts))),
                        B=runif(50, min = 0, max = 2*min(autoc$counts)))
  
  return(st)
}

```


### Exponential fits on subparts of the autoc 

```{r}
rm(list = ls(pattern = "^curr*"))
rm(p,st,autocsm,unit)


exp_fit_autoc <- function(autocsm,unit){
  curr_peak <- dip_fun(autocsm,unit);
  
  # Fit from the 1st peak till the end
  curr_mono <- dplyr::filter(autocsm,mids>=curr_peak$peak.first.lat);
  
  # Repeat the fit 50 times with random starting values
  st <- f.start_val(curr_mono);
  
  curr_mono_fit <- data.frame();
  for (ifit in 1:nrow(st)){
    curr_st <- st[ifit,];
    curr_mono_fit <- rbind(curr_mono_fit,exp_fit_fun(curr_mono,unit,curr_st));
  }
  # # Keep the best fit
  # curr_mono_fit <- curr_mono_fit[which.min(curr_mono_fit$fit.err),];
  # 
  # if(nrow(curr_mono_fit)==0){ #if all iteration did not converge
  #   curr_mono_fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA,fit="mono");
  # }
  
  curr_mono_fit$dof <- nrow(curr_mono);
  curr_mono_fit$fit <- "mono";
  
  curr_fit <- curr_mono_fit;
  
  ## Fit from 1st peak to dip
  if(!is.na(curr_peak$dip)){
    # Fit from the 1st peak till the end
    curr_first_peak_to_dip <- dplyr::filter(autocsm,between(mids,curr_peak$peak.first.lat,curr_peak$dip.lat));
    
    # Repeat the fit 50 times with random starting values
  st <- f.start_val(curr_first_peak_to_dip);
    
    curr_first_peak_to_dip_fit <- data.frame();
    for (ifit in 1:nrow(st)){
      curr_st <- st[ifit,];
      curr_first_peak_to_dip_fit <- rbind(curr_first_peak_to_dip_fit,exp_fit_fun(curr_first_peak_to_dip,unit,curr_st));
    }
    # # Keep the best fit
    # curr_first_peak_to_dip_fit <- curr_first_peak_to_dip_fit[which.min(curr_first_peak_to_dip_fit$fit.err),];
    # 
    # if(nrow(curr_first_peak_to_dip_fit)==0){ #if all iteration did not converge
    #   curr_first_peak_to_dip_fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA,fit="dip1");
    # }
    curr_first_peak_to_dip_fit$dof <- nrow(curr_first_peak_to_dip);
    curr_first_peak_to_dip_fit$fit <- "dip1";
  }else{
    curr_first_peak_to_dip_fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA,fit="dip1",dof=NA);
  }
  
  
  curr_fit <- rbind(curr_fit,curr_first_peak_to_dip_fit);
  
  
  ## Fit from 2nd peak to end
  if(!is.na(curr_peak$peak.second)){
    # Fit from the 1st peak till the end
    curr_second_peak_to_end <- dplyr::filter(autocsm,mids>=curr_peak$peak.second.lat);
    
    # Repeat the fit 50 times with random starting values
    st <- f.start_val(curr_second_peak_to_end);
    
    curr_second_peak_to_end_fit <- data.frame();
    for (ifit in 1:nrow(st)){
      curr_st <- st[ifit,];
      curr_second_peak_to_end_fit <- rbind(curr_second_peak_to_end_fit,exp_fit_fun(curr_second_peak_to_end,unit,curr_st));
    }
    # # Keep the best fit
    # curr_second_peak_to_end_fit <- curr_second_peak_to_end_fit[which.min(curr_second_peak_to_end_fit$fit.err),];
    # 
    # if(nrow(curr_second_peak_to_end_fit)==0){ #if all iteration did not converge
    #   curr_second_peak_to_end_fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA,fit="dip2");
    # }
    curr_second_peak_to_end_fit$fit <- "dip2";
    curr_second_peak_to_end_fit$dof <- nrow(curr_second_peak_to_end);
  }else{
    curr_second_peak_to_end_fit <- data.frame("unit"=unit,"TAU" = NA,"TAU.se"=NA,"TAU.p"=NA,"A" = NA,"B"=NA,"fit.err"=NA,fit="dip2",dof=NA);
  }
  
  
  
  curr_fit <- rbind(curr_fit,curr_second_peak_to_end_fit);
  
  
  curr_fit$peak.lat <- curr_peak$peak.first.lat;
  curr_fit$peak <- curr_peak$peak.first;
  curr_fit$dip.lat <- curr_peak$dip.lat;
  curr_fit$dip <- curr_peak$dip;
  curr_fit$peak.second.lat <- curr_peak$peak.second.lat;
  curr_fit$peak.second <- curr_peak$peak.second;
  
  return(curr_fit)
}
```

## Fit every unit

```{r,message=F,warning=F}
ls.unit <- neurons_info$unit;
Nunit <- length(ls.unit);

rm(list = ls(pattern = "^ISI*"))
rm(list = ls(pattern = "^Autoc*"))

Autoc <- data.frame();
ISI <- data.frame();
Autoc.metrics <- data.frame();

for (iunit in 1:Nunit){
  
  rm(list = ls(pattern = "^curr*"))
  rm(spike01,unit,filename,p)
  
  fprintf("Unit %i out of %i\n",iunit,Nunit)
  
  
  filename <- paste0(path,"Spike timeseries/",iunit,".txt")
  unit <- read.csv(filename, header = FALSE)
  
  spike01 <- unit$V1/100000  # transform back in seconds
  
  curr_ls <- compute_autoc(spike01,iunit);
  curr_autoc <- curr_ls$Autoc;
  
  Autoc <- rbind(Autoc,curr_autoc);
  ISI <- rbind(ISI,curr_ls$ISI);
  
  #Use the smoothed normalized value for the analysis (it needs to be named 'counts')
  curr_autoc$counts <- curr_autoc$countsnsmn;
  
  Autoc.metrics <- rbind(Autoc.metrics,exp_fit_autoc(curr_autoc,iunit));
}

```

## Post-hoc treatment of fits

```{r}
# Post-hoc treatment
# Only keep valid fits. In order to be valid, A, B and TAU must be positive (like Matthieu)
Autoc.metrics <- Autoc.metrics %>% dplyr::filter(TAU<1000 & TAU>0 & A >0 & B>0);

# For a given unit x fit type (mono,dip,etc), only keep the fit with the minimal error
Autoc.metrics <- Autoc.metrics %>% group_by(unit,fit) %>%
  slice(which.min(fit.err)) %>% ungroup();

# If all fits are valid for a given unit, compare dip and mono errors ; if the mono is valid but the other two, validate the mono, if none: none
## When comparing the units, use RMSE that is the square root of the sum of the errors divided by the sum of the dof 

dum <- Autoc.metrics %>% tidyr::pivot_wider(id_cols=unit,names_from=fit,values_from = c(fit.err,dof)) %>% group_by(unit) %>%
  dplyr::mutate(rmse_mono=fit.err_mono/sqrt(dof_mono),rmse_dip=sqrt(sum((c(fit.err_dip1,fit.err_dip2)^2))/(dof_dip1+dof_dip2))) %>%
  dplyr::mutate(bestfit=case_when(
    rmse_dip < rmse_mono ~"dip",
    !is.na(rmse_dip) & is.na(rmse_mono) ~"dip",
    !is.na(rmse_mono) ~"mono",
    TRUE ~"none"))
                                                                                                     
Autoc.metrics <- dum %>% dplyr::select(unit,bestfit) %>%
  left_join(Autoc.metrics,.,by="unit");

Autoc.mono <- Autoc.metrics %>% dplyr::filter(fit=="mono" & bestfit =="mono");

Autoc.mono <- neurons_info %>% dplyr::select(unit,area,SWg) %>%
  left_join(Autoc.mono,.,by="unit",copy=F);

Autoc.mono$area <- factor(Autoc.mono$area);
levels(Autoc.mono$area) <- list(MCC="1",LPFC="2");

Autoc.dip <- Autoc.metrics %>% dplyr::filter(is.element(fit,c("dip1","dip2")) & bestfit =="dip");

ls.dip <- unique(Autoc.dip$unit);
ls.mono <- unique(Autoc.mono$unit);

ls.none <- setdiff(unique(neurons_info$unit),c(ls.dip,ls.mono));


Ndip <- length(ls.dip);
Nmono <- length(ls.mono);
Nnone <- length(ls.none);
Nunits <- length(unique(neurons_info$unit));
```


# Spike autocorrelogram features between MCC and LPFC (figure 2c)


## TAU distribution

MCC units have longer TAU

```{r,fig.cap="true autoc TAU distributions"}

Autoc.area.med <- Autoc.mono %>% group_by(area) %>%
  dplyr::summarise(TAU=median(TAU,na.rm=T),peak.lat=median(peak.lat,na.rm=T),n=n()) %>% ungroup();

Autoc.area.sw.med <- Autoc.mono %>% group_by(area,SWg) %>%
  dplyr::summarise(TAU.sd=sd(TAU,na.rm=T),TAU=median(TAU,na.rm=T),lat.sd=sd(peak.lat,na.rm=T),peak.lat=median(peak.lat,na.rm=T),n=n()) %>% ungroup();

my.binsize <- 25;

df <- Autoc.mono %>% 
  dplyr::mutate(my.bin = fun.bin(TAU,my.binsize)) %>%
    group_by(area,SWg,my.bin) %>% 
    dplyr::summarise (n = n()) %>%
    dplyr::mutate(freq = n / sum(n)) %>% ungroup();


ggplot(df,aes(x=my.bin,y=freq,fill=area))+
  geom_vline(data=Autoc.area.sw.med,aes(xintercept=TAU,colour=area), size=1, linetype='dashed',show.legend = F)+
  geom_bar(stat="identity",position="identity",alpha=.6,)+
  scale_fill_manual(values=colArea)+
  scale_colour_manual(values=colArea)+
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),aspect.ratio = 1,legend.position="none")+ 
  labs(x="TAU",y="Frequency")+
  facet_grid(cols=vars(SWg),labeller = labeller(SWg=c(FS="Fast Spiking",RS="Regular Spiking")))


# TAU and peak latency BLOM transformation for stat test
N <- nrow(Autoc.mono);
Autoc.mono <- Autoc.mono %>% dplyr::mutate(tau.rank = rank(TAU),lat.rank=rank(peak.lat)) %>% ungroup() %>% group_by(unit) %>%
  dplyr::mutate(tau.blom=blom(tau.rank,!!N),lat.blom=blom(lat.rank,!!N)) %>% ungroup();


shapiro.test(Autoc.mono$tau.blom)
summary(lm(tau.blom ~ area*SWg, data=Autoc.mono))

```


## Autocorr peak latency distribution

```{r,fig.cap="true autoc peak latency distributions"}
my.binsize <- 7.5;

df <- Autoc.mono %>% 
  dplyr::mutate(my.bin = fun.bin(peak.lat,my.binsize)) %>%
    group_by(area,SWg,my.bin) %>% 
    dplyr::summarise (n = n()) %>%
    dplyr::mutate(freq = n / sum(n)) %>% ungroup();


ggplot(df,aes(x=my.bin,y=freq,fill=area))+
  geom_bar(stat="identity",position="identity",alpha=.6)+
  geom_vline(data=Autoc.area.sw.med,aes(xintercept=peak.lat,colour=area), size=1, linetype='dashed',show.legend = F)+
  scale_fill_manual(values=colArea)+
  scale_colour_manual(values=colArea)+
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),aspect.ratio = 1,legend.position = "none")+ 
  xlim(0,300)+
  labs(x="Autocorr peak latency (ms)",y="Frequency")+
  facet_grid(.~SWg)


shapiro.test(Autoc.mono$lat.blom)
summary(lm(lat.blom ~ area*SWg, data=Autoc.mono))

```

TAU and autocorr peak latency medians for each group of units.

```{r,echo=F,results="asis"}
kable(Autoc.area.med,caption="Spike-rate autoc area medians")

kable(Autoc.area.sw.med,caption="Spike-rate autoc area/spike shape medians")
```






