---
title: "Plotting a isovoltage topoplots. This is adapted from some code that Matt Craddock posted on his blog"
output: html_notebook
fontsize: 11pt
geometry: margin=1in
---

##Loading and formatting the data

Load some libraries.

```{r}
library(tidyverse)
library(akima) #for interpolation
library(scales)
library(mgcv) #another interpolation option
library(gridExtra)
library(png)
library(grid)
```

Load data from individual ERPs. One could load directly from a grandaveraged ERP, but since we used individual ERPs for the waveforms, we can use them again here.

* Set paths and load data

```{r set_path_to_data_and_load, results = 'hide', warning=FALSE, message=FALSE, error=FALSE}
data.base.path <- "/Users/jsprouse/Desktop/"
erp.path <- "exported ERPs"
erp.files <- file.path(data.base.path, erp.path) %>% 
  dir(full.names = TRUE)
erp.labels <- file.path(data.base.path, erp.path) %>% 
  dir() %>%
  stringr::str_replace(".txt", "")

## Load the data
erp.data.list <- map(erp.files, read_tsv)
```

Preprocessing
=============

* load the data and preprocess it (to long format)

```{r preprocess_data}
## -------------------
## Preprocess the data
## -------------------
## 1. Bind the individual subject data
## 2. Create Subject and Condition labels
## 3. Restructure the data to long format
## 4. Filter out time points that won't be used

erp.data <- bind_rows(erp.data.list) %>%
  dplyr::mutate(Subject_Cond = rep(erp.labels, 
                                   each = 32)) %>% #this is the number of channels
  tidyr::separate(col = Subject_Cond, 
           into = c("subject", "condition"), sep = "_") %>%
  tidyr::gather(key = timet, value = mV, 
         -time, -subject, -condition) %>%
  dplyr::filter(!is.na(mV)) %>%
  dplyr::rename(channel = time, time.ms = timet) %>%
  dplyr::mutate(mV = as.double(mV)) %>%
  dplyr::mutate(time.ms = as.double(time.ms)) %>%
  arrange(condition) %>%
  dplyr::mutate(channel = as.factor(channel), subject = as.factor(subject), condition = as.factor(condition))
```

##Electrode locations

We need electrode coordinates (Cartesian) too. Here is a giant list, formed into a data frame, taken from Matt's code. These are not the same as the coordinates we use in the topographic waveform plot (taken from the erpR package).

```{r}
elec.coord = data.frame(channel=c("Fp1",	"AF7",	"AF3",	"F1",	"F3",	"F5",	"F7",	"FT7",	"FC5",	"FC3",	"FC1",	"C1",	"C3",	"C5",	"T7",	"TP7",	"CP5",	"CP3",	"CP1",	"P1",	"P3",	"P5",	"P7",	"P9",	"PO7",	"PO3",	"O1",	"Iz",	"Oz",	"POz",	"Pz",	"CPz",	"Fpz",	"Fp2",	"AF8",	"AF4",	"AFz",	"Fz",	"F2",	"F4",	"F6",	"F8",	"FT8",	"FC6",	"FC4",	"FC2",	"FCz",	"Cz",	"C2",	"C4",	"C6",	"T8",	"TP8",	"CP6",	"CP4",	"CP2",	"P2",	"P4",	"P6",	"P8",	"P10",	"PO8",	"PO4",	"O2",	"LO1",	"LO2",	"SO1",	"IO1",	"M1",	"M2", "FT9", "FT10"), x=c(-0.158507951196727,	-0.306220792940307,	-0.160894602703784,	-0.111231689343979,	-0.221253907813213,	-0.327718412610652,	-0.426753882970127,	-0.505736594237533,	-0.381956848927333,	-0.25547983525988,	-0.127945796584781,	-0.13319,	-0.26669,	-0.3999,	-0.53318,	-0.505742347629983,	-0.381951819153741,	-0.255491477354236,	-0.127956990743876,	-0.111218290253606,	-0.221267739734655,	-0.327693903532661,	-0.426737596055039,	-0.529417172201933,	-0.306206021969086,	-0.160887810154621,	-0.158542158546356,	7.75813747259848E-17,	6.20516286659972E-17,	4.65292304868045E-17,	3.10301005967956E-17,	1.55064777708038E-17,	0,	0.158507951196727,	0.305957920131658,	0.160894602703784,	0,	0,	0.111132030087498,	0.220965559664675,	0.327483020901138,	0.426495945228168,	0.505736594237533,	0.381956848927333,	0.25547983525988,	0.127945796584781,	0,	0,	0.13348,	0.26667,	0.3999,	0.53318,	0.505559973212276,	0.381951819153741,	0.255491477354236,	0.127956990743876,	0.111118643002187,	0.220979398075575,	0.327458529427387,	0.426512248599182,	0.529417172201933,	0.305972696504555,	0.160887810154621,	0.158542158546356,	-0.438655260304613,	0.438655260304613,	-0.246525910621202,	-0.315269162639131,	-0.735005330945149,	0.735005330945149, -0.735005330945149,	0.735005330945149), y=c(0.489989723879405,	0.423151810667294,	0.389182995158822,	0.255900863979949,	0.264176032935205,	0.280839046672563,	0.311045993013766,	0.164831379442813,	0.144085733011642,	0.133419721838956,	0.128281197126804,	8.1555353589218E-18,	1.63300527432304E-17,	2.44868127489513E-17,	3.26478590184693E-17,	-0.164813725804961,	-0.144099065733086,	-0.133397426509471,	-0.128270031261287,	-0.255906687706798,	-0.264164447745561,	-0.280867644429768,	-0.311068337368441,	-0.39447167981708,	-0.423162499531643,	-0.389185803240111,	-0.489978656742783,	-0.6335,	-0.50669,	-0.37994,	-0.25338,	-0.12662,	0.50669,	0.489989723879405,	0.423317241804193,	0.389182995158822,	0.37994,	0.25338,	0.255671586783967,	0.264299964892311,	0.280637327206245,	0.311382616091412,	0.164831379442813,	0.144085733011642,	0.133419721838956,	0.128281197126804,	0.12662,	0,	8.17329273750943E-18,	1.63288280964312E-17,	2.44868127489513E-17,	3.26478590184693E-17,	-0.165340139063697,	-0.144099065733086,	-0.133397426509471,	-0.128270031261287,	-0.255677405292983,	-0.264288394800371,	-0.280665904422416,	-0.311360284389113,	-0.39447167981708,	-0.423306561600138,	-0.389185803240111,	-0.489978656742783,	0.487176021789961,	0.487176021789961,	0.528676521412089,	0.61875057065737,	-0.135163946310443,	-0.135163946310443, 0.135163946310443,	0.135163946310443))
            
```

Find the electrodes in our dataset:
```{r}
coords = elec.coord %>% 
  filter(channel %in% levels(erp.data$channel)) %>%
  droplevels()
```

Combine the electrode locations to the data using left_join():

```{r}

#first we make sure that the order of the levels of each dataframe are identical. This eliminates any warnings or errors.
erp.data$channel=factor(erp.data$channel, levels=sort(levels(erp.data$channel)))

coords$channel=factor(coords$channel, levels=sort(levels(coords$channel)))

#now we combine using left_join
full.locs <- erp.data %>% left_join(coords, by = "channel")
```

This assumes that we have three bins: the two conditions, and a bin that is the difference between them. If you don't have a difference bin, you need to calcluate a difference. There is code a bit further down to do that. But we can skip it here.
```{r}
all = filter(full.locs, condition == "congruent" | condition == "incongruent" | condition == "N400")
all = droplevels(all)
```

Choose a subset of time that includes the words we want
```{r}
plottable = all %>%
  filter(time.ms>= -200 & time.ms<=1000)
```

Do you need to re-baseline? If so, here is some code to do it right here in R.
```{r}
#calculate the baseline value
# baseline =  plottable %>%
#   filter(time.ms >= -200 & time.ms<=0) %>%
#   group_by(subject, channel, condition) %>%
#   summarise(microvolts = mean(mV)) %>%
#   arrange(subject,channel, condition) %>%
#   rename(baseline = microvolts)
# 
# #left_join the baseline value to the original dataset
# plottable = left_join(plottable, baseline)
# 
# #create the re-baselined dataset
# plottable = plottable %>%
#   mutate(mV.original = mV, mV = mV.original - baseline)

```

Choose a subset to work on. In this case, we only have one subset, so this seems like a waste. But if we had more than one subset, we could choose then one we want here.
```{r}
plotting = plottable
```


Here is some code to calculate difference effects. If you already have this in your data set, you don't need this. 

```{r}
# #split by conditions
# split.conditions = plotting %>%
#   split(f=plotting$condition)
# 
# #calculate differences in voltage
# N400.difference = split.conditions$incongruent$mV - split.conditions$congruent$mV
# 
# #recreate the dataset
# N400.effect = data.frame(select(split.conditions$incongruent, -condition, -mV), condition=rep("difference", times=length(N400.difference)), mV = N400.difference)
# 
# #now bind it to the original data set
# plotting2 = rbind(plotting, select(N400.effect, channel, subject, condition, time.ms, mV, x, y))

```

Calculate mean values in the windows that we want to plot. This also labels the windows.

```{r}

startList = c(100, 200, 300)
endList = c(150, 250, 500)

for(i in 1:length(startList)){
  #select a time
  temp = plotting %>%
    filter(time.ms >= startList[i] & time.ms<=endList[i]) %>%
    group_by(condition, channel, y, x) %>%
    summarise(microvolts = mean(mV))
  
  if(startList[i] == 100 | startList[i] == 600 | startList[i] == 1100){
    temp2 = data.frame(window = factor(rep("N1", times=nrow(temp))), temp)}
  
  if(startList[i] == 200 | startList[i] == 700 | startList[i] == 1200){
    temp2 = data.frame(window = factor(rep("P2", times=nrow(temp))), temp)}
  
  if(startList[i] == 300 | startList[i] == 800 | startList[i] == 1300){
    temp2 = data.frame(window = factor(rep("300-500", times=nrow(temp))), temp)}
  
  
  if(i==1){
    plottable=temp2
  }else{
    plottable = rbind(plottable, temp2)
  }
}


```

##Draw the head, and choose some colors

Do we want to define our own colors?

```{r}
#this is matlab's jet palette
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

#this is the 9 color diverging palette from colorbrewer
#it comes out a little pale 
jon1=colorRampPalette(c("#2166ac", "#4393c3", "#92c5de", "#d1e5f0", "#f7f7f7", "#fddbc7", "#f4a582", "#d6604d", "#b2182b"))

#This is 13 color sequence that I made up by taking colors from two sequential palettes
jon2=colorRampPalette(c("#08306b","#08519c","#2171b5","#4292c6","#6baed6","#9ecae1", "#f7f7f7","#fc9272", "#fb6a4a", "#ef3b2c", "#cb181d","#a50f15","#67000d"))

#This is 7 of the 13 colors above
jon3=colorRampPalette(c("#08306b","#2171b5","#6baed6", "#f7f7f7", "#fb6a4a",  "#cb181d","#67000d"))

#We can also use colorbrewer palettes directly
```

Draw the head (taken directly from Matt's code):

```{r}
theme_topo <- function(base_size = 12)
  {
  theme_bw(base_size = base_size) %+replace%
      theme(
            rect             = element_blank(),
            line             = element_blank(),
            axis.text = element_blank(),
            axis.title = element_blank()
           )
}

circleFun <- function(center = c(0,0),diameter = 1, npoints = 100) {
  r = diameter / 2
  tt <- seq(0,2*pi,length.out = npoints)
  xx <- center[1] + r * cos(tt)
  yy <- center[2] + r * sin(tt)
  return(data.frame(x = xx, y = yy))
}

headShape <- circleFun(c(0, 0), round(max(coords$x)), npoints = 100) # 0
nose <- data.frame(x = c(-0.075,0,.075),y=c(.495,.575,.495))

ggplot(headShape,aes(x,y))+
  geom_path()+
  geom_text(data = coords,
            aes(x, y, label = channel))+
  geom_line(data = nose,
            aes(x, y))+
  theme_topo()+
  coord_equal()
```



##Interpolation

The trick here is that we want to interpolate based on each grouping factor (effect and time). I can't get dplyr pipelines to work here. I don't know why. So I am just going to loop through the dataframe and re-assemble it


```{r}
gridRes <- 150 # Specify the number of points for each grid dimension i.e. the resolution/smoothness of the interpolation

#loop through effects, and through times

#first, define a results dataframe. We are going to rbind() into this. 

results = data.frame()

for(i in 1:length(levels(plottable$condition))){
  for(j in 1:length(levels(plottable$window))){
   
    tempData = filter(plottable, condition == levels(plottable$condition)[i] & window == levels(plottable$window)[j])
   
    #need to check if there is data here, because the final word in the sentence doesn't have a 500-800 window
    if(nrow(tempData)>1){
    tmpTopo <- with(tempData,
                interp(x = x, y = y, z = microvolts,
                       xo = seq(min(x)*2,
                                max(x)*2,
                                length = gridRes),
                       yo = seq(min(y)*2,
                                max(y)*2,
                                length = gridRes),
                       linear = FALSE,
                       extrap = TRUE)
                ) 

  interpTopo <- data.frame(x = tmpTopo$x, tmpTopo$z)

  names(interpTopo)[1:length(tmpTopo$y)+1] <- tmpTopo$y

  interpTopo <- gather(interpTopo,
                    key = y,
                    value = microvolts,
                    -x,
                    convert = TRUE)

  interpTopo$incircle <- sqrt(interpTopo$x^2 + interpTopo$y^2) < .7 # mark grid elements that are outside of the plotting circle

  interpTopo <- interpTopo[interpTopo$incircle,] #remove the elements outside the circle

  tempInterpTopo = data.frame(condition = rep(as.character(levels(plottable$condition)[i]), times=nrow(interpTopo)), window=rep(as.character(levels(plottable$window)[j]), times=nrow(interpTopo)), interpTopo)
  
  #add to results file
  results = rbind(results, tempInterpTopo, row.names=NULL)
  } else {next} #move to next value of window loop if there is no data to work with
  
  }
}
```

Here we change some things for aesthetic reasons.

```{r}
#reorder the levels of the effect factor

results$condition = factor(results$condition, levels=c("congruent", "incongruent", "N400"))

results$window = factor(results$window, levels=c("N1", "P2", "300-500"))

#define a list of points for the channel locations in the plot
dots = filter(plottable, channel !="FT9" & channel !="FT10" & channel !="M1" & channel !="M2")

```


```{r}
maskRing <- circleFun(diameter = 1.42) #create a circle round the outside of the plotting area to mask the jagged edges of the interpolation

akimaPlot <- ggplot(results,
                    aes(x = x, y = y, fill = microvolts)
                    ) +
  facet_grid(condition ~ window) +
  geom_raster() +
  stat_contour(aes(z = microvolts),
               colour = "black",
               binwidth = 0.5,
               alpha=0.2) + #this sets the opacity/transparency of the lines in the plot
  #theme_topo()+
  theme_minimal()+
  theme(line=element_blank(), axis.text = element_blank(),axis.title = element_blank())+
  scale_fill_gradientn(colours = jon1(10),
                       limits = c(-2,2),
                       guide = "colourbar",
                       oob = squish) + 
  #theme(legend.position = "bottom") +
  #guides(fill = guide_colorbar(direction = "horizontal")) +
  geom_path(data = maskRing,
            aes(x, y, z = NULL, fill =NULL),
            colour = "white",
            size = 6)+
  geom_point(data = dots,
             aes(x, y),
             size = 1)+
  geom_path(data = headShape,
            aes(x, y, z = NULL, fill = NULL),
            size = 1.5)+
  geom_path(data = nose,
            aes(x, y, z = NULL, fill = NULL),
            size = 1.5)+
  coord_equal()

akimaPlot
```

plot it and save it

```{r}
quartz(width = 9, height = 9)
akimaPlot
ggsave("isovoltage.N400.pdf")
```

