library("plyr")
library("readr")
library("ggplot2")

setwd("~/Documents/Behavioral Analyses/csv")

master <- read.csv("Vid_metadata.csv")
master_names <- master$File

mussel_names <- vector(mode = "list", length = length(master_names))

#####historical bar graph

historical <- read.csv("freq_barplot.csv")

date <- as.numeric(historical$Date)
number <- as.numeric(historical$N)
lure <- as.factor(historical$Lure)
colors <- c("gray","orange","gray","orange","gray","orange","gray","orange","gray","orange")
fills <-  c("gray","orange","gray","orange","gray","orange","gray","orange","gray","orange")

data <- data.frame(date,lure,number)

ggplot(data, aes(fill=lure, y=number, x=date)) +
		theme_bw() + 
    geom_bar(position="stack", stat="identity", fill = fills) +
    ylab("Number") +
    xlab("Year") +
    geom_abline(intercept=0,slope=0) +
    scale_x_continuous(breaks=100)

chidata <- read.csv("chi.csv",header=TRUE)
chisq.test(chidata$historical,chidata$contemporary)



##initialize calculation loop
setwd("mussels")
for(i in 1:length(master_names)){
		#load sample data
	x <- paste(master_names[i], "data.csv", sep = "_")
	master_names[i] <- as.character(master_names[i])
	individual_mussel <- read.csv(x)
	individual_mussel_dat <- individual_mussel[3:nrow(individual_mussel),]


	#creates logical vectors for event data
	left = individual_mussel_dat[,1]=="l"
	both = individual_mussel_dat[,1]=="b"
	right = individual_mussel_dat[,1]=="r"

	left_num = sum(left)
	right_num = sum(right)
	both_num = sum(both)

	left_prop = left_num/(left_num+both_num)
	print(left_prop)

	right_prop = right_num/(right_num+both_num)
	print(right_prop)

	both_prop = both_num/(left_num+right_num+both_num)
	print(both_prop)

	left_index = left+both
	right_index = right+both
	##both_index = both :)

	left_sync_prop = left_index/both
	right_sync_prop = right_index/both


	#creates empty matrix for L/R data
	l_matrix = matrix(0,nrow(individual_mussel_dat),2)
	r_matrix = matrix(0,nrow(individual_mussel_dat),2)

	#Turns test_dat into a "matrix"
	individual_mussel_dat = as.matrix(individual_mussel_dat)


	#Seperates number of left beats from right beats
	for(j in 1:nrow(individual_mussel_dat)){
		if(left_index[j]==1){
			l_matrix[j,1] = individual_mussel_dat[j,2]
			l_matrix[j,2] = individual_mussel_dat[j,3]

		}

		if(right_index[j]==1){
			r_matrix[j,1] = individual_mussel_dat[j,2]
			r_matrix[j,2] = individual_mussel_dat[j,3]
		}
		
		
	}
	l_matrix[l_matrix==0]  <- NA
	l_matrix = na.omit(l_matrix)


	r_matrix[r_matrix==0] <- NA
	r_matrix = na.omit(r_matrix)


	#calculates speed of undulations	
	left_frame_time = as.numeric(l_matrix[,2])-as.numeric(l_matrix[,1])
	right_frame_time = as.numeric(r_matrix[,2])-as.numeric(r_matrix[,1])
	left_frame_time_fps = left_frame_time/120
	right_frame_time_fps = right_frame_time/120

	#calculates time between undulations
	left_start = as.numeric(l_matrix[,1])
	left_interim_fps = rep(0,length(left_start)-1)
	for(k in 1:(length(left_interim_fps))){
		left_interim_fps[k]=(left_start[k+1]-left_start[k])/120


	}
	
	right_start = as.numeric(r_matrix[,1])
	right_interim_fps = rep(0,length(right_start)-1)
	for(g in 1:(length(right_interim_fps))){
		right_interim_fps[g]=(right_start[g+1]-right_start[g])/120
	}

	

	#reports data
	master[i,2] = mean(left_frame_time_fps)
	master[i,3] = mean(right_frame_time_fps)
	master[i,4] = mean(left_interim_fps)
	master[i,5] = mean(right_interim_fps)
	master[i,6] = sd(left_interim_fps)
	master[i,7] = sd(right_interim_fps)
	master[i,8] = both_prop

}



#############Plot Data

#Plot data

boxplot(master$freq_left~master$Lure, ylab = "Average interval between movements of left mantle flap (s)", xlab = "")
boxplot(master$freq_right~master$Lure, ylab = "Average beat interval of right mantle flap (s)", xlab = "")
boxplot(master$st_dev_left~master$Lure, ylab = "Standard deviation of movement interval of left mantle flap (s)", xlab = "")
boxplot(master$st_dev_right~master$Lure, ylab = "Standard deviation of beat interval of right mantle flap (s)", xlab = "")
boxplot(master$avg_len_left~master$Lure, ylab = "Average duration of left mantle flap movement (s)", xlab = "")
boxplot(master$avg_len_right~master$Lure, ylab = "Average  of right mantle flap (s)", xlab = "")
boxplot(master$prop_sync~master$Lure, ylab = "Proportion of movements that are left-right synchronized", xlab = "")


kruskal.test(master$freq_left~master$Lure)
pairwise.wilcox.test(master$freq_left,master$Lure, p.adjust.method="BH")

kruskal.test(master$freq_right~master$Lure)
pairwise.wilcox.test(master$freq_right,master$Lure, p.adjust.method="BH")

kruskal.test(master$st_dev_left~master$Lure)
pairwise.wilcox.test(master$st_dev_left,master$Lure, p.adjust.method="BH")

kruskal.test(master$st_dev_right~master$Lure)
pairwise.wilcox.test(master$st_dev_right,master$Lure, p.adjust.method="BH")

plot(master$freq_left~master$Temp, xlab = "Temperature (c)", ylab = "Average beat frequency of left mantle flap (s)")
cor.test(master$freq_left,master$Temp, method="spearman")

plot(master$avg_len_left~master$Temp, ylab="Average interval between movements of left mantle flap (s)", xlab="Temperature (c)")

#####historical bar graph

historical <- read.csv(frequency_data.csv)

leech_ex <- read.csv("GH010579_data.csv")
leech_ex <- leech_ex[3:nrow(leech_ex),]

o.side <- as.character(leech_ex[,1])
L.o <- o.side == "l"
R.o <- o.side == "r"
B.o <- o.side == "b"

o.start <- as.numeric(as.character(leech_ex[,2])) - 5000
o.end <- as.numeric(as.character(leech_ex[,3])) - 5000
o.gait <- rep("black", length(o.start))
for (j in 1:length(o.gait)){
	if(o.side[j]=="b"){
		o.gait[j]="red"
	}
}

#Gait diagrams

orange_ex <- read.csv("GH010579_data.csv")
orange_ex <- orange_ex[3:nrow(orange_ex),]

o.side <- as.character(orange_ex[,1])
L.o <- o.side == "l"
R.o <- o.side == "r"
B.o <- o.side == "b"

o.start <- as.numeric(as.character(orange_ex[,2])) - 5000
o.end <- as.numeric(as.character(orange_ex[,3])) - 5000
o.gait <- rep("black", length(o.start))
for (j in 1:length(o.gait)){
	if(o.side[j]=="b"){
		o.gait[j]="red"
	}
}

o.values <- rep(0,length(o.start))

for (x in (1:length(o.start))){
	if(o.side[x] == "l"){
		o.values[x] = 1
	} else{
		if(o.side[x] == "r"){
		o.values[x] = -1
	} else{
		o.values[x] = 0
	}
}
}


darter_ex <- read.csv("GH010077_data.csv")
darter_ex <- darter_ex[3:nrow(darter_ex),]

d.side <- as.character(darter_ex[,1])
L.d <- d.side == "l"
R.d <- d.side == "r"
B.d <- d.side == "b"

d.start <- as.numeric(as.character(darter_ex[,2])) - 5000
d.end <- as.numeric(as.character(darter_ex[,3])) - 5000
d.gait <- rep("black", length(d.start))
for (j in 1:length(d.gait)){
	if(d.side[j]=="b"){
		d.gait[j]="red"
	}
}


d.values <- rep(0,length(d.start))

for (x in (1:length(d.start))){
	if(d.side[x] == "l"){
		d.values[x] = 1
	} else{
		if(d.side[x] == "r"){
		d.values[x] = -1
	} else{
		d.values[x] = 0
	}
}
}


#Cardium

cardium_ex <- read.csv("GH010060_data.csv")
cardium_ex <- cardium_ex[3:nrow(cardium_ex),]

c.side <- as.character(cardium_ex[,1])
L.c <- c.side == "l"
R.c <- c.side == "r"
B.c <- c.side == "b"

c.start <- as.numeric(as.character(cardium_ex[,2])) - 5000
c.end <- as.numeric(as.character(cardium_ex[,3])) - 5000
c.gait <- rep("black", length(c.start))
for (j in 1:length(c.gait)){
	if(c.side[j]=="b"){
		c.gait[j]="red"
	}
}


c.values <- rep(0,length(c.start))

for (x in (1:length(c.start))){
	if(c.side[x] == "l"){
		c.values[x] = 1
	} else{
		if(c.side[x] == "r"){
		c.values[x] = -1
	} else{
		c.values[x] = 0
	}
}
}


 

#plotting gait analysis

par(mfrow=c(3,1),mar=c(2, 25, 4, 2), xpd=NA)

plot(y=d.values, x=d.start,xlim=c(min(d.start),max(10000)),pch=2,cex=.001,ylab="",xlab="",axes=F, ylim=c(-5,5))

axis(side=1, at=c(0, 2000, 4000, 6000,8000,10000))
for (j in 1:length(d.start)){
  segments(d.start[j],d.values[j],d.end[j],d.values[j], lend=1, lwd=20,col=d.gait[j])
}
abline(h=0, lty=2)

plot(y=o.values, x=o.start,xlim=c(min(o.start),max(10000)),pch=2,cex=.001,ylab="",xlab="",axes=F, ylim=c(-5,5))

axis(side=1, at=c(0, 2000, 4000, 6000,8000,10000))
for (j in 1:length(o.start)){
  segments(o.start[j],o.values[j],o.end[j],o.values[j],lend=1, lwd=20,col=o.gait[j])
}
abline(h=0, lty=2)

plot(y=c.values, x=c.start,xlim=c(min(c.start),max(10000)),pch=2,cex=.001,ylab="",xlab="",axes=F, ylim=c(-5,5))

axis(side=1, at=c(0, 2000, 4000, 6000,8000,10000))
for (j in 1:length(c.start)){
  segments(c.start[j],c.values[j],c.end[j],c.values[j],lend=1, lwd=20,col=c.gait[j])
}
abline(h=0, lty=2)

#creates function for combining columns of different lengths, substituting NA
cbindPad <- function(...){
args <- list(...)
n <- sapply(args,nrow)
mx <- max(n)
pad <- function(x, mx){
    if (nrow(x) < mx){
        nms <- colnames(x)
        padTemp <- matrix(NA, mx - nrow(x), ncol(x))
        colnames(padTemp) <- nms
        if (ncol(x)==0) {
          return(padTemp)
        } else {
        return(rbind(x,padTemp))
          }
    }
    else{
        return(x)
    }
}
rs <- lapply(args,pad,mx)
return(do.call(cbind,rs))
}

df <- c("ID_vector","Phenotype","Site","Species","Left_move_duration", "Left_rest","Right_move_duration","Right_rest")
boot_df <- c("ID_vector","Phenotype","Site","Species","Left_move_duration", "Left_rest","Right_move_duration","Right_rest","Left_interval","Right_interval")

#for loop to create the main dataframe
for(i in 1:length(names)){
		#load sample data
	x <- paste(names[i], "data.csv", sep = "_")
	names[i] <- as.character(names[i])
	individual_mussel <- read.csv(x)
	individual_mussel_dat <- individual_mussel[3:nrow(individual_mussel),]

	ID <- as.character(names[i])
	morph <- as.character(individual_mussel[1,1])
	site <- as.character(individual_mussel[1,2])
	species <- if (morph=="cardium"){
		species = "cardium" } else {
			species = "fasciola"
		}
	

	m.side <- as.character(individual_mussel_dat[,1])

	L.m <- m.side == "l"
	B.m <- m.side == "b"
	R.m <- m.side == "r"

	Left_index <- L.m + B.m
	Right_index <- R.m + B.m

	l_matrix = matrix(0,nrow(individual_mussel_dat),2)
	r_matrix = matrix(0,nrow(individual_mussel_dat),2)

	individual_mussel_dat <- as.matrix(individual_mussel_dat)

	for(j in 1:nrow(individual_mussel_dat)){
		if(Left_index[j]==1){
			l_matrix[j,1] = individual_mussel_dat[j,2]
			l_matrix[j,2] = individual_mussel_dat[j,3]
		}

		if (Right_index[j]==1){
			r_matrix[j,1] = individual_mussel_dat[j,2]
			r_matrix[j,2] = individual_mussel_dat[j,3]
		}
	}

	l_matrix[l_matrix==0] <- NA
	l_matrix = na.omit(l_matrix)

	r_matrix[r_matrix==0] <- NA
	r_matrix = na.omit(r_matrix)

	left_start = as.numeric(l_matrix[,1])
	left_end <- as.numeric(l_matrix[,2])

	left_move_duration = rep(0,length(left_start))
	for(k in 1:(length(left_move_duration))){
		left_move_duration[k]=(left_end[k]-left_start[k])/120

	}
	
	left_rest = rep(0,length(left_start)-1)
	for(k in 1:(length(left_rest))){
		left_rest[k]=(left_start[k+1]-left_end[k])/120
	}

	right_start = as.numeric(r_matrix[,1])
	right_end <- as.numeric(r_matrix[,2])

	right_move_duration= rep(0,length(right_start))
	for(g in 1:(length(right_move_duration))){
		right_move_duration[g]=(right_end[g]-right_start[g])/120
	}

	right_rest = rep(0,length(right_start)-1)
	for(k in 1:(length(right_rest))){
		right_rest[k]=(right_start[k+1]-right_end[k])/120
	}

	left_interval <- rep(0,length(left_start)-1)
	for(y in 1:(length(left_interval))){
		left_interval[y]=(left_start[y+1]-left_start[y])/120
	}

	right_interval <- rep(0,length(right_start)-1)
	for(z in 1:(length(right_interval))){
		right_interval[z]=(right_start[z+1]-right_start[z])/120
	}

	boot_L <- as.numeric(sample(left_move_duration,1000,replace=TRUE))
	boot_Lrest <- as.numeric(sample(left_rest,1000,replace=TRUE))
	boot_R <- as.numeric(sample(right_move_duration,1000,replace=TRUE))
	boot_Rrest <- as.numeric(sample(left_rest,1000,replace=TRUE))
	boot_Linterval <- as.numeric(sample(left_interval,1000,replace=TRUE))
	boot_Rinterval <- as.numeric(sample(right_interval,1000,replace=TRUE))

	if(length(left_move_duration) > length(right_move_duration)){
		ID_vector <- rep(ID,length(left_move_duration))
		Phenotype <- rep(morph,length(left_move_duration))
		Site <- rep(site,length(left_move_duration))
		Species <- rep(species,length(left_move_duration))
		} else {
			ID_vector <- rep(ID,length(right_move_duration))
			Phenotype <- rep(morph,length(right_move_duration))
			Site <- rep(site,length(right_move_duration))
			Species <- rep(species,length(right_move_duration))
		}

	boot_ID <- rep(ID,1000)
	boot_Phenotype <- rep(morph,1000)
	boot_Site <- rep(site,1000)
	boot_Species <- rep(species,1000)


	df <- rbind(df,cbindPad(as.data.frame(ID_vector), as.data.frame(Phenotype), as.data.frame(Site), as.data.frame(Species),as.data.frame(left_move_duration), as.data.frame(left_rest),as.data.frame(right_move_duration),as.data.frame(right_rest)))
	boot_df <-rbind(boot_df,cbind(as.data.frame(boot_ID), as.data.frame(boot_Phenotype), as.data.frame(boot_Site),as.data.frame(boot_Species) ,as.data.frame(boot_L),as.data.frame(boot_Lrest),as.data.frame(boot_R),as.data.frame(boot_Rrest),as.data.frame(boot_Linterval),as.data.frame(boot_Rinterval)))
}

df <- df[-1,]
boot_df <- boot_df[-1,]

Lmove <- as.numeric(df$left_move_duration)
Lrest <- as.numeric(df$left_rest)
Rmove <- as.numeric(df$right_move_duration)
Rrest <- as.numeric(df$right_rest)

type <- df$Phenotype
id <- df$ID_vector

lmm_left_move <- lmer(Lmove~type + (1 | id), data = df)
lmm_left_rest <- lmer(Lrest~type + (1 | id), data = df)
lmm_right_move <- lmer(Rmove~type + (1 | id), data = df)
lmm_right_rest <- lmer(Rrest~type + (1 | id), data = df)



boot_Lmove <- as.numeric(boot_df$boot_L)
boot_Lrest <- as.numeric(boot_df$boot_Lrest)
boot_Rmove <- as.numeric(boot_df$boot_R)
boot_Rrest <- as.numeric(boot_df$boot_Rrest)
boot_Linterval <- as.numeric(boot_df$boot_Linterval)
boot_Rinterval <- as.numeric(boot_df$boot_Rinterval)

boot_type <- as.factor(boot_df$boot_Phenotype)
boot_id <- as.factor(boot_df$boot_ID)
boot_site <- as.factor(boot_df$boot_Site)

lmm_left_move_boot <- lmer(boot_Lmove~boot_type + (1 | boot_id))
lmm_left_rest_boot <- lmer(boot_Lrest~boot_type + (1 | boot_id))
lmm_right_move_boot <- lmer(boot_Rmove~boot_type + (1 | boot_id))
lmm_right_rest_boot <- lmer(boot_Rrest~boot_type + (1 | boot_id))
lmm_left_interval_boot <- lmer(boot_Linterval~boot_type + (1 | boot_id))
lmm_right_interval_boot <- lmer(boot_Rinterval~boot_type + (1 | boot_id))

#p values for models

coefs.left_rest_boot <- data.frame(coef(summary(lmm_left_rest_boot)))

#histographs for all individuals for each lure type

darter <- boot_df[boot_df$boot_Phenotype=="darter",]
hist(as.numeric(darter$boot_Lrest), breaks=100)

orange <- boot_df[boot_df$boot_Phenotype=="orange",]
hist(as.numeric(orange$boot_Lrest), breaks=100)

cardium <- boot_df[boot_df$boot_Phenotype=="cardium",]
hist(as.numeric(cardium$boot_Lrest), breaks=100)