Skip to content

Commit

Permalink
tidying up
Browse files Browse the repository at this point in the history
  • Loading branch information
@vincentvanhees
vincentvanhees committed Jul 18, 2024
1 parent 8525afd commit fdeb4ae
Showing 1 changed file with 153 additions and 119 deletions.
272 changes: 153 additions & 119 deletions R/visualreport.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,65 +9,149 @@ visualreport = function(metadatadir = c(),

# Declare functions:
panelplot = function(mdat, ylabels_plot2, Nlevels, selfreport_vars, binary_vars,
BCN, BCC, title = "", dayid) {
BCN, BCC, title = "", NhoursPerRow = NULL) {
window_duration = mdat$timenum[nrow(mdat)] - mdat$timenum[1]

# create vector with color blind friendly colors
mycolors = c("#E69F00","#56B4E9","#009E73","#F0E442",
"#0072B2","#D55E00", "#999999", #"#CC79A7"
"#222255", "black")
mycolors = rep(mycolors, 3)

signalcolor = "grey50"
signalcolor = "black"
mycolors = grDevices::adjustcolor(col = mycolors, alpha.f = 0.6)
# mygreys = rep(c("darkblue", "lightblue"), 20)
mygreys = rep(c("#009E73","#F0E442","#56B4E9", "#D55E00"), 20)
mygreys = grDevices::adjustcolor(col = mygreys, alpha.f = 0.6)


grey_for_lux = gray.colors(n = 100, start = 1, end = 0.1)
myred = grDevices::adjustcolor(col = "#CC79A7", alpha.f = 0.6) #"red"

Ymax = pmax(800, max(mdat$ACC, na.rm = TRUE))
Ymin = 0 #min(mdat$ACC, na.rm = TRUE)

# Prepare time tick points
date = paste0(as.numeric(format(mdat$timestamp, "%d")), " ", format(mdat$timestamp, "%b"))
hour = as.numeric(format(mdat$timestamp, "%H"))
min = as.numeric(format(mdat$timestamp, "%M"))
sec = as.numeric(format(mdat$timestamp, "%S"))
ticks = which(min == 0 & sec == 0 & hour %in% seq(0, 24, by = 1))
atTime = mdat$timestamp[ticks]
datLIM = as.Date(range(mdat$timestamp, na.rm = TRUE), tz = desiredtz)
if (datLIM[1] == datLIM[2]) datLIM = datLIM + c(0, 1)
XLIM = as.POSIXct(c(paste0(datLIM[1], " 00:00:00"),
paste0(datLIM[2], " 12:00:00")), tz = desiredtz)
#----- LUX
if (title == "") {
date = paste0(as.numeric(format(mdat$timestamp[1], "%d")), " ", format(mdat$timestamp[1], "%b"))
title = paste0("Day ", dayid, " | ", weekdays(mdat$timestamp[1]), " | ", date)
}
#----- Acceleration and main non-overlapping behavioural classes
par(mar = c(0, 0, 2, 4.5))
meanAcc = diff(range(mdat$ACC)) / 2
plot(mdat$timestamp, (mdat$ACC / 2) + meanAcc, type = "l",
ylim = c(0, Ymax), xlim = XLIM, col = signalcolor, cex = 0.5, xaxt = 'n', axes = FALSE,
xlab = "", ylab = "", cex.lab = 0.6,
datLIM = as.Date(min(mdat$timestamp, na.rm = TRUE), tz = desiredtz)
XLIM = as.POSIXct(paste0(datLIM[1], " 00:00:00"), tz = desiredtz)
XLIM[2] = XLIM[1] + NhoursPerRow * 3600

# if (datLIM[1] , NhoursPerRow = NULL== datLIM[2]) datLIM = datLIM + c(0, NhoursPerRow/24)
# XLIM = as.POSIXct(c(paste0(datLIM[1], " 00:00:00"),
# paste0(datLIM[2], " 12:00:00")), tz = desiredtz)

par(mar = c(2.5, 0, 1.5, 4.5))
#============================================
# Acceleration and angle signal
# scale 0 - 100, where bottom half is used for angle and top half for acceleration
Ymax = 100
Ymin = 0 #min(mdat$ACC, na.rm = TRUE)
accy = (mdat$ACC / 25) + 50
angy = (((mdat$angle + 90) * 50) / 180)
luxy = ceiling(pmin(mdat$lightpeak + 1, 20000) / 400) * 2
plot(mdat$timestamp, accy, type = "l",
ylim = c(0, Ymax), xlim = XLIM, col = signalcolor,
cex = 0.5, xaxt = 'n', axes = FALSE,
xlab = "", ylab = "", cex.lab = 0.6, cex.axis = 0.7,
main = "", lwd = 0.3, cex.main = 0.9)
lines(mdat$timestamp, angy, type = "l",
col = signalcolor, lwd = 0.3)

lines(mdat$timestamp, (-mdat$ACC / 2) + meanAcc, type = "l", col = signalcolor, cex = 0.5, lwd = 0.3)
text(x = mdat$timestamp[1], y = Ymax * 0.5, labels = "Acceleration",
text(x = mdat$timestamp[1], y = 55, labels = "Acceleration",
pos = 4, cex = 0.7, col = signalcolor, font = 2)
text(x = mdat$timestamp[1], y = 5, labels = "Angle-z",
pos = 4, cex = 0.7, col = signalcolor, font = 2)

title(main = title, adj = 0)

if (NhoursPerRow <= 36) {
cex_axis = 0.5
} else {
cex_axis = 0.5
}

changepoints = function(x) {
Nvalues = length(x)
if (Nvalues > 3) {
changes = which(x[2:Nvalues] != x[1:(Nvalues - 1)])
changes = sort(unique(c(1, changes, changes + 1, Nvalues)))
if (length(changes) == 0) {
changes = 1:Nvalues
}
} else {
changes = 1:Nvalues
}
return(changes)
}

# assign timestamp axis:
mtext(text = hour[ticks], side = 1, line = 0, cex = 0.3, at = atTime)

# with window numbers
windowNow = mdat$window[ticks]
changes = changepoints(windowNow)
mtext(text = windowNow[changes], side = 1, line = 0.75, cex = 0.3,
at = atTime[changes])
# with dates
dateNow = paste0(as.numeric(format(mdat$timestamp[ticks], "%d")), " ",
format(mdat$timestamp[ticks], "%b"))
if (NhoursPerRow <= 36) {
cex_mtext = 0.3
} else {
cex_mtext = 0.25
}

changes = changepoints(dateNow)
mtext(text = dateNow[changes], side = 1, line = 1.5, cex = cex_mtext,
at = atTime[changes])

mtext(text = "hour", side = 1, line = 0, cex = 0.3,
at = atTime[length(atTime)] + 5500, adj = 0)
mtext(text = "window", side = 1, line = 0.75, cex = 0.3,
at = atTime[length(atTime)] + 5500, adj = 0)
mtext(text = "date", side = 1, line = 1.5, cex = 0.3,
at = atTime[length(atTime)] + 5500, adj = 0)


# Add grid of dots
abline(v = atTime, col = "black", lty = "1F", lwd = 0.5)
COL = mycolors[1:Nlevels[1]]
yticks = Ymax * seq(from = 0.1, to = 0.9, length.out = Nlevels[1])

#============================================
# add rectangular blocks to reflect classes
Nlevels = sum(Nlevels) - 1 # -1 because invalidepoch will not be a rectangle
if ("lightpeak" %in% colnames(mdat)) Nlevels = Nlevels + 1
COL = mycolors[1:Nlevels]
yticks = Ymax * seq(from = 0.05, to = 0.95, length.out = Nlevels)
yStepSize = min(diff(yticks)) * 0.5
for (gi in 1:length(yticks)) {
axis(side = 4, at = yticks[gi], lwd = 2,
labels = BCN[gi], col = mycolors[gi], las = 2, cex.axis = 0.5,
line = 0)
}
lev = 1
# Binary classes
for (labi in 1:length(binary_vars)) {
freqtab = table(mdat[,binary_vars[labi]])
if (length(freqtab) > 1 || names(freqtab)[1] == "1") {
t0 = mdat$timestamp[which(diff(c(0, mdat[,binary_vars[labi]])) == 1)]
t1 = mdat$timestamp[which(diff(c(mdat[, binary_vars[labi]], 0)) == -1)]
if (binary_vars[labi] != "invalidepoch") {
y0 = yticks[lev] + yStepSize
y1 = yticks[lev] - yStepSize
col = mycolors[lev]
rect(xleft = t0, xright = t1, ybottom = y0, ytop = y1 , col = col, border = FALSE)
}
}
if (binary_vars[labi] != "invalidepoch") lev = lev + 1
}
# Diary-based variables
for (si in 1:length(selfreport_vars)) {
tempi = which(mdat$selfreport == selfreport_vars[si])
if (length(tempi) > 0) {
A = rep(0, nrow(mdat))
A[tempi] = 1
t0 = mdat$timestamp[which(diff(c(0, A)) == 1)]
t1 = mdat$timestamp[which(diff(c(A, 0)) == -1)]
y0 = yticks[lev] + yStepSize
y1 = yticks[lev] - yStepSize
rect(xleft = t0, xright = t1, ybottom = y0, ytop = y1 , col = mycolors[lev], border = FALSE)
}
lev = lev + 1
}

# Behavioural classes
for (clai in 1:length(BCC)) {
tempi = which(mdat$class_id == BCC[clai])
if (length(tempi) > 0) {
Expand All @@ -77,10 +161,18 @@ visualreport = function(metadatadir = c(),
t1 = mdat$timestamp[which(diff(c(A, 0)) == -1) + 1]
y0 = yticks[lev] - yStepSize * 0.8
y1 = yticks[lev] + yStepSize * 0.8
rect(xleft = t0, xright = t1, ybottom = y0, ytop = y1 , col = COL[clai], border = FALSE) #boarder = FALSE
rect(xleft = t0, xright = t1, ybottom = y0, ytop = y1 , col = mycolors[lev], border = FALSE) #boarder = FALSE
}
lev = lev + 1
}

if ("lightpeak" %in% colnames(mdat)) {
lines(mdat$timestamp, rep(yticks[lev], nrow(mdat)), type = "p", pch = 15,
col = grey_for_lux[luxy], cex = 0.8, lwd = 1)
lev = lev + 1
}

# Highlight invalid epochs as hashed area on top of all rects
if ("invalidepoch" %in% binary_vars) {
freqtab = table(mdat$invalidepoch)
if (length(freqtab) > 1 || names(freqtab)[1] == "1") {
Expand All @@ -92,93 +184,34 @@ visualreport = function(metadatadir = c(),
}
}

# onset/wake lines:
window_edges = which(diff(mdat$SleepPeriodTime) != 0)
if (length(window_edges) > 0) {
abline(v = mdat$timestamp[window_edges], col = "black", lwd = 1)
}

#----- Angle and overlapping classes and self-reported classes:
par(mar = c(2, 0, 0, 4.5), bty = "n")
XLAB = ""
Ymax = 70
angleChange = abs(diff(c(mdat$angle, 0)))
plot(mdat$timestamp, (angleChange / 2), type = "l",
ylim = c(-Ymax, Ymax), xlim = XLIM, col = signalcolor, cex = 0.5, xaxt = 'n', axes = FALSE,
xlab = XLAB, ylab = "", cex.lab = 0.6, lwd = 0.3)

lines(mdat$timestamp, (-angleChange / 2), type = "l",
col = signalcolor, cex = 0.5, lwd = 0.3)
text(x = mdat$timestamp[1], y = Ymax * 0.5, labels = "Angle-z change",
pos = 4, cex = 0.7, col = signalcolor, font = 2)
LUX_scale_hundred = ceiling(pmin(mdat$lightpeak + 1, 20000) / 400) * 2
CL = rep("yellow", 100)
for (ci in 1:100) {
CL[ci] = grDevices::adjustcolor(col = CL[ci], alpha.f = 1/ci)
}
CL = rev(CL)
CL[1:2] = "white"
CL = CL[LUX_scale_hundred]
stepticks = (Ymax*2) / Nlevels[2]
yticks = seq(-Ymax + (Ymax/Nlevels[2]/2), Ymax, by = stepticks)
# Add color labels
tickLabels = c(ylabels_plot2, BCN)
if ("lightpeak" %in% colnames(mdat)) tickLabels = c(tickLabels, "Lux")
tickLabels = tickLabels[which(tickLabels != "invalid")]
for (gi in 1:length(yticks)) {
if (ylabels_plot2[gi] == "invalid") {
col = myred
if (tickLabels[gi] != "Lux") {
col = mycolors[gi]
} else {
col = mygreys[gi]
col = "grey"
}
axis(side = 4, at = yticks[gi], lwd = 2,
labels = ylabels_plot2[gi], col = col, las = 2, cex.axis = 0.5)
}
abline(v = atTime, col = "black", lty = "1F", lwd = 0.5)
# assign timestamp axis:
labTime = paste0(hour[ticks], "H")
axis(side = 1, at = atTime,
labels = labTime, las = 1, cex.axis = 0.5)

buffer = stepticks / 2
lev = 1
for (si in 1:length(selfreport_vars)) {
tempi = which(mdat$selfreport == selfreport_vars[si])
if (length(tempi) > 0) {
A = rep(0, nrow(mdat))
A[tempi] = 1
t0 = mdat$timestamp[which(diff(c(0, A)) == 1)]
t1 = mdat$timestamp[which(diff(c(A, 0)) == -1)]
y0 = yticks[lev] + buffer
y1 = yticks[lev] - buffer
rect(xleft = t0, xright = t1, ybottom = y0, ytop = y1 , col = mygreys[lev], border = FALSE)
}
lev = lev + 1
labels = tickLabels[gi], col = col, las = 2, cex.axis = 0.5,
line = 0)
}

for (labi in 1:length(binary_vars)) {
freqtab = table(mdat[,binary_vars[labi]])
if (length(freqtab) > 1 || names(freqtab)[1] == "1") {
t0 = mdat$timestamp[which(diff(c(0, mdat[,binary_vars[labi]])) == 1)]
t1 = mdat$timestamp[which(diff(c(mdat[, binary_vars[labi]], 0)) == -1)]
if (binary_vars[labi] == "invalidepoch") {
col = myred
rect(xleft = t0, xright = t1, ybottom = -Ymax, ytop = Ymax,
col = col, density = 20, border = FALSE)
# onset/wake lines:
window_edges = which(diff(mdat$SleepPeriodTime) != 0)
if (length(window_edges) > 0) {
abline(v = mdat$timestamp[window_edges], col = "black", lwd = 0.7)
for (wei in 1:length(window_edges)) {
# onset and wake
if (mdat$SleepPeriodTime[window_edges[wei]] == 1) {
toptext = "wake-up"
} else {
y0 = yticks[lev] + buffer
y1 = yticks[lev] - buffer
col = mygreys[lev]
rect(xleft = t0, xright = t1, ybottom = y0, ytop = y1 , col = col, border = FALSE)
toptext = "onset"
}
# print(paste0("y for rect bi: ", y0, " ", y1))

mtext(text = toptext, side = 3, line = 0, cex = 0.5,
at = mdat$timestamp[window_edges[wei]])
}
lev = lev + 1
}

# plot lux
lines(mdat$timestamp, rep(Ymax*1.07, nrow(mdat)), type = "p", pch = 15, col = CL, cex = 0.8, lwd = 1)

# onset/wake lines:
if (length(window_edges) > 0) {
abline(v = mdat$timestamp[window_edges], col = "black", lwd = 1)
}
}

Expand Down Expand Up @@ -231,7 +264,7 @@ visualreport = function(metadatadir = c(),
selfreport_vars = c("nap", "nonwear", "sleeplog")
binary_vars = c("SleepPeriodTime", "sibdetection", "invalidepoch")
Nlevels[2] = length(binary_vars) + length(selfreport_vars)
ylabels_plot2 = c(selfreport_vars, binary_vars)
ylabels_plot2 = c(binary_vars, selfreport_vars)
ylabels_plot2 = gsub(pattern = "invalidepoch", replacement = "invalid", x = ylabels_plot2)
ylabels_plot2 = gsub(pattern = "SleepPeriodTime", replacement = "acc_spt", x = ylabels_plot2)
ylabels_plot2 = gsub(pattern = "sibdetection", replacement = "acc_sib", x = ylabels_plot2)
Expand All @@ -256,12 +289,13 @@ visualreport = function(metadatadir = c(),
mdat$SleepPeriodTime == 0)


NhoursPerRow = 48
midnightsi = which(format(mdat$timestamp, "%H") == "00" &
format(mdat$timestamp, "%M") == "00" &
format(mdat$timestamp, "%S") == "00")
subploti = c(1, midnightsi + 1)
subploti = cbind(subploti,
c(midnightsi + 720, nrow(mdat)))
c(midnightsi + (NhoursPerRow - 24) * 60, nrow(mdat)))

invalid = which(mdat$invalidepoch == 1)
# Skip windows without naps?
Expand All @@ -274,13 +308,13 @@ visualreport = function(metadatadir = c(),
subploti[which(subploti[,2] > nrow(mdat)), 2] = nrow(mdat)

NdaysPerPage = 7
par(mfrow = c(NdaysPerPage * 2, 1), mgp = c(2, 0.8, 0), omi = c(0, 0, 0, 0), bty = "n")
par(mfrow = c(NdaysPerPage, 1), mgp = c(2, 0.8, 0), omi = c(0, 0, 0, 0), bty = "n")
if (nrow(subploti) > 0) {
for (ani in 1:nrow(subploti)) {

panelplot(mdat[(subploti[ani, 1] + 1):subploti[ani, 2], ],
ylabels_plot2, Nlevels, selfreport_vars, binary_vars,
BCN, BCC, title = "", dayid = ani)
BCN, BCC, title = "", NhoursPerRow = NhoursPerRow)
}
}
dev.off()
Expand Down

0 comments on commit fdeb4ae

Please sign in to comment.