## Density Sep search analysis
## N Mizumoto
## 10/24/2019 initiate
## 10/27/2019 major update
## 12/13/2019 5FPS major revision
## 02/25/2020 Format for publication
## 05/28/2020 Revision start

## Environment
{
  RPROJ <- list(PROJHOME = normalizePath(getwd()))
  attach(RPROJ)
  rm(RPROJ)
}

## setting
{ # skip to main
  options(warn=0)
  
  today <- Sys.Date()
  
  ## packages
  library(ggplot2)
  library(lme4)
  library(car)
  library(multcomp)
  library(survival)
  library(multcomp)
  library("survminer")
  
  ## housekeeping
  rm(list = ls())
  
  ## fuction
  se  <-  function(x){
    y  <-  x[!is.na(x)]  #  remove  the  missing  values
    sqrt(var(as.vector(y))/length(y))
  }
  
  angle_cal <- function(X, Y, Length){
    Ax <- (X[3:Length-1] - X[3:Length-2])
    Bx <- (X[3:Length] - X[3:Length-1])
    Ay <- (Y[3:Length-1] - Y[3:Length-2])
    By <- (Y[3:Length] - Y[3:Length-1])
    hugo <- (Ax * By - Ay * Bx + 0.000001)/abs(Ax * By - Ay * Bx + 0.000001)
    cos <- round((Ax * Bx + Ay * By) / ((Ax^2 + Ay^2)*(Bx^2 + By^2))^0.5,14)
    return(acos(cos)*hugo)
  }
  
  tandem.smoothing <- function(vec, xvec, min.sec){
    if(sum(vec)>0){
      timing <- which(vec)[c(T, diff(which(vec))>1)]
      end <- which(vec)[c(diff(which(vec))>1,T)]
      for(fi in 1:length(timing)){
        if(length( vec[timing[fi]:end[fi]]) < min.sec ){
          vec[timing[fi]:end[fi]] <- F
        }
      }
    }
    return(vec)
  }
  
  tandem.smoothing2 <- function(vec, xvec, min.sec){
    if(sum(vec)>0){
      timing <- which(vec)[c(T, diff(which(vec))>1)]
      end <- which(vec)[c(diff(which(vec))>1,T)]
      for(fi in 1:length(timing)){
        if(any(xvec[timing[fi]:end[fi]])){
          vec[timing[fi]:end[fi]] <- T
        } else if(length( vec[timing[fi]:end[fi]]) < min.sec ){
          vec[timing[fi]:end[fi]] <- F
        }
      }
    }
    return(vec)
  }
  
  logistic <- function(x){ 1/(1+exp(-x)) }
  
  data_summrize <- function(d, variable, cname){
    l <- length(cname)
    if(l>=2){
      df <- data.frame(as.vector(tapply(d[,cname][,1], d[,cname], unique)))
      for(i in 2:l){
        df <- data.frame(df,
                         as.vector(tapply(d[,cname][,i], d[,cname], unique))
        )
      }
    } else {
      df <- data.frame(as.vector(tapply(d[,cname], d[,cname], unique)))
    }
    df <- data.frame(df,
                     as.vector(tapply(d[,variable], d[,cname], mean)),
                     as.vector(tapply(d[,variable], d[,cname], sd)),
                     as.vector(tapply(d[,variable], d[,cname], se)))
    colnames(df) <- c(cname, paste0(variable, ".", c("mean", "sd", "se")))
    return(df)
  }
  
  ## data
  TandemRdata <- file.path ( PROJHOME , "data", "tandem", "rdata") # location of 5FPS data (formatted in .Rdata file)
  files = list.files (TandemRdata)
}

## parameters
{
  options(warn=2)

  analFPS = 5
  analDuration = 10 #minutes
  inter_thresh = 9.75 # threshold for tandem interactions 7.5*1.3
  noninter_thresh = inter_thresh*2 # threshold for no tandem interactions
  minsec <- 5 # minimum time for tandem runs
  minsepsec <- 2 # minimum time for separation event
  thresh_dis = 30 # minumum moved distance for tandem

  preobs <- 5
  postobs <- 10
  extraobs <- 120*analFPS
  
  good.4.col <- c("#ff2800","#fab01b","#66ccff","#0041ff") #c("#ff2800","#35a16b","#66ccff","#0041ff")
  good.4.col2 <- c("#cf1816ff","#f9b811ff","#66ccff","#0041ff")
  plot.traj = F
  bg.set <- c("white","white","white","white")
  pch.set <- c(16,21,24,17)
}

## Main calcuration for separation events
{
  res = res2 = res.ind = res.tan.time = res.sep.time = res.fem.dir = res.ma.dir <- NULL
  event.count = event.count2 = 0
  for(v in 1:length(files)){
    name <- substr(files[v],  1, 5)
    treat <- substr(files[v],  1, 2)
    
    print(paste(v, name))
    
    load(file.path(TandemRdata, files[v]))
    d <- d.video
    d <- d[d[,1,1]<=60*analDuration,,]
    
    dims = dim(d)
    frames <- dims[1]
    num_ind <- dims[3]
    
    ####### 1. identifying tandem #######
    dist.mat <- array(0, c(dims[1], num_ind, num_ind))
    for(i in 1:num_ind){
      dist.mat[,,i] <- sqrt((d[,2,1:num_ind] - d[,2,i])^2 + (d[,3,1:num_ind] - d[,3,i])^2)
    }
    
    # 1-1. two individuals need to be closer
    # interaction: < 10 mm
    interaction.mat = dist.mat < inter_thresh & dist.mat > 0
    interaction.mat.both = dist.mat < inter_thresh & dist.mat > 0
    for(i in 1:num_ind){
      if(i <= num_ind/2){
        interaction.mat[,1:(num_ind/2),i] <- F
      } else {
        interaction.mat[,(num_ind/2+1):num_ind,i] <- F
      }
    }
    interaction.each = apply(interaction.mat, c(1,3), any)
    interaction.each.both = apply(interaction.mat.both, c(1,3), any)
    interaction.num = apply(interaction.mat, c(1,3), sum)
    
    
    interaction.partner=interaction.mat*1
    for(i in 1:num_ind){
      interaction.partner[,i,] = interaction.mat[,i,] * i
    }
    interaction.partner.mat = interaction.partner
    interaction.partner <- apply(interaction.partner, c(1,3), max)
    
    tandem = interaction.each
    
    # 1-2-1. not non interaction  & less than minsepsec seconds = tandem
    noninteraction.mat = dist.mat >= noninter_thresh & dist.mat != 0
    for(i in 1:num_ind){
      if(i <= num_ind/2){
        noninteraction.mat[,1:(num_ind/2),i] <- T
      } else {
        noninteraction.mat[,(num_ind/2+1):num_ind,i] <- T
      }
    }
    noninteraction.each = apply(noninteraction.mat, c(1,3), all)
    for(i in 1:num_ind){
      if(sum(!tandem[,i])>0){
        tandem[,i] = !tandem.smoothing2(!tandem[,i], noninteraction.each[,i], minsepsec*analFPS)
      }
    }
    
    ## 1-2-2. smoothing (tandem needs to last 5 seconds)
    for(i in 1:num_ind){
      if(sum(tandem[,i])>0){
        tandem[,i] = tandem.smoothing(tandem[,i], noninteraction.each[,i], minsec*analFPS)
      }
    }
    
    ## 1-2-3. tandem: both need to move more than thresh_dis (30 mm) during interaction
    moved.dis.each <- sqrt(diff(d[,2,1:num_ind])^2+diff(d[,3,1:num_ind])^2)
    moved.dis.each2 <- rbind(0, moved.dis.each)
    for(i in 1:num_ind){
      if(sum(tandem[,i])>0){
        tan.timing <- which(tandem[,i])[c(T, diff(which(tandem[,i]))>1)]
        tan.end <- which(tandem[,i])[c(diff(which(tandem[,i]))>1,T)]
        for(j in 1:length(tan.timing)){
          if(sum(moved.dis.each2[tan.timing[j]:tan.end[j],i])< thresh_dis){
            tandem[tan.timing[j]:tan.end[j],i] <- F
          }
          if(sum(diag(moved.dis.each2[tan.timing[j]:tan.end[j],interaction.partner[tan.timing[j]:tan.end[j],i]]))< thresh_dis){
            tandem[tan.timing[j]:tan.end[j],i] <- F
          }
        }
      }
    }
    
    # fix inconsistenct caused by above smoothing process
    for(i in 1:num_ind){
      for(j in 2:length(tandem[,1])){
        if(tandem[j,i] & interaction.partner[j,i]==0){
          interaction.partner[j,i] = interaction.partner[j-1,i]
        }
      }
    }
    
    #########################################
    
    ####### 2. get scheme #########
    # grab data for 1 minutes after tandem interrupt
    scheme = duration = tandemduration <- tandem
    scheme[interaction.each] <- "i" # interaction
    scheme[!interaction.each] <- "r" # random search (others)
    scheme[tandem] <- "t" # tandem run
    
    for(i in 1:num_ind){
      if(sum(tandem[,i])>0){
        
        tan.end <- which(tandem[,i])[c(diff(which(tandem[,i]))>1,T)]
        if(tan.end[1]==frames){ next; }
        sep.begin <- tan.end[tan.end<frames]+1
        
        ## sep until next interaction with the others for more than 1sec
        if(sum(interaction.each.both[,i])>0){
          interaction.each.both[,i] = tandem.smoothing(interaction.each.both[,i], noninteraction.each[,i], 1*analFPS)
        }
        tan.timing <- which(interaction.each.both[,i])[c(T, diff(which(interaction.each.both[,i]))>1)]
        if(length(tan.timing)==1){ sep.end <- NULL} else {
          sep.end <- tan.timing[2:(length(tan.timing))]-1}
        if(length(sep.begin) - length(sep.end) == 1){ sep.end <- c(sep.end,frames)}
        for(j in 1:length(sep.begin)){
          if(sep.begin[j]+1*analFPS <= frames){
            if(sum(!interaction.each.both[sep.begin[j]:(sep.begin[j]+1*analFPS),i]) == 0){
              next;
            }
          }
          for(k in 1:length(sep.end)){
            if(sep.end[k]>=sep.begin[j]){
              scheme[sep.begin[j]:sep.end[k],i] <- "s" # sepration search
              break;
            }
          }
        }
      }
    }
    ##############################
    
    ##### angle #####
    d2 <- d
    angle.mat <-  moved.dis.each2
    for(i in 1:num_ind){
      angle.mat[,i] = c(angle_cal(c(d[1,2,i],d2[,2,i]), c(d[1,3,i],d2[,3,i]), dim(d2)[1]+1), NA)
    }
    #angle.mat = abs(angle.mat)
    
    ###### Plot trajectories? ######
    if(plot.traj){
      png ( file.path ( PROJHOME , "Plots" , "Trajectories", paste0(name, ".png" )), 
            pointsize=8, width=9600, height=4800, res=1200)
      par(mfrow=c(1,2), pin=c(3,3))
      plot( d[,2,1], d[,3,1], type="l", col=2, xlim=c(0,140), ylim=c(0,140), main=name,
            xlab="x(mm)", ylab="y(mm)")
      for(i in 2:num_ind){
        if(i < num_ind/2+1){lines( d[,2,i], d[,3,i], col=2)} else {lines( d[,2,i], d[,3,i], col=1)}
      }
      #pmoved.dis.each2 = moved.dis.each2
      #pmoved.dis.each2[scheme!="t"] <- NA
      
      sep.begin <- NULL
      for(i in 1:num_ind){
        if(sum(tandem[,i])>0){
          tan.end <- which(tandem[,i])[c(diff(which(tandem[,i]))>1,T)]
          if(tan.end[1]==frames){ next; }
          sep.begin <- c(sep.begin, tan.end[tan.end<frames]+1)
        }
      }
      sep.begin <- unique(sep.begin)
      matplot(moved.dis.each2, type="o", pch=1, ylim=c(0,200), col=rep(c(2,5),each=num_ind/2), 
              main=paste(v,name, "tandem:", sum(scheme=="t"), "sep:", sum(scheme=="s")), ylab="Speed (mm/sec)", xlab="Time (sec)", las=1)
      if(sum(scheme=="s")){
        arrows(sep.begin, rep(70, length(sep.begin)), sep.begin, rep(40,length(sep.begin)), length=0.2)
      }
      dev.off()
    }
    ##############################
    
    ###### Output data ######
    sex <- rep(c("f","m"), each=num_ind/2)
    ## For separation events
    {
      Video = Ind = Pair = Sex = Event = Scheme = Time = RealTime = Speed = Angle = Suround = X = Y <- NULL
      for(i in 1:num_ind){
        if(sum(scheme[,i]=="s")>0){
          sep.timing <- which(scheme[,i]=="s")[c(T, diff(which(scheme[,i]=="s"))>1)]
          sep.end <- which(scheme[,i]=="s")[c(diff(which(scheme[,i]=="s"))>1,T)]
          tan.timing <- which(tandem[,i])[c(T, diff(which(tandem[,i]))>1)]
          
          for(j in 1:length(sep.timing)){
            event.count <- event.count + 1
            Video = c(Video, rep(name, sep.end[j]-tan.timing[j]+1) )
            Ind = c(Ind, rep(paste(name,i,sep="-"), sep.end[j]-tan.timing[j]+1) )
            Pair = c(Pair, rep(num_ind/2, sep.end[j]-tan.timing[j]+1) )
            Sex = c(Sex, rep(sex[i], sep.end[j]-tan.timing[j]+1) )
            Event = c(Event, rep(event.count, sep.end[j]-tan.timing[j]+1) )
            Scheme = c(Scheme, scheme[tan.timing[j]:sep.end[j],i])
            Time = c(Time, ((tan.timing[j]-sep.timing[j]):(sep.end[j]-sep.timing[j]))/analFPS)
            RealTime = c(RealTime, (tan.timing[j]:sep.end[j])/analFPS)
            Speed = c(Speed, moved.dis.each2[tan.timing[j]:sep.end[j],i])
            Angle = c(Angle, angle.mat[tan.timing[j]:sep.end[j],i])
            Suround = c(Suround, interaction.num[tan.timing[j]:sep.end[j],i])
            X = c(X, d[tan.timing[j]:sep.end[j],2,i])
            Y = c(Y, d[tan.timing[j]:sep.end[j],3,i])
          }
        }
      }
      df <- data.frame(Video, Ind, Pair, Sex, Event, Scheme, Time, RealTime, Speed, Angle, Suround, X, Y)
      res <- rbind(res, df)
    }
    
    ## tandem survival time
    {
      Video = Ind = Pair = Sex = Tan.time = Cens <- NULL
      if(sum(scheme=="t")>0){
        for(i in 1:num_ind){
          if(sum(scheme[,i]=="t")>0){
            tan.timing <- which(tandem[,i])[c(T, diff(which(tandem[,i]))>1)]
            tan.end <- which(tandem[,i])[c(diff(which(tandem[,i]))>1,T)]
            for(j in 1:length(tan.timing)){
              Tan.time = c(Tan.time, tan.end[j] - tan.timing[j] + 1)
              if(tan.timing[j]==1 || tan.end[j]==frames){
                Cens = c(Cens, 0)
              } else {
                Cens = c(Cens, 1)
              }
              Video = c(Video, name)
              Ind = c(Ind, paste(name,i,sep="-"))
              Pair = c(Pair, num_ind/2)
              Sex = c(Sex, sex[i])
            }
          }
        }
        df <- data.frame(Video, Ind, Pair, Sex, Tan.time=Tan.time/analFPS, Cens)
        res.tan.time <- rbind(res.tan.time, df)
      }
    }
    
    ## sep search survival time
    {
      if(sum(scheme=="s")>0){
        Video = Ind = Pair = Sex = Sep.time = Cens = Change <- NULL
        for(i in 1:num_ind){
          if(sum(scheme[,i]=="s")>0){
            sep.timing <- which(scheme[,i]=="s")[c(T, diff(which(scheme[,i]=="s"))>1)]
            sep.end <- which(scheme[,i]=="s")[c(diff(which(scheme[,i]=="s"))>1,T)]
            for(j in 1:length(sep.timing)){
              Sep.time = c(Sep.time, sep.end[j] - sep.timing[j] + 1)
              if(sep.timing[j]==1 || sep.end[j]==frames){
                Cens = c(Cens, 0)
              } else {
                Cens = c(Cens, 1)
              }
              if(interaction.partner[sep.end[j]+1,i]>0){
                tan.bef <- interaction.partner[(sep.timing[j]-10):(sep.timing[j]-1),i]
                tan.bef.partner <- as.numeric(names(rev(sort(table(tan.bef))))[1])
                if(interaction.partner[sep.end[j]+1,i] == tan.bef.partner){
                  Change = c(Change, 0)
                } else {
                   Change = c(Change, 1)
                }
              } else {
                Change = c(Change, NA)
              }
              Video = c(Video, name)
              Ind = c(Ind, paste(name,i,sep="-"))
              Pair = c(Pair, num_ind/2)
              Sex = c(Sex, sex[i])
            }
          }
        }
        df <- data.frame(Video, Ind, Pair, Sex, Sep.time=Sep.time/analFPS, Cens, Change)
        res.sep.time <- rbind(res.sep.time, df)
      }
    }
    
    ## Female direction from male upon separation
    if(sum(scheme=="s")>0){
      Video = Ind = Pair = Sex = Sep.time = FemDir <- NULL
      for(i in (num_ind/2+1):num_ind){
        if(sum(scheme[,i]=="s")>0){
          sep.timing <- which(scheme[,i]=="s")[c(T, diff(which(scheme[,i]=="s"))>1)]
          mx = d[sep.timing,2,i]
          my = d[sep.timing,3,i]
          mpx = d[sep.timing-1,2,i]
          mpy = d[sep.timing-1,3,i]
          theta2 <- atan2(my-mpy, mx-mpx)
          for(j in 1:length(mx)){
            fx = d[sep.timing,2,interaction.partner[sep.timing-1,i][j]][j]
            fy = d[sep.timing,3,interaction.partner[sep.timing-1,i][j]][j]
            theta1 <- atan2( fy-my[j], fx-mx[j])
            FemDir <- c(FemDir, theta2[j]-theta1)
            Video = c(Video, name)
            Ind = c(Ind, paste(name,i,sep="-"))
            Pair = c(Pair, num_ind/2)
            Sex = c(Sex, sex[i])
          }
        }
      }
      df <- data.frame(Video, Ind, Pair, Sex, FemDir)
      res.fem.dir <- rbind(res.fem.dir, df)
    }
    
    ## Male direction from female upon separation
    if(sum(scheme=="s")>0){
      Video = Ind = Pair = Sex = Sep.time = MaDir <- NULL
      for(i in 1:(num_ind/2)){
        if(sum(scheme[,i]=="s")>0){
          sep.timing <- which(scheme[,i]=="s")[c(T, diff(which(scheme[,i]=="s"))>1)]
          fx = d[sep.timing,2,i]
          fy = d[sep.timing,3,i]
          fpx = d[sep.timing-1,2,i]
          fpy = d[sep.timing-1,3,i]
          theta2 <- atan2(fy-fpy, fx-fpx)
          for(j in 1:length(fx)){
            mx = d[sep.timing,2,interaction.partner[sep.timing-1,i][j]][j]
            my = d[sep.timing,3,interaction.partner[sep.timing-1,i][j]][j]
            theta1 <- atan2( my-fy[j], mx-fx[j])
            MaDir <- c(MaDir, theta2[j]-theta1)
            Video = c(Video, name)
            Ind = c(Ind, paste(name,i,sep="-"))
            Pair = c(Pair, num_ind/2)
            Sex = c(Sex, sex[i])
          }
        }
      }
      df <- data.frame(Video, Ind, Pair, Sex, MaDir)
      res.ma.dir <- rbind(res.ma.dir, df)
    }
  }
}

## Analysis
## Speed time development
{
  ## 1 pair vs multiple pairs
  if(T){
    yrange = c(0,7)
    options(warn=0)
    res.anal <- res
    res.anal$Pair <- (res.anal$Pair>1) * 1
    
    ## data
    rmeanf = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], mean)[,,"f"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rmeanm = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], mean)[,,"m"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rsef = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], se)[,,"f"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rsem = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], se)[,,"m"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rsef = rsef * 2
    rsem = rsem * 2
    
    par(mfrow=c(2,1), pin=c(4.854,3))
    for(i in 1:2){
      plot(0,type="n",ylim=yrange, xlim=c(-preobs,postobs),
           xlab="Time (sec)", ylab="Dispacement (mm/0.2sec)", las=1, main="")
      abline(v=-0)
      par(new=T)
      matplot(round(seq(-preobs,postobs,1/analFPS),1), rmeanf[,i], type="o",ylim=yrange, xlim=c(-preobs,postobs),
              col=good.4.col[1], axes=F, xlab="", ylab="",
              pch=pch.set, bg=bg.set, las=1, cex=0.5)
      polygon(c(round(seq(-preobs,postobs,1/analFPS),1),round(seq(postobs,-preobs,-1/analFPS),1)),
              c((rmeanf+rsef)[,i],rev((rmeanf-rsef)[,i])),
                col=adjustcolor(good.4.col[1], 0.5), border=NA)
      par(new=T)
      matplot(round(seq(-preobs,postobs,1/analFPS),1), rmeanm[,i], type="o",ylim=yrange, xlim=c(-preobs,postobs),
              col=good.4.col[4], axes=F, xlab="", ylab="",
              pch=pch.set, bg=bg.set, las=1, cex=0.5)
      polygon(c(round(seq(-preobs,postobs,1/analFPS),1),round(seq(postobs,-preobs,-1/analFPS),1)), 
              c((rmeanm+rsem)[,i],rev((rmeanm-rsem)[,i])),
                col=adjustcolor(good.4.col[4], 0.5), border=NA)
    }
  }
  
  ## each pair
  if(T){
    yrange = c(0,10)
    options(warn=0)
    res.anal <- res
    
    ## data
    rmeanf = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], mean)[,,"f"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rmeanm = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], mean)[,,"m"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rsef = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], se)[,,"f"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rsem = na.omit(tapply(res.anal$Speed, res.anal[,c("Time","Pair","Sex")], se)[,,"m"])[paste(round(seq(-preobs,postobs,1/analFPS),1)),]
    rsef = rsef*2
    rsem = rsem*2
    
    {
      par(mfrow=c(2,2), pin=c(4.854,3)*0.75)
      for(i in 1:4){
        plot(0,type="n",ylim=yrange, xlim=c(-preobs,postobs),
             xlab="Time (sec)", ylab="Speed (mm/sec)", las=1, main="")
        abline(v=-0.1)
        par(new=T)
        matplot(round(seq(-preobs,postobs,1/analFPS),1), rmeanf[,i], type="o",ylim=yrange, xlim=c(-preobs,postobs),
                col=good.4.col[1], axes=F, xlab="", ylab="",
                pch=pch.set, bg=bg.set, las=1, cex=0.5)
        polygon(c(round(seq(-preobs,postobs,1/analFPS),1),round(seq(postobs,-preobs,-1/analFPS),1)),
                c((rmeanf+rsef)[,i],rev((rmeanf-rsef)[,i])),
                col=adjustcolor(good.4.col[1], 0.5), border=NA)
        par(new=T)
        matplot(round(seq(-preobs,postobs,1/analFPS),1), rmeanm[,i], type="o",ylim=yrange, xlim=c(-preobs,postobs),
                col=good.4.col[4], axes=F, xlab="", ylab="",
                pch=pch.set, bg=bg.set, las=1, cex=0.5)
        polygon(c(round(seq(-preobs,postobs,1/analFPS),1),round(seq(postobs,-preobs,-1/analFPS),1)), 
                c((rmeanm+rsem)[,i],rev((rmeanm-rsem)[,i])),
                col=adjustcolor(good.4.col[4], 0.5), border=NA)
      }
    }
  }
  
  ## Speed comparison just after separation
  iF(T){
    res.just.after = res[res$Time>0 & res$Time<= 2,]
    res.just.after.event = data.frame(res.just.after[res.just.after$Time==0.2,1:5], 
          Speed = na.omit(as.vector(tapply(res.just.after$Speed, res.just.after[,c("Sex","Pair","Event")], mean))))
    
    res.just.after.event$Group <- res.just.after.event$Pair>1
    
    # male
    df <- res.just.after.event[res.just.after.event$Sex=="m",]
    df$Pair <- as.factor(df$Pair)
    r <- lmer(Speed~Pair+(1|Ind), data=df)
    Anova(r)
    #Analysis of Deviance Table (Type II Wald chisquare tests)
    #Response: Speed
    #Chisq Df Pr(>Chisq)   
    #Pair 15.571  3   0.001388 **
    multicomparison<-glht(r,linfct=mcp(Pair="Tukey"))
    summary(multicomparison)
    #Simultaneous Tests for General Linear Hypotheses
    #Multiple Comparisons of Means: Tukey Contrasts
    #Fit: lmer(formula = Speed ~ Pair + (1 | Ind), data = df)
    #Linear Hypotheses:
    #  Estimate Std. Error z value Pr(>|z|)    
    #2 - 1 == 0  1.905419   0.706429   2.697   0.0332 *  
    #4 - 1 == 0  1.685694   0.578550   2.914   0.0179 *  
    #8 - 1 == 0  1.900430   0.487561   3.898   <0.001 ***
    #4 - 2 == 0 -0.219726   0.664400  -0.331   0.9869    
    #8 - 2 == 0 -0.004989   0.586876  -0.009   1.0000    
    #8 - 4 == 0  0.214737   0.424384   0.506   0.9557  
    
    # female
    df <- res.just.after.event[res.just.after.event$Sex=="f",]
    df$Pair <- as.factor(df$Pair)
    r <- lmer(Speed~Pair+(1|Ind), data=df)
    Anova(r)
    #Analysis of Deviance Table (Type II Wald chisquare tests)
    #Response: Speed
    #Chisq Df Pr(>Chisq)
    #Pair 0.0614  3      0.996
    
    errbar_lims <- data_summrize(res.just.after.event, "Speed", colnames(res.just.after.event)[3:4])
    ggplot() +
      geom_point(data=errbar_lims, aes(x=as.factor(Pair), y=Speed.mean), size=3) +
      geom_errorbar(aes(x=as.factor(errbar_lims$Pair), 
                        ymax=errbar_lims$Speed.mean+errbar_lims$Speed.se*2, 
                        ymin=errbar_lims$Speed.mean-errbar_lims$Speed.se*2), stat='identity',
                    width=.1)+
      theme_bw() + theme(aspect.ratio = 1.8) +
      ylim(0,7) + ylab("Displacements (mm)") + xlab("Pair")
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-SpeedJustAfter.pdf")), width=3)
  }
  
  ## Speed comparison later after separation
  iF(T){
    res.just.after = res[res$Time>2 & res$Time<= 10,]
    res.just.after.event = data.frame(res.just.after[res.just.after$Time==2.2,1:5], 
                                      Speed = na.omit(as.vector(tapply(res.just.after$Speed, res.just.after[,c("Sex","Pair","Event")], mean))))
    
    res.just.after.event$Group <- res.just.after.event$Pair>1
    
    # male
    df <- res.just.after.event[res.just.after.event$Sex=="m",]
    df$Pair <- as.factor(df$Pair)
    r <- lmer(Speed~Pair+(1|Ind), data=df)
    Anova(r)
    #Analysis of Deviance Table (Type II Wald chisquare tests)
    #Response: Speed
    #Chisq Df Pr(>Chisq)   
    #Pair 15.571  3   0.001388 **
    
    # female
    df <- res.just.after.event[res.just.after.event$Sex=="f",]
    df$Pair <- as.factor(df$Pair)
    r <- lmer(Speed~Pair+(1|Ind), data=df)
    Anova(r)
    #Analysis of Deviance Table (Type II Wald chisquare tests)
    #Response: Speed
    #Chisq Df Pr(>Chisq)
    #Pair 0.0614  3      0.996
    
    errbar_lims <- data_summrize(res.just.after.event, "Speed", colnames(res.just.after.event)[3:4])
    ggplot() +
      geom_point(data=errbar_lims, aes(x=as.factor(Pair), y=Speed.mean), size=3) +
      geom_errorbar(aes(x=as.factor(errbar_lims$Pair), 
                        ymax=errbar_lims$Speed.mean+errbar_lims$Speed.se*2, 
                        ymin=errbar_lims$Speed.mean-errbar_lims$Speed.se*2), stat='identity',
                    width=.1)+
      theme_bw() + theme(aspect.ratio = 1.8) +
      ylim(0,8) + ylab("Displacements (mm)") + xlab("Pair")
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-SpeedJustAfter.pdf")), width=3)
  }
}

## Survival analysis
{
    ## Tandem
    table(res.tan.time[res.tan.time$Cens==1,4:3])
    #Pair
    #Sex   1   2   4   8
    #f  96  68 130 621
    #m  96  70 146 714
    
    df<-survfit(Surv(Tan.time,Cens)~Pair, type = "kaplan-meier", data=res.tan.time[res.tan.time$Cens==1,])
    r <- coxph(Surv(Tan.time,Cens)~Pair*Sex, data=res.tan.time[res.tan.time$Cens==1,])
    Anova(r)
    #Analysis of Deviance Table (Type II tests)
    #LR Chisq Df Pr(>Chisq)    
    #Pair       84.261  1     <2e-16 ***
    #Sex         1.100  1     0.2943    
    #Pair:Sex    0.091  1     0.7624
    summary(r)
    
    tapply(res.tan.time[res.tan.time$Cens==1,"Tan.time"], res.tan.time[res.tan.time$Cens==1,"Pair"], median)
    #1    2    4    8 
    #63.4 48.8 30.3 22.4
    
    # ggplot
    df<-survfit(Surv(Tan.time,Cens)~Pair, type = "kaplan-meier", data=res.tan.time[res.tan.time$Cens==1,])
    g <- ggsurvplot(fit = df, data = res.tan.time,
               pval = F, pval.method = TRUE,
               risk.table = F, conf.int = T,
               ncensor.plot = F, size = 0.2, 
               xlab="Time (sec)", palette = good.4.col)
    g$plot + theme_classic() + theme(panel.border = element_rect(colour = "black", fill=NA), aspect.ratio = 3/4, legend.position = 'none')
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-TandemSurv.pdf")), width=3)
    
    
    ## SepSearch
    table(res.sep.time[res.sep.time$Cens==1,4:3])
    #Pair
    #Sex   1   2   4   8
    #f 112  72 142 595
    #m 112  68 146 668
    
    df<-survfit(Surv(Sep.time,Cens)~Pair, type = "kaplan-meier", data=res.sep.time[res.sep.time$Cens==1,])
    r <- coxph(Surv(Sep.time,Cens)~Pair*Sex, data=res.sep.time[res.sep.time$Cens==1,])
    Anova(r)
    #Analysis of Deviance Table (Type II tests)
    #LR Chisq Df Pr(>Chisq)    
    #Pair      116.316  1  < 2.2e-16 ***
    #Sex        37.757  1  8.013e-10 ***
    #Pair:Sex    1.141  1     0.2854 
    
    errbar_lims <- data_summrize(res.sep.time, "Sep.time", colnames(res.sep.time)[4:3])
    
    par(mfrow=c(1,1))
    XonFig <- barplot(errbar_lims$Sep.time.mean, beside=T, ylim=c(0,25),
                      col=adjustcolor(good.4.col[rep(c(1,4),4)], 0.5),
                      xlab="Pair", ylab="Time (sec)", las=1)
    arrows(XonFig , errbar_lims$Sep.time.mean, XonFig , errbar_lims$Sep.time.mean + errbar_lims$Sep.time.se*2, angle=90, length=0.1)
    box()
    
    for(i in c(1,2,4,8)){
    print(survdiff(Surv(Sep.time,Cens)~Sex, rho=1, data=res.sep.time[res.sep.time$Pair==i,]))
    }
    
    #Call:1
    #Chisq= 0  on 1 degrees of freedom, p= 1 
    #Call:2
    #Chisq= 0  on 1 degrees of freedom, p= 0.9
    #Call:4
    #Chisq= 11.7  on 1 degrees of freedom, p= 6e-04 
    #Call:8
    #Chisq= 19.3  on 1 degrees of freedom, p= 1e-05 
}

## Number of changing partners
{
  par(mfrow=c(1,2),pin=c(45,40)/15)
  switch.table <- table(res.sep.time[,c("Change","Pair","Sex")])
  print(switch.table)
  
  switch.prob <- (switch.table[2,,] / apply(switch.table,c(2,3),sum))
  par(pin=c(45,40)/15)
  matplot(cbind(c(1,2,4,8)-0.1,c(1,2,4,8)+0.1), switch.prob, ylim=c(0,0.8), pch=c(21,16), col=good.4.col[c(1,4)],
          xlab="Number of pairs", ylab="Probability to switch partner", las=1)
  legend(1,0.8,c("Female","Male"),col=good.4.col[c(1,4)],pch=c(21,16))
}

## Angle
{
  library(CircStats)
  rho = rho.low = rho.upp = pair = Sex = time.angle <- NULL
  for(p in c(1,2,4,8)){
    for(sex in c("f","m")){
      for(i in seq(-preobs,postobs,0.2)){
        d.angle = na.omit(res[res$Time==round(i,1) & res$Sex == sex & res$Pair == p, "Angle"])
        r <- wrpcauchy.ml(d.angle, 0, 0, acc=1e-015)
        
        r.per <- rep(0,1000)
        for(j in 1:1000){
          d.angle2 <- sample(d.angle, length(d.angle), replace=T)
          r.res <- wrpcauchy.ml(d.angle2, 0, 0, acc=1e-015)
          r.per[j] <- r.res$rho
        }
        
        
        time.angle <- c(time.angle, i) 
        rho <- c(rho, r$rho)
        rho.low <- c(rho.low, quantile(r.per,0.025))
        rho.upp <- c(rho.upp, quantile(r.per,0.975))
        pair <- c(pair, p)
        Sex <- c(Sex, sex)
      }
    }
  }
  df <- data.frame(pair, Sex, time.angle, rho, rho.low, rho.upp)
  res.angle <- df
  #save(res.angle, file="200530_res.angle.rda")
  
  load("200530_res.angle.rda")
  par(pin=c(5,4)/2, mfrow=c(2,2))
  yrange <- c(0.3,1)
  for(i in c(1,2,4,8)){
    df <- res.angle[res.angle$pair==i,]
    df.m <- df[df$Sex=="m",]
    df.f <- df[df$Sex=="f",]
    
    plot(0,type="n",ylim=yrange, xlim=c(-preobs,postobs),
         xlab="Time (sec)", ylab="Dispacement (mm/0.2sec)", las=1, main="")
    abline(v=-0)
    par(new=T)
    matplot(round(seq(-preobs,postobs,1/analFPS),1), df.f$rho, type="o",ylim=yrange, xlim=c(-preobs,postobs),
            col=good.4.col[1], axes=F, xlab="", ylab="",
            pch=pch.set, bg=bg.set, las=1, cex=0.5)
    polygon(c(round(seq(-preobs,postobs,1/analFPS),1),round(seq(postobs,-preobs,-1/analFPS),1)),
            c((df.f$rho.upp),rev((df.f$rho.low))),
            col=adjustcolor(good.4.col[1], 0.5), border=NA)
    par(new=T)
    matplot(round(seq(-preobs,postobs,1/analFPS),1), df.m$rho, type="o",ylim=yrange, xlim=c(-preobs,postobs),
            col=good.4.col[4], axes=F, xlab="", ylab="",
            pch=pch.set, bg=bg.set, las=1, cex=0.5)
    polygon(c(round(seq(-preobs,postobs,1/analFPS),1),round(seq(postobs,-preobs,-1/analFPS),1)), 
            c(df.m$rho.upp,rev(df.m$rho.low)),
            col=adjustcolor(good.4.col[4], 0.5), border=NA)
  }
  
  res.angle[res.angle$time.angle>=0,]
  df <- res.angle[res.angle$time.angle>0 & res.angle$time.angle<=2,]
  tapply(df$rho, df[,c("pair","Sex")], mean)
  df <- res.angle[res.angle$time.angle>2 & res.angle$time.angle<=10,]
  tapply(df$rho, df[,c("pair","Sex")], mean)
  df <- res.angle[res.angle$time.angle<0,]
  tapply(df$rho, df[,c("pair","Sex")], mean)
  
}

rho = rho.low = rho.upp = pair = Sex = time.angle <- NULL
for(p in c(1,2,4,8)){
  for(sex in c("f","m")){
    #d.angle = na.omit(res[res$Time>-5 & res$Time<=-2 & res$Sex == sex & res$Pair == p, "Angle"])
    d.angle = na.omit(res[res$Time>2 & res$Time<=10 & res$Sex == sex & res$Pair == p, "Angle"])
    r <- wrpcauchy.ml(d.angle, 0, 0, acc=1e-015)
    r.per <- rep(0,1000)
    for(j in 1:1000){
      d.angle2 <- sample(d.angle, length(d.angle), replace=T)
      r.res <- wrpcauchy.ml(d.angle2, 0, 0, acc=1e-015)
      r.per[j] <- r.res$rho
    }
    #time.angle <- c(time.angle, "tandem") 
    time.angle <- c(time.angle, "after1") 
    rho <- c(rho, r$rho)
    rho.low <- c(rho.low, quantile(r.per,0.025))
    rho.upp <- c(rho.upp, quantile(r.per,0.975))
    pair <- c(pair, p)
    Sex <- c(Sex, sex)
  }
}
df <- data.frame(pair, Sex, time.angle, rho, rho.low, rho.upp)
#df.tandem <- df
#df.after1 <- df     
#df.after2 <- df
df.rho <- rbind(df.tandem, df.after1, df.after2)
#save(df.rho, file="200602_res.rho.rda")
df.rho.male <- df.rho[df.rho$Sex=="m",]
posn.d <- position_dodge(width = 0.5) # ずらすのに必要
ggplot(df.rho.male, aes(x=time.angle, group=pair, y=rho)) + 
  geom_point( position = posn.d, cex=2) + 
  geom_errorbar(aes(ymin=rho.low, ymax=rho.upp), width=0.1, position = posn.d) +
  ylim(c(0.5,0.85)) + 
  theme_bw() + theme(aspect.ratio = 1.5)
ggsave(file.path(PROJHOME,paste0("Plots/",today,"-rhoMale.pdf")), width=3, height=4)

## Heading direction
{
  ## Male heading direction to females
  {
    res.fem.dir$FemDir[res.fem.dir$FemDir<0] <- res.fem.dir$FemDir[res.fem.dir$FemDir<0] + 2*pi
    
    ggplot(res.fem.dir[res.fem.dir$Pair==1,], aes(x=FemDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-5,30))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-FemaleDir-1p.pdf")), width=3, height=4)
    
    ggplot(res.fem.dir[res.fem.dir$Pair==2,], aes(x=FemDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-5,30))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-FemaleDir-2p.pdf")), width=3, height=4)
    
    ggplot(res.fem.dir[res.fem.dir$Pair==4,], aes(x=FemDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-6,40))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-FemaleDir-4p.pdf")), width=3, height=4)
    
    ggplot(res.fem.dir[res.fem.dir$Pair==8,], aes(x=FemDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-20,200))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-FemaleDir-8p.pdf")), width=3, height=4)
    
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==1,"FemDir"]))
    # Test Statistic:  0.191 
    # P-value:  0.0168 
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==2,"FemDir"]))
    # Test Statistic:  0.5619 
    # P-value:  0 
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==4,"FemDir"]))
    # Test Statistic:  0.4675 
    # P-value:  0 
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==8,"FemDir"]))
    # Test Statistic:  0.48 
    # P-value:  0 
  }
  ## Female heading direction to males
  {
    res.ma.dir$MaDir[res.ma.dir$MaDir<0] <- res.ma.dir$MaDir[res.ma.dir$MaDir<0] + 2*pi
    
    ggplot(res.ma.dir[res.ma.dir$Pair==1,], aes(x=MaDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-5,40))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-MaleDir-1p.pdf")), width=3, height=4)
    
    ggplot(res.ma.dir[res.ma.dir$Pair==2,], aes(x=MaDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-5,30))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-MaleDir-2p.pdf")), width=3, height=4)
    
    ggplot(res.ma.dir[res.ma.dir$Pair==4,], aes(x=MaDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-6,40))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-MaleDir-4p.pdf")), width=3, height=4)
    
    ggplot(res.ma.dir[res.ma.dir$Pair==8,], aes(x=MaDir))+
      geom_histogram(fill='blue', color='black', breaks=seq(0,pi*2,length=11)) +
      coord_polar(start=pi/2*3, direction=-1) +
      scale_x_continuous("", limits = c(0, pi*2), 
                         breaks = seq(0,pi*7/4,pi/4),
                         labels = c(0, "pi/4", "pi/2", "3pi/4", "pi", "-3pi/4", "-pi/2", "-pi/4")) +
      theme_bw() + theme(aspect.ratio = 1) + scale_y_continuous("",limits=c(-20,200))
    ggsave(file.path(PROJHOME,paste0("Plots/",today,"-MaleDir-8p.pdf")), width=3, height=4)
    
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==1,"FemDir"]))
    # Test Statistic:  0.191 
    # P-value:  0.0168 
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==2,"FemDir"]))
    # Test Statistic:  0.5619 
    # P-value:  0 
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==4,"FemDir"]))
    # Test Statistic:  0.4675 
    # P-value:  0 
    rayleigh.test(circular(res.fem.dir[res.fem.dir$Pair==8,"FemDir"]))
    # Test Statistic:  0.48 
    # P-value:  0 
    
    rayleigh.test(circular(res.ma.dir[res.ma.dir$Pair==1,"MaDir"]))
    # Test Statistic:  0.5855 
    # P-value:  0 
    rayleigh.test(circular(res.ma.dir[res.ma.dir$Pair==2,"MaDir"]))
    # Test Statistic:  0.5186 
    # P-value:  0 
    rayleigh.test(circular(res.ma.dir[res.ma.dir$Pair==4,"MaDir"]))
    # Test Statistic:  0.4349 
    # P-value:  0 
    rayleigh.test(circular(res.ma.dir[res.ma.dir$Pair==8,"MaDir"]))
    # Test Statistic:  0.4396 
    # P-value:  0 
  }
}


## Parameters for simulations
iF(T){
  ## Speed
  res.just.after = res[res$Time>0 & res$Time<= 2,]
  res.just.after.event = data.frame(res.just.after[res.just.after$Time==0.2,1:5], 
                                    Speed = na.omit(as.vector(tapply(res.just.after$Speed, res.just.after[,c("Sex","Pair","Event")], mean))))
  res.just.after.event$Group <- res.just.after.event$Pair > 1 
  tapply(res.just.after.event$Speed, res.just.after.event[,c("Group","Sex")], mean)*5
  # Speed of males in a higher density after separation: 27.02146
  # Speed of males in a lower density after separation: 17.04721
  
  res.later.after = res[res$Time>2 & res$Time<= 10,]
  res.later.after.event = data.frame(res.later.after[res.later.after$Time==2.2,1:5], 
                                    Speed = na.omit(as.vector(tapply(res.later.after$Speed, res.later.after[,c("Sex","Pair","Event")], mean))))
  tapply(res.later.after.event$Speed, res.later.after.event[,c("Sex")], mean)*5
  # Speed of males later after separation: 29.375606
  
  res.after = res[res$Time>0 & res$Time<= 10,]
  res.after.event = data.frame(res.after[res.after$Time==0.2,1:5], 
                                    Speed = na.omit(as.vector(tapply(res.after$Speed, res.after[,c("Sex","Pair","Event")], mean))))
  tapply(res.after.event$Speed, res.after.event[,c("Sex")], mean)*5
  # Speed of females after separation: 6.772733
  
  res.just.tandem = res[res$Time>-preobs & res$Time<= -2,]
  res.just.tandem.event = data.frame(res.just.tandem[res.just.tandem$Time==-2.2,1:5], 
                                     Speed = na.omit(as.vector(tapply(res.just.tandem$Speed, res.just.tandem[,c("Sex","Pair","Event")], mean))))
  res.just.tandem.event$Group <- res.just.tandem.event$Pair > 1 
  tapply(res.just.tandem.event$Speed, res.just.tandem.event[,c("Group")], mean)*5
  # Other individuals in a higher density: 12.8750
  
  ## Sinuousity
  res.angle.after.separation <- res.angle[res.angle$time.angle>0,]
  res.angle.just.after.separation <- res.angle.after.separation[res.angle.after.separation$time.angle<=2,]
  tapply(res.angle.just.after.separation$rho, res.angle.just.after.separation[,c("Sex")], mean)
  # Sinuousity (rho) when time <= 2
  # male: 0.6412996
  # female: 0.7179393
  
  res.angle.later.after.separation <- res.angle.after.separation[res.angle.after.separation$time.angle>2,]
  tapply(res.angle.later.after.separation$rho, res.angle.later.after.separation[,c("Sex")], mean)
  # Sinuousity (rho) when time > 2
  # male: 0.7291529
  # female: 0.6867064
  
  res.tandem.sinuousity <- res[res$Scheme=="t",]
  r <- wrpcauchy.ml(na.omit(res.tandem.sinuousity$Angle), 0, 0, acc=1e-015)
  r$rho
  # 0.7799844
}
