---
title: "meaRtools: Basics of correlating spike trains"
author: "Stephen J Eglen"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{meaRtools: Basics of correlating spike trains}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r rsetup-knitr,eval=TRUE,include=FALSE}
require(knitr)
options(width=60)
opts_chunk$set(cache=TRUE)
```


# Correlating spike trains

We use the method proposed by (Cutts and Eglen, J Neurosci 2014) to
compute the correlation between a pair of spike trains.  This measure
is called the "Spike Time Tiling Coefficient" (sttc).  To illustrate
these ideas, we first generate three related synthetic spike trains.
A is a Poisson train at 1 Hz; B simply shifts each spike in A
uniformly by up to +/- 0.1 seconds.  Finally C is simply B shifted
rightwards by 0.8 seconds.

(These correlation routines assume that the spike trains are sorted, smallest
time first.)


```{r}
require(meaRtools)
```

```{r}
poisson_train <- function(n = 1000, rate = 1, beg = 0) {
  ## Generate a Poisson spike train with N spikes and firing rate RATE.
  ## BEG is time of start of recording
  ## Check that the histogram looks exponentially distributed.
  ## hist( diff(poisson_train()))
  x <- runif(n)
  isi <- log(x) / rate
  spikes <- beg - cumsum(isi)
  spikes
}

a = poisson_train()
b = sort(a + runif(length(a), -0.1, 0.1))
c = b + 0.8
allspikes = list(a=a, b=b, c=c)
```

To compute the correlation between A and B we need three extra bits of
information.  First, dt is the coincidence window (typically between
10 ms and 1 s); we then provide the beginning and end time in seconds
of the two recordings.

```{r}
range = range(unlist(allspikes))
beg = range[1]
end = range[2]
sttc(a, b, dt=0.01, rec_time = c(beg, end))
```

To compute all pairwise correlations we can put all the spikes into a
list.  A upper-diagonal matrix is then returned with all pairwise
correlations; the leading diagonal should be 1.0 for autocorrelations.

```{r}
sttc_allspikes1(allspikes, 0.01, beg, end)
```

Since the STTC method was published, we have extended it to produce a
correlogram-like profile. For a given value of tau, we cross correlate
the spike trains A and B by adding tau to the time of each spike in B.
This allows us to detect correlations when there are temporal delays.
For this calculation, we simply also need to state the maximum
absolute tau value (tau_max), and the step size for tau (tau_step) in
sttcp.  In the following plot, we compare our spike trains when using
a range of coincidence windows, dt:


```{r, fig.width = 6, fig.height=6}
par(mfrow=c(2,2), las=1)
dt = 1.0;  plot(sttcp(a, b, dt = dt, beg = beg, end = end), main=sprintf("a x b: dt = %.3f",dt))
dt = 0.5;  plot(sttcp(a, c, dt = dt, beg = beg, end = end), main=sprintf("a x c: dt = %.3f",dt))
dt = 0.1;  plot(sttcp(a, c, dt = dt, beg = beg, end = end), main=sprintf("a x c: dt = %.3f",dt))
dt = 0.05; plot(sttcp(a, c, dt = dt, beg = beg, end = end), main=sprintf("a x c: dt = %.3f",dt))
```

It clearly shows that a and c are ofset by about 0.8 seconds.  We also
see that for small values of dt, the cross-correlation peak drops
significantly below 1.0.

<!--This file can be rendered using rmarkdown::render("sttc.Rmd")   -->
