strike behavioral analysis, force analysis

strike.probe.r

## strike.probe.r
## a software for the analysis of strike data obtained from mantis shrimp
## https://gist.github.com/4180275

library(R.matlab)
library(lattice)
library(timsac)
library(ggplot2)

working.dir <- "~/Dropbox/Analysis2/R_programming/strike_probe"

## functions for high sampling rate (500k Hz) strike force data
## class: strike

read.strike <- function(strike.dir=NA, animal.id=NA, species=NA, force.scale=NA){
  ## read strike data into "strikes" in the current workspace
  ## ex) read.strike(strike.dir=N.oerstedii.4.dirs[2], animal.id="Y4",
  ##                     species="N.oerstedii", force.scale=15.03 (N/V))
  previous.working.dir <- getwd()
  
  if (is.na(strike.dir)){
    cat("Enter the directory path containing force data. \n")
    strike.dir <- scan( what="" )
  }
  setwd(strike.dir)
  
  files <- list.files()
  files <- files[order(as.numeric(gsub('.mat', '', gsub('strike', '', files))))]
  
  if (is.na(animal.id)){
    cat("Enter animal id and push return key. \n")
    animal.id <- scan( what="" )
  }
  if (is.na(species)){
    cat("Enter animal species and push return key. \n")
    species <- scan( what="" )
  }
  if (is.na(force.scale)){
    cat("Enter Newton/Voltage and push return key. \n")
    force.scale <- as.numeric(scan( what="" ))
  }
  strikes <- list()
  n = 1

  for (each.file in files) {
    cat(".")
    one.strike <- readMat(each.file)
    
    is.rewarded <- one.strike[[1]][[1]][1,1]
    strike.time <- one.strike[[1]][[2]][1,1]
    voltage <- one.strike[[1]][[3]][,1]
    offset <- one.strike[[1]][[4]][1,1]
    voltage <- voltage + offset
    force <- ( voltage + offset ) * force.scale
    
    if ( n == 1 ) {
      duration <- one.strike[[1]][[5]][1,1]
      strike.threshold <- one.strike[[1]][[6]][1,1]
      reward.threshold <- one.strike[[1]][[7]][1,1]
      time <- seq(0, by = 1/500, length.out = length(force))
      strikes <- list(header = list(
                         animal.id = animal.id,
                         species = species,
                         strike.dir = strike.dir,
                         force.scale = force.scale,
                         duration = duration,
                         strike.threshold = strike.threshold,
                         reward.threshold = reward.threshold,
                         time = time ))
    }
    strikes <- c( strikes, list(each.strike = list(
                                  is.rewarded = is.rewarded,
                                  strike.time = strike.time,
                                  voltage = voltage,
                                  force = force )))
    n = n + 1
  }
  cat("\n")
  setwd(previous.working.dir)
  class(strikes) <- "strikes"
  return(strikes)
}

select.strikes <- function(strikes){
  num <- length(strikes)
  act <- c()
  time <- strikes[[1]]$time
  for (i in 2:num) {
    x <- strikes[[i]]
    # browser()
    plot(time, x$force, type="l")
    midline <- diff(range(time))/2
    text(max(time), min(x$force),
         labels="click here, \n if this is an actual strike", pos=2)
    abline(v=midline, col="red")
    
    is.act <- locator(1)$x
    if(is.act > midline) {
      cat("strike", paste(i-1), '\n')
      act <- c(act, i)
    } else {
      cat("\n")
    }
  }
  new.strikes <- list()
  new.strikes$header <- strikes$header
  for (i in act) {
    new.strikes <- c(new.strikes, each.strike=list(strikes[[i]]))
  }
  class(new.strikes) <- "strikes"
  return(new.strikes)
}

print.strikes <- function(strikes){
  strike.num <- length(strikes) - 1
  cat(paste(strike.num, 'strikes data. The header information is shown below.\n'))
  cat(' ', 'animal id: ')
  cat(strikes$header$animal.id, '\n')
  cat(' ', 'species: ')
  cat(strikes$header$species, '\n')
  cat(' ', 'directory of the data: ')
  cat(strikes$header$strike.dir, '\n')
  cat(' ', 'force scale: ')
  cat(strikes$header$force.scale, '\n')
  cat(' ', 'duration: ')
  cat(strikes$header$duration, '\n')
  cat(' ', 'strike threshold: ')
  cat(strikes$header$strike.threshold, '\n')
  cat(' ', 'reward.threshold: ')
  cat(strikes$header$reward.threshold, '\n')
}

export.as.text <- function(strikes){
  strike.n <- length(strikes) - 1
  force.df <- data.frame()
  mkdirs("exported_strikes")
  setwd("./exported_strikes")
  for ( i in 1:strike.n ){
    strike.ind <- i + 1
    voltage <- strikes[[strike.ind]]$voltage
    offset <- mean(strikes[[strike.ind]]$voltage[1:(length(voltage)/4)])
    voltage <- voltage - offset
    duration <- strikes$header$duration
    voltage.df <- data.frame(time = seq(duration/length(voltage),
                               duration, duration/length(voltage)),
                           voltage=voltage, voltage=voltage, voltage=voltage)
    write.table(voltage.df, file=paste("strike", i, ".txt", sep=""),
                col.names=FALSE,
                row.names=FALSE)
  }
  setwd("../")
}


plot.strikes <- function(my.strikes, which.strikes=NA,
                         cummulative=FALSE) {
  time <- my.strikes$header$time
  if (is.na(which.strikes)) {
    which.strikes <- 1:(length(my.strikes)-1) + 1
  }
  strike.num <- length(my.strikes) - 1
  d <- data.frame()
  time.stamps.rewarded <- c()
  rewarded.strike.ind <- c()
  for (i in which.strikes) {
    d <- rbind(d, data.frame(
                    strike.n = i,
                    force = my.strikes[[i]]$force,
                    time = time))
    if (my.strikes[[i]]$is.rewarded) {
      time.stamps.rewarded <- c(time.stamps.rewarded, my.strikes[[i]]$strike.time)
      rewarded.strike.ind <- c(rewarded.strike.ind, i - 1) # because of header
    } 
  }

  ## col.seq <- rep("black", strike.num)
  ## col.seq[rewarded.strike.ind] <- rep("red", length(rewarded.strike.ind))
  
  if (cummulative==TRUE) {
    plot(time.stamps.rewarded, 1:length(time.stamps.rewarded), type="s",
         xlim=c(0, 2000), ylim=c(0, 50),
         xlab="Time (sec)", ylab="Cummulative # of actions")
  } else {
    xyplot(force~time|strike.n,
           d, col="black", type="l",
           xlab="Time (ms)", ylab = "Force (N)", strip=FALSE)
  }
}

## make data frame of animal 04
make.df04 <- function(){
  strikes <- read.strike(strike.dir=readline(prompt="enter dir: "),
                         animal.id="animal04", species="unknown", force.scale=15.03)
  new.strikes <- select.strikes(strikes)
  d <- summary(new.strikes)
  session <- as.numeric(readline(prompt="enter session number: "))
  condition <- readline(prompt="enter condition (rewarded or nonrewarded?): ")
  is.rewarded <- rep(0, length(new.strikes)-1)
  is.rewarded[d$rewarded.strike.ind] <- 1
  df04 <<- rbind(df04,
                 data.frame(session=session, condition=condition,
                            event=1:length(d$peaks),
                            time=sort(c(d$rewarded.time, d$non.rewarded.time)),
                            is.rewarded=is.rewarded,
                            peak=d$peaks))
}

summary.strikes <- function(my.strikes) {
  num.strikes <- length(my.strikes) - 1
  which.strikes <- 1:num.strikes + 1
  time.stamps.rewarded <- c()
  time.stamps.nonrewarded <- c()
  rewarded.strike.ind <- c()
  for (i in which.strikes) {
    if (my.strikes[[i]]$is.rewarded) {
      time.stamps.rewarded <- c(time.stamps.rewarded, my.strikes[[i]]$strike.time)
      rewarded.strike.ind <- c(rewarded.strike.ind, i - 1) # because of header
    } else {
      time.stamps.nonrewarded <- c(time.stamps.nonrewarded,
                                   my.strikes[[i]]$strike.time)
    }
  }
  time.stamps <- sort(c(time.stamps.rewarded, time.stamps.nonrewarded))
  num.rewarded.strikes <- length(time.stamps.rewarded)

  peaks <- c()
  for (i in which.strikes) {
    peaks <- c(peaks, max(my.strikes[[i]]$force))
  }
  
  summary.d <- list(
                 animal.id=my.strikes$header$animal.id,
                 strike.dir=my.strikes$header$strike.dir,
                 rewarded.strike.ind=rewarded.strike.ind,
                 rewarded.time=time.stamps.rewarded,
                 non.rewarded.time=time.stamps.nonrewarded,
                 peaks=peaks)
  return(summary.d)
}

aggregate.action.counts <- function(dirs=N.oerstedii.4.dirs,
                                    animal.id="Y4", species="N.oerstedii",
                                    force.scale=15.03){
  df <- data.frame()
  session <- 1
  for (i in dirs) {
    cat('reading: ', paste(i), '\n')
    strikes <- read.strike(strike.dir=i, animal.id=animal.id,
                           species=species, force.scale=force.scale)
    df <- rbind(df,
                data.frame(session=session,
                           rewarded.time=summary(strikes)$rewarded.time,
                           action.counts=1:length(summary(strikes)$rewarded.time)))
    session <- session + 1
  }
  return(df)
}




## functions for physiological data (text data exported from chart)
## class: phys

read.phys <- function(phys.dir=NA, file.name="strike"){
  wd <- getwd()
  setwd(phys.dir)
  files <- list.files()
  files <- files[order(as.numeric(gsub('.txt', '', gsub(file.name, '', files))))]
  phys <- list()
  for (i in files) {
    one.phys <- read.table(i)
    phys <- c( phys, list(one.phys = list(strike.id=i, phys=one.phys)))
    cat(i, "\n")
  }
  setwd(wd)
  class(phys) <- "phys"
  return(phys)
}

print.phys <- function(phys){
  cat(length(phys), 'physiological data.', '\n')
}

plot.phys <- function(phys, ch=4){
  strike.num <- length(phys)
  ys <- data.frame()
  for (i in 1:strike.num) {
    y <- phys[[i]]$phys[,ch]
    time <- phys[[i]]$phys$V1 - min(phys[[i]]$phys$V1)
    ys <- rbind(ys, data.frame(strike.id=i, time=time, y=y))
  }
  xyplot(y~time | strike.id, ys, strip=FALSE, type="l")
}

strike2df <- function(strike){
  event.n <- length(strike) - 1
  time <- strike[[1]]$time
  phys.df <- data.frame()
  for (i in 2:(event.n+1)) {
    event <- i - 1
    y <- strike[[i]]$voltage
    phys.df <- rbind(phys.df, data.frame(event=event, time=time, y=y, peak=max(y)))
  }
  return(phys.df)
}

align.phys.peak <- function(phys, ch, w=c(-0.003, 0.003), bw=c(-0.003, -0.002)){
  n <- length(phys)
  phys2 <- data.frame()
  for (i in 1:n) {
    y <- phys[[i]]$phys[, ch]
    t <- phys[[i]]$phys[, 1] # time should be channel one
    peak <- max(y)
    t <- t - t[y==peak][1]
    base.line.y <- y[t >= bw[1] & t <= bw[2]]
    base.line.t <- t[t >= bw[1] & t <= bw[2]]
    y <- y - mean(base.line.y)
    base.line.y <- base.line.y - mean(base.line.y)
    peak <- max(y)

    phys2 <- rbind(phys2, data.frame(
                            event = i,
                            time=t[t >= w[1] & t <= w[2]],
                            y=y[t >= w[1] & t <= w[2]],
                            peak=peak))
  }
  p <- xyplot(y~time | event, phys2, type="l", strip=FALSE)
  print(p)
  return(phys2)
}

select.emg.onset <- function(phys){
  strike.num <- length(phys)
  emg.onsets <- data.frame()
  for (i in 1:strike.num) {
    emg <- phys[[i]]$phys$V4
    time <- phys[[i]]$phys$V1
    plot(time, emg, type="l")
    emg.onset <- locator(1)$x
    emg.onsets <- c(emg.onsets, emg.onset)
    time <- time - emg.onset
    browser()
  }
}

## read Spike2 data which is exported wave marks
read.diff.strike <- function(){
## read differentiated strike data into "diff.strikes" in the current workspace
  previous.working.dir <- setwd(tk_choose.dir())
  strike.dir <- getwd()
  cat(list.files(), sep="\n")
  diff.strikes.name <- file.choose()
  cat("Enter animal id and push return key. \n")
  animal.id <- scan( what="" )
  diff.strikes.raw <- readLines( diff.strikes.name )
  wave.initial.lines <- c()
  strike.time <- c()
  sampling.rate <- c()
  strike.code <- c()
  for ( i in 1:( length( diff.strikes.raw ))) {
    if ( grepl( "WAVMK", diff.strikes.raw[i] )){
      wave.info <- unlist( strsplit( diff.strikes.raw[i], split = "\t" ))
      strike.time <- c(strike.time, as.numeric( wave.info[2] ))
      if (length(sampling.rate) == 0){
        sampling.rate <- as.numeric( wave.info[3] )
      }
      strike.code <- c( strike.code, as.numeric( wave.info[4] ))
      wave.initial.lines <- c( wave.initial.lines, i+1 )
    }
  }

  wave.length <- diff( wave.initial.lines )[1] - 3
  diff.strikes <<- data.frame()
  for ( i in 1:length( wave.initial.lines )){
    diff.strikes <<- rbind( diff.strikes,        
                           data.frame(
                             animal.id = animal.id,
                             session = diff.strikes.name,
                             strike.id = i,
                             strike.code = strike.code[i],
                             strike.time = strike.time[i],
                             time = seq(0, by = sampling.rate, length.out = wave.length),
                             wave = as.numeric( diff.strikes.raw[ wave.initial.lines[i]:
                               ( wave.initial.lines[i] + wave.length - 1)] )))
  }
  setwd(previous.working.dir)
  cat("Differentiated strike data were read into 'diff.strikes' in the current work space.\n")
}


#### EMG analysis

align.phys.onset <- function(phys, ch=4, r=10, bl=200, pdf=FALSE) {
  starts <- c()
  rec.starts <- c()
  for (i in 1:length(phys)) {
    bouts <- bout.detect(phys[[i]], ch, r, bl)
    starts <- c(starts, bouts$start[1])
    rec.starts <- c(rec.starts, min(bouts$time))
  }
  rec.start <- max(rec.starts)
  bout.set <- list()
  for (i in 1:length(phys)) {
    phys[[i]]$phys[,1] <- phys[[i]]$phys[,1] - starts[i]
    ## trim phys by rec.start ##
    bouts2 <- bout.detect(phys[[i]], ch, r, bl, event=i, pdf=pdf)
    bout.set <- c(bout.set, list(bouts2))
   }
  return(bout.set)
}

bout.detect <- function(y=phys[[1]], ch=4, r=10, bl=200, pdf=FALSE, event=1){
  y0 <- y$phys[,ch]
  t <- y$phys[,1]
  yl <- length(y0)
  l <- yl %% r
  if (l!=0){ 
    for (i in 1:l) {
      y0 <- y0[-1]
      t <- t[-1]
    }
  }
  yl <- length(y0)
  m <- yl / r
  if (m %% 2 == 1){
    y0 <- y0[-((length(y0):(length(y0) - r + 1)))]
    t <- t[-((length(t):(length(t) - r + 1)))]
    yl <- yl - r
    m <- m - 1
  }
  yy <- rep(0, m)
  tt <- rep(0, m/2)
  for (i in 1:m){
    yy[i] <- y0[i*r]
  }
  for (i in 1:m/2) {
    tt[i] <- t[i*r*2] # tvvar
  }
  z <- tvvar(yy, 2, 6.6e-06, 1.0e-06, FALSE)
  z1 <- ngsmth(z$ts, 2, 2.6e-04,, 2, 1.644934, k=200, plot=FALSE)
  
  ## detect 'on' and 'off' timing of the activity from abrupt positive change
  zt <- z1$trend[,4]
  zd <- diff(zt)
  b <- mean(zd[1:5])
  bsd <- sd(zd[1:5])
  on <- which(zd > b + 100*bsd)
  first <- c()
  middles <- c()
  last <- c()
  cluster.i <- which(diff(on) > 10) # 10 points silence in trend
  if (length(cluster.i) >= 1) {
    first.c <- 1:cluster.i[1]
    first <- which(zd == max(zd[on[first.c]]))
    last.c <- (cluster.i[length(cluster.i)] + 1):length(on)
    last <- which(zd == max(zd[on[last.c]]))
    if (length(cluster.i) >= 2) {
      middles <- c()
      for ( i in 2:length(cluster.i)) {
        if (i >= 2) {
          cluster  <- (cluster.i[i-1] + 1):cluster.i[i]
          middles <- c(middles, which(zd == max(zd[on[cluster]])))
        }
      }
    }
    on2 <- c(first, middles, last)
  } else {
    on2 <- on
  }

  ## 'off'
  off <- which(zd < b - 100*bsd)
  cluster.i <- which(diff(off) > 10) # 10 points silence in trend
  if (length(cluster.i) >= 1) {
    first.c <- 1:cluster.i[1]
    first <- which(zd == min(zd[off[first.c]]))
    last.c <- (cluster.i[length(cluster.i)] + 1):length(off)
    last <- which(zd == min(zd[off[last.c]]))
    if (length(cluster.i) >= 2) {
      middles <- c()
      for ( i in 2:length(cluster.i)) {
        if (i >= 2) {
          cluster  <- (cluster.i[i-1] + 1):cluster.i[i]
          middles <- c(middles, which(zd == min(zd[off[cluster]])))
        }
      }
    }
    off2 <- c(first, middles, last)
  } else {
    off2 <- off
  }
  
  on3 <- on2 * 2 * r + l
  off3 <- off2  * 2 * r + l

  off4 <- c()
  for ( i in 1:length(on3)) {
    if (i == length(on3)) {
      if (length(off3[off3 > on3[i]]) >= 1) {
        off4 <- c(off4, min(off3[off3 > on3[i]]))
      } else {
        on3 <- on3[-length(on3)]
      }
    } else {
      off4 <- c(off4, min(off3[off3 > on3[i] & off3 < on3[i+1]]))
    }
  }
  
  par(mfrow=c(3,1))
  plot(t, y0, type="l", xlab="time (sec)", ylab="EMG (voltage)", main="bout detection")

  for (i in 1:length(on3)) {
    rect(t[on3[i]], min(y0), t[off4[i]], max(y0), col="gray", border="transparent")
    lines(t, y0)
  }
  
  plot(tt, z$ts, type="l", xlab="time (sec)", ylab="logarithm of variance)")
  lines(x=tt, y=zt, col="gray", lwd=2)

  plot(tt[-length(tt)], zd, type="b", xlab="time (sec)", ylab="")
  abline(v=t[on3], col="red")
  abline(v=t[off4], col="blue")

  if (pdf==TRUE){
    dev.copy2pdf(file=paste("bout_", event, ".pdf", sep=""))
  }
  return(list(time=t, y=y0, start=t[on3], termination=t[off4], z=z))
}

my.f <- function(){
  source("strike.probe.r")
}