Commit 780f3b52 authored by Björn Brembs's avatar Björn Brembs
Browse files

Adjusted for new arena data

Can now read analog motor data as well as light-guides position data
Position data rescaled to +/-1800
Bugfixes and cosmetic improvements.
parent 7bea25a1
No preview for this file type
......@@ -12,9 +12,383 @@ install.packages("data.table")
library("XML")
flyData <- xmlParse(file.choose())
project.file <- file.choose()
devtools::install_github("ACCLAB/dabestr")
install.packages(devtools)
install.packages("devtools")
devtools::install_github("ACCLAB/dabestr")
devtools::install_github("ACCLAB/dabestr")
install.packages("dabestr")
install.packages(c("assertthat", "backports", "BayesFactor", "BH", "bit", "cli", "cluster", "coda", "colorspace", "cowplot", "DBI", "digest", "dplyr", "dygraphs", "effsize", "evaluate", "ggplot2", "ggpubr", "ggsignif", "glue", "gtable", "gtools", "highr", "htmlwidgets", "jsonlite", "knitr", "lazyeval", "lubridate", "markdown", "MASS", "Matrix", "mgcv", "mime", "munsell", "mvtnorm", "nlme", "pbapply", "pillar", "pkgconfig", "purrr", "R6", "Rcpp", "RcppEigen", "RCurl", "rlang", "rmarkdown", "rpart", "RSQLite", "scales", "seewave", "stringi", "stringr", "survival", "tibble", "tidyr", "tidyselect", "tuneR", "utf8", "XML", "xtable", "xts", "yaml", "zoo"))
install.packages(c("nlme", "rpart"))
install.packages(c("nlme", "rpart"))
devtools::install_github("ACCLAB/dabestr")
install.packages("processx")
devtools::install_github("ACCLAB/dabestr")
install.packages("ps")
devtools::install_github("ACCLAB/dabestr")
install.packages("pkgload")
install.packages("dabestr")
devtools::install_github("ACCLAB/dabestr")
install.packages("fs")
devtools::install_github("ACCLAB/dabestr")
install.packages("plyr")
devtools::install_github("ACCLAB/dabestr")
install.packages("datasets")
install.packages("datasets")
install.packages("datasets")
install.packages("datasets")
install.packages("datasets")
install.packages("datasets")
devtools::install_github("ACCLAB/dabestr")
load("D:/data/course_tlearn/evaluations/.RData")
View(rawdata)
View(rawdata)
View(rawdata)
rawdata$a_pos = rawdata$a_pos-1800
View(rawdata)
View(rawdata)
View(rawdata)
View(experiment)
binsize = as.integer(experiment$sample_rate/20)
real_sample_rate = nrow(rawdata)/(rawdata$time[nrow(rawdata)]/1000)
experiment$sample_rate = real_sample_rate/binsize
View(experiment)
seq(1,length(rawdata$time),binsize)
abs(max(rawdata$period[1002:1005]) - min(rawdata$period[1002:1005]))
abs(max(rawdata$period[1000:2000]) - min(rawdata$period[1000:2000]))
abs(max(rawdata$period[5000:7000]) - min(rawdata$period[5000:7000]))
dumydata = rawdata$period[5000:7000]
table(dummydata)
table(dumydata)
dumydata$1
view(dumydata)
View(dumydata)
ntestnumber=table(dumydata)$1
testnumber=table(dumydata)[1,1]
testnumber=table(dumydata)[1,1]
testtable<-table(dumydata)
testtable
testtable$1
testtable[1.1]
testtable[1.1]>testtable[1.2]
testtable[1.2]
testtable[1.1]
testtable[2.1]
testtable[2]
testtable[1]
testtable[1]>testtable[2]
testtable[2]>testtable[1]
table(dumydata)[1]>table(dumydata)[2]
table(dumydata)[2]>table(dumydata)[1]
View(experiment)
# create the vectors in which to save the downsampled data
timeDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$time)/binsize))
a_posDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$a_pos)/binsize))
torqueDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$torque)/binsize))
periodDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$period)/binsize))
# downsampling time
for (index in seq(1,length(rawdata$time),binsize)) {
if(index < (length(rawdata$time)-binsize)) { # check whether we reached the end of the data; if not:
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:length(rawdata$time)])/length(rawdata$time[index:length(rawdata$time)])) # average over the remaining values and save the result
}
}
# downsampling position
for (index in seq(1,length(rawdata$a_pos),binsize)) {
if(index < (length(rawdata$a_pos)-binsize)) { # check whether we reached the end of the data; if not:
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:length(rawdata$a_pos)])/length(rawdata$a_pos[index:length(rawdata$a_pos)])) # average over the remaining values and save the result
}
}
# downsampling torque
for (index in seq(1,length(rawdata$torque),binsize)) {
if(index < (length(rawdata$torque)-binsize)) { # check whether we reached the end of the data; if not:
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:length(rawdata$torque)])/length(rawdata$torque[index:length(rawdata$torque)])) # average over the remaining values and save the result
}
}
# downsampling period
for (index in seq(1,length(rawdata$period),binsize)) {
if (abs(max(rawdata$period[index:(index+binsize-1)]) - min(rawdata$period[index:(index+binsize-1)])) = 0){ #check if there is a period switch in the bin, if there is:
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index]
}
else #if there are two different period values in the bin
{
if (table(rawdata$period[index:(index+binsize-1)])[1]>table(rawdata$period[index:(index+binsize-1)])[2]) #if the majority of values indicates old period
{
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index] #make the period vale that of the old period
}
else {periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index+binsize-1] #if the majority of values indicates new period, set the value to the new period
}
}
}
# bind the downsampled vectors into one dataframe
rawdataDown <- data.frame("time" = timeDownsampled, "a_pos" = a_posDownsampled, "torque" = torqueDownsampled, "period" = periodDownsampled)
# create the vectors in which to save the downsampled data
timeDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$time)/binsize))
a_posDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$a_pos)/binsize))
torqueDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$torque)/binsize))
periodDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$period)/binsize))
# downsampling time
for (index in seq(1,length(rawdata$time),binsize)) {
if(index < (length(rawdata$time)-binsize)) { # check whether we reached the end of the data; if not:
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:length(rawdata$time)])/length(rawdata$time[index:length(rawdata$time)])) # average over the remaining values and save the result
}
}
# downsampling position
for (index in seq(1,length(rawdata$a_pos),binsize)) {
if(index < (length(rawdata$a_pos)-binsize)) { # check whether we reached the end of the data; if not:
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:length(rawdata$a_pos)])/length(rawdata$a_pos[index:length(rawdata$a_pos)])) # average over the remaining values and save the result
}
}
# downsampling torque
for (index in seq(1,length(rawdata$torque),binsize)) {
if(index < (length(rawdata$torque)-binsize)) { # check whether we reached the end of the data; if not:
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:length(rawdata$torque)])/length(rawdata$torque[index:length(rawdata$torque)])) # average over the remaining values and save the result
}
}
# downsampling period
for (index in seq(1,length(rawdata$period),binsize)) {
if (abs(max(rawdata$period[index:(index+binsize-1)]) - min(rawdata$period[index:(index+binsize-1)])) == 0){ #check if there is a period switch in the bin, if there is:
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index]
} else { #if, instead, there are two different period values in the bin...
if (table(rawdata$period[index:(index+binsize-1)])[1]>table(rawdata$period[index:(index+binsize-1)])[2]) #if the majority of values indicates old period
{
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index] #make the period vale that of the old period
} else {
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index+binsize-1] #if the majority of values indicates new period, set the value to the new period
}
}
}
# bind the downsampled vectors into one dataframe
rawdataDown <- data.frame("time" = timeDownsampled, "a_pos" = a_posDownsampled, "torque" = torqueDownsampled, "period" = periodDownsampled)
View(rawdataDown)
View(rawdataDown)
# create the vectors in which to save the downsampled data
timeDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$time)/binsize))
a_posDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$a_pos)/binsize))
torqueDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$torque)/binsize))
periodDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$period)/binsize))
# downsampling time
for (index in seq(1,length(rawdata$time),binsize)) {
if(index < (length(rawdata$time)-binsize)) { # check whether we reached the end of the data; if not:
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:length(rawdata$time)])/length(rawdata$time[index:length(rawdata$time)])) # average over the remaining values and save the result
}
}
# downsampling position
for (index in seq(1,length(rawdata$a_pos),binsize)) {
if(index < (length(rawdata$a_pos)-binsize)) { # check whether we reached the end of the data; if not:
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:length(rawdata$a_pos)])/length(rawdata$a_pos[index:length(rawdata$a_pos)])) # average over the remaining values and save the result
}
}
# downsampling torque
for (index in seq(1,length(rawdata$torque),binsize)) {
if(index < (length(rawdata$torque)-binsize)) { # check whether we reached the end of the data; if not:
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:length(rawdata$torque)])/length(rawdata$torque[index:length(rawdata$torque)])) # average over the remaining values and save the result
}
}
# downsampling period
for (index in seq(1,length(rawdata$period),binsize)) {
if (abs(max(rawdata$period[index:(index+binsize-1)]) - min(rawdata$period[index:(index+binsize-1)])) == 0){ #check if there is a period switch in the bin, if there is:
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index]
} else if (table(rawdata$period[index:(index+binsize-1)])[1]>table(rawdata$period[index:(index+binsize-1)])[2]) #if the majority of values indicates old period
{
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index] #make the period vale that of the old period
} else {
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index+binsize-1] #if the majority of values indicates new period, set the value to the new period
}
}
# bind the downsampled vectors into one dataframe
rawdataDown <- data.frame("time" = timeDownsampled, "a_pos" = a_posDownsampled, "torque" = torqueDownsampled, "period" = periodDownsampled)
length(rawdata$torque)
rawdata$period[length(rawdata$period))
rawdata$period[length(rawdata$period)]
# create the vectors in which to save the downsampled data
timeDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$time)/binsize))
a_posDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$a_pos)/binsize))
torqueDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$torque)/binsize))
periodDownsampled <- vector(mode = "numeric", length = ceiling(length(rawdata$period)/binsize))
# downsampling time
for (index in seq(1,length(rawdata$time),binsize)) {
if(index < (length(rawdata$time)-binsize)) { # check whether we reached the end of the data; if not:
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
timeDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$time[index:length(rawdata$time)])/length(rawdata$time[index:length(rawdata$time)])) # average over the remaining values and save the result
}
}
# downsampling position
for (index in seq(1,length(rawdata$a_pos),binsize)) {
if(index < (length(rawdata$a_pos)-binsize)) { # check whether we reached the end of the data; if not:
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
a_posDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$a_pos[index:length(rawdata$a_pos)])/length(rawdata$a_pos[index:length(rawdata$a_pos)])) # average over the remaining values and save the result
}
}
# downsampling torque
for (index in seq(1,length(rawdata$torque),binsize)) {
if(index < (length(rawdata$torque)-binsize)) { # check whether we reached the end of the data; if not:
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:(index+binsize-1)])/binsize) # average all data in the bin and save it in the right slot of the downsampled vector
} else { # in case we reached the end
torqueDownsampled[((index-1)/binsize)+1] <- round(sum(rawdata$torque[index:length(rawdata$torque)])/length(rawdata$torque[index:length(rawdata$torque)])) # average over the remaining values and save the result
}
}
# downsampling period
for (index in seq(1,length(rawdata$period),binsize)) {
if(index < (length(rawdata$period)-binsize)) { # check whether we reached the end of the data; if not:
if (abs(max(rawdata$period[index:(index+binsize-1)]) - min(rawdata$period[index:(index+binsize-1)])) == 0){ #check if there is a period switch in the bin, if there is:
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index]
} else if (table(rawdata$period[index:(index+binsize-1)])[1]>table(rawdata$period[index:(index+binsize-1)])[2]) #if the majority of values indicates old period
{
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index] #make the period vale that of the old period
} else {
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[index+binsize-1] #if the majority of values indicates new period, set the value to the new period
}
} else { # in case we reached the end
periodDownsampled[((index-1)/binsize)+1] <- rawdata$period[length(rawdata$period)]
}
}
# bind the downsampled vectors into one dataframe
rawdataDown <- data.frame("time" = timeDownsampled, "a_pos" = a_posDownsampled, "torque" = torqueDownsampled, "period" = periodDownsampled)
View(rawdataDown)
table(rawdata$period)
table(rawdataDown$period)
timeownsampled <- timeDownsampled-timeDownsampled[1]
timeDownsampled <- timeDownsampled-timeDownsampled[1]
rawdataDown <- data.frame("time" = timeDownsampled, "a_pos" = a_posDownsampled, "torque" = torqueDownsampled, "period" = periodDownsampled)
View(rawdataDown)
rawdata$diff <- diff(rawdata$time, 1, 1)
difference <- diff(rawdata$time, 1, 1)
view(difference)
View(difference)
plot(difference)
ggplot(difference)
library(ggplot2)
ggplot(difference)
ggplot(as.data.frame(difference))
lines(difference)
lines(difference)#
lines(difference)#
lines(difference)
plot.new#
plot.new
lines(difference)
plot(difference)
View(rawdata)
real_sample_rate
range=vector(min=sapply(rawdata$torque,min),max=sapply(rawdata$torque,max))
range=data.frame(min=sapply(rawdata$torque,min),max=sapply(rawdata$torque,max))
View(range)
rawdata$torque,min
setwd("B:/GitHub/DTSevaluations")
source('B:/GitHub/DTSevaluations/HTML_DTS_project.R', echo=TRUE)
setwd("B:/GitHub/DTSevaluations")
source('B:/GitHub/DTSevaluations/HTML_DTS_project.R', echo=TRUE)
View(pooled.data)
View(pooled.data)
knitr::opts_chunk$set(echo = TRUE)
exp_groups <- as.data.frame(t(plyr::ldply(exp_groups, rbind)))
View(exp_groups)
exp_groups <- list()
exp_groups[[x]] <- c(grp_title, grp_description, xml_list)
View(exp_groups)
plot(exp_groups)
print(exp_groups)
exp_groups <- as.data.frame(t(plyr::ldply(exp_groups, rbind)))
View(exp_groups)
unlist(exp_groups[1,])
colnames(exp_groups) = exp_groups[1, ]
exp_groups = exp_groups[-1,]
colnames(exp_groups) = exp_groups[1, ]
colnames(exp_groups) = exp_groups[1,]
exp_groups <- list()
exp_groups[[x]] <- c(grp_title, grp_description, xml_list)
View(exp_groups)
exp_groups <- as.data.frame(t(plyr::ldply(exp_groups, rbind)))
View(exp_groups)
colnames(exp_groups) = exp_groups[1, ]
View(exp_groups)
backup<-exp_groups
View(backup)
exp_groups = exp_groups[-1, ]
exp_groups<-backup
View(exp_groups)
View(exp_groups)
exp_groups = exp_groups[-2, ]
exp_groups<-backup
exp_groups = exp_groups[-7, ]
exp_groups<-backup
print(exp_groups)
exp_groups = exp_groups[-1, ]
print(exp_groups)
exp_groups<-backup
typeof(exp_groups)
exp_groups <- list()
exp_groups[[x]] <- c(grp_title, grp_description, xml_list)
typeof(exp_groups)
exp_groups <- as.data.frame(t(plyr::ldply(exp_groups, rbind)))
typeof(exp_groups)
exp_groups <- as.data.frame(t(plyr::ldply(exp_groups, rbind))) #make dataframe from list
typeof(exp_groups)
exp_groups <- list()
exp_groups[[x]] <- c(grp_title, grp_description, xml_list)
do.call(rbind.data.frame, exp_groups)
exp_groups<-do.call(rbind.data.frame, exp_groups)
typeof(exp_groups)
View(exp_groups)
typeof(backup)
exp_groups<-backup
as.data.frame(exp_groups)
exp_groups<-as.data.frame(exp_groups)
typeof(backup)
typeof(exp_groups)
View(exp_groups)
class(exp_groups)
View(exp_groups)
rownames(exp_groups) <- NULL
View(exp_groups)
kable(exp_groups)
class(exp_groups)
exp_groups <- exp_groups[-1, ]
class(exp_groups)
print(exp_groups)
kable(exp_groups)
as.data.frame(exp_groups)
exp_groups<-backup
kable(exp_groups)
exp_groups$test = c(1,2,3,4,5,6)
exp_groups$test = c(7,6,5,4,3,2,12)
kable(exp_groups)
exp_groups <- exp_groups[-1, ]
kable(exp_groups)
exp_groups = exp_groups[-1, drop=F]
kable(exp_groups)
exp_groups<-backup
exp_groups = exp_groups[-1, ,drop=F]
kable(exp_groups)
rownames(exp_groups) <- NULL
kable(exp_groups)
View(PIprofile)
View(sequence)
setwd("B:/GitHub/DTSevaluations")
source('B:/GitHub/DTSevaluations/HTML_DTS_project.R', echo=TRUE)
View(exp_groups)
View(PIprofile)
View(grouped.PIprofiles)
grouped.PIprofiles[[x]] = PIprofile
setwd("B:/GitHub/DTSevaluations")
source('B:/GitHub/DTSevaluations/HTML_DTS_project.R', echo=TRUE)
meanspec(rawdata$torque, f = samplerate, wl = 512)
......@@ -8,17 +8,17 @@ trq_pos_traces <- function(temp)
pos <- temp$a_pos
for(p in 2:nrow(temp))
{
if (abs(as.numeric(pos[p]))>1950 && sign(as.numeric(pos[p]))!=sign(as.numeric(pos[p-1]))) {pos[p-1]="NA"}
if (abs(as.numeric(pos[p]))>1770 && sign(as.numeric(pos[p]))!=sign(as.numeric(pos[p-1]))) {pos[p-1]="NA"}
}
par(mar=c(5, 4, 4, 4) + 0.1)
traces <- plot(x = temp$time, y = temp$a_pos, type = "n", axes=FALSE, xaxs = "i", yaxs = "i", ylim=c(-2047,2048), ylab = "",xlab="")
traces <- plot(x = temp$time, y = temp$a_pos, type = "n", axes=FALSE, xaxs = "i", yaxs = "i", ylim=c(-1800,1800), ylab = "",xlab="")
#shade quadrants
rect(temp$time[1],-1538,temp$time[nrow(temp)],-512, col = "grey95")
rect(temp$time[1],512,temp$time[nrow(temp)],1538, col = "grey95")
rect(temp$time[1],-1350,temp$time[nrow(temp)],-450, col = "grey95")
rect(temp$time[1],450,temp$time[nrow(temp)],1350, col = "grey95")
par(new=TRUE)
plot(x = temp$time, y = pos, type = "l", col="red3", xaxs = "i", yaxs = "i", ylim=c(-2047,2048), ylab = "position[arb.units]",xlab="time [ms]")
plot(x = temp$time, y = pos, type = "l", col="red3", xaxs = "i", yaxs = "i", ylim=c(-1800,1800), ylab = "position[arb.units]",xlab="time [ms]")
par(new=TRUE)
plot(x = temp$time, temp$torque, type = "l", col="blue", xaxs = "i", ylim=c(-500,500), main = paste("Fly Traces", flyname, "Period", i), axes=F, xlab=NA, ylab="")
plot(x = temp$time, temp$torque, type = "l", col="blue", xaxs = "i", ylim=maxtorque, main = paste("Fly Traces", flyname, "Period", i), axes=F, xlab=NA, ylab="")
lines(c(temp$time[1],temp$time[nrow(temp)]),c(0,0),type="l",lty=1,lwd=1, col="black")
axis(4)
mtext("torque [arb_units]", side = 4, line = 3)
......@@ -32,9 +32,9 @@ dytraces <- function(rawdata)
dySeries("a_pos", label = "position", color = "darkred") %>%
dySeries("torque", axis = 'y2', color = "blue") %>%
dyAxis("x", drawGrid = FALSE) %>%
dyAxis("y", label = "Position [arb_units]", valueRange = c(-2048,2048)) %>%
dyAxis("y", label = "Position [arb_units]", valueRange = c(-1800,1800)) %>%
dyOptions(gridLineColor = "lightgrey") %>%
dyAxis("y2", label = "Torque [arb_units]", independentTicks = TRUE, valueRange = c(-700,700), drawGrid = FALSE) %>%
dyAxis("y2", label = "Torque [arb_units]", independentTicks = TRUE, valueRange = maxtorque, drawGrid = FALSE) %>%
dyOptions(includeZero = TRUE) %>%
dyRangeSelector()
......
################## An R-script to read YAML DTS project files, visualize and statistically evaluate data
#####################################################################################################################################
################## R-script to read YAML DTS project files, visualize and statistically evaluate data. Reports in HTML ##############
#####################################################################################################################################
library(ggplot2)
library(tidyr)
......@@ -6,6 +8,7 @@ library(dygraphs)
library(grid)
library(reshape2)
library(dplyr)
library(plyr)
library(gridExtra)
library(yaml)
library(ggsignif)
......@@ -72,6 +75,8 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
NofPeriods = singleflydata[[5]]
sequence <- singleflydata[[6]]
samplerate = as.numeric(as.character(singleflydata[[4]]$sample_rate))
real_sample_rate = as.numeric(as.character(singleflydata[[11]]))
down_sample_rate = as.numeric(as.character(singleflydata[[12]]))
##extract fly meta-data
fly <- singleflydata[[3]]
......@@ -83,13 +88,16 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
##extract the rawdata
rawdata <- singleflydata[[9]]
torquerange = singleflydata[[10]]
#calculate max torque values for axes when plotting
maxtorque = c(-round_any(max(abs(torquerange)), 100, f=ceiling),round_any(max(abs(torquerange)), 100, f=ceiling))
#create/empty plot lists
poshistos <- list()
trqhistos <- list()
#### call RMarkdown for single fly evaluations ################################################
rmarkdown::render('b:/GitHub/DTSevaluations/single_fly.Rmd',
rmarkdown::render(paste(start.wd,"/single_fly.Rmd", sep=""),
output_file = paste(flyname,"descr_anal.html", sep="_"),
output_dir = evaluation.path)
#### end RMarkdown for single fly evaluations ################################################
......@@ -108,20 +116,17 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
########### plot graphs for all experiments #####################
##pool all data except "optomotor" by period
##pool all data by period
pooled.data<-list()
for(i in 1:NofPeriods)
{
period.data<-data.frame()
if(sequence$type[i]!="optomotor")
{
for (l in 1:length(xml_list))
for (l in 1:length(xml_list))
{
period.data <- rbind(period.data, grouped.data[[l]][[i]])
}
}
pooled.data[[i]] <- period.data
} #for number of periods
......@@ -135,7 +140,7 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
trqhistos[[i]] <- ggplot(data=temp, aes_string(temp$torque)) +
geom_histogram(binwidth = 3, fill = sequence$histocolor[i]) +
labs(x="torque [arb units]", y="frequency") +
xlim(-600,600) +
xlim(maxtorque) +
ggtitle(paste("Period", i))
#position
......@@ -144,7 +149,7 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
poshistos[[i]] <- ggplot(data=temp, aes_string(temp$a_pos)) +
geom_histogram(binwidth=10, fill = sequence$histocolor[i]) +
labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
xlim(-1800,1800) +
ggtitle(paste("Period", i))
}
}
......@@ -160,14 +165,14 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
trqhistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$torque)) +
geom_histogram(binwidth=3) +
labs(x="torque [arb units]", y="frequency") +
xlim(-600,600) +
xlim(maxtorque) +
ggtitle("Pooled Torque Histogram")
#position
poshistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$a_pos)) +
geom_histogram(binwidth=10) +
labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
xlim(-1800,1800) +
ggtitle("Pooled Position Histogram")
} else {print("You have selected files with differing metadata")}
......@@ -222,7 +227,7 @@ if(NofGroups>2){}
###### continue for all projects with two groups
#### call RMarkdown for project evaluations ################################################
rmarkdown::render('b:/GitHub/DTSevaluations/project.Rmd',
rmarkdown::render(paste(start.wd,"/project.Rmd", sep=""),
output_file = paste(project.data$name,"html", sep = "."),
output_dir = evaluation.path)
#### end RMarkdown for project evaluations #################################################
......
......@@ -24,8 +24,7 @@ dir.create(evaluation.path, showWarnings = FALSE)
setwd(evaluation.path)
#start evaluating
if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group have the identical experimental design
{
for (l in 1:length(xml_list))
{
xml_name=xml_list[l]
......@@ -48,7 +47,7 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
poshistos <- list()
trqhistos <- list()
#start writing to PDF
#start writing to single-fly PDF
filename = paste(flyname,"descr_anal.pdf", sep="_")
pdf(file=filename, paper="a4r", pointsize=14, width = 0, height = 0)
......@@ -66,11 +65,14 @@ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group ha
grid.table(periods)
##### analyze data from each period separately #######
#add columns for PIs to sequence data
sequence$lambda <- NA
for(i in 1:NofPeriods){
#save colors for later plotting
if(sequence$outcome[i]==0){sequence$color[i]="lightyellow"} else {sequence$color[i]="orange"}
if(sequence$outcome[i]==0){sequence$histocolor[i]="darkgreen"} else {sequence$histocolor[i]="orange"}