From 94ad5990de7cab50a9d374ff4a30a031c5db9f91 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 16 Oct 2023 18:12:55 +0200 Subject: [PATCH] speeding up detect_nonwear_clipping for #945 --- R/detect_nonwear_clipping.R | 38 +++++++++++++++------------- tests/testthat/test_detect_nonwear.R | 12 ++++----- 2 files changed, 26 insertions(+), 24 deletions(-) diff --git a/R/detect_nonwear_clipping.R b/R/detect_nonwear_clipping.R index 75083cea8..dd28e7442 100644 --- a/R/detect_nonwear_clipping.R +++ b/R/detect_nonwear_clipping.R @@ -11,6 +11,11 @@ detect_nonwear_clipping = function(data = c(), windowsizes = c(5, 900, 3600), sf CWav = NWav = rep(0, nmin) crit = ((window/window2)/2) + 1 + if (is.numeric(data[,1]) == FALSE) { + data[,1] = as.numeric(data[,1]) + data[,2] = as.numeric(data[,2]) + data[,3] = as.numeric(data[,3]) + } if (nonwear_approach %in% c("2013", "2023")) { # define windows to check: for (h in 1:nmin) { #number of windows @@ -26,8 +31,8 @@ detect_nonwear_clipping = function(data = c(), windowsizes = c(5, 900, 3600), sf hoc1 = 1 hoc2 = window } else if (h >= (nmin - crit)) { - hoc1 = (nmin - crit)*window2 - hoc2 = nmin*window2 #end of data + hoc1 = (nmin - crit) * window2 + hoc2 = nmin * window2 #end of data } else if (h > crit & h < (nmin - crit)) { hoc1 = (((h - 1) * window2) + window2 * 0.5 ) - window * 0.5 hoc2 = (((h - 1) * window2) + window2 * 0.5 ) + window * 0.5 @@ -38,7 +43,7 @@ detect_nonwear_clipping = function(data = c(), windowsizes = c(5, 900, 3600), sf NWflag = h:(h + window/window2 - 1) if (NWflag[length(NWflag)] > nmin) NWflag = NWflag[-which(NWflag > nmin)] # window to check (not aggregated values) - hoc1 = h*window2 - window2 + 1 + hoc1 = h * window2 - window2 + 1 hoc2 = hoc1 + window - 1 if (hoc2 > nrow(data)) { hoc2 = nrow(data) @@ -53,35 +58,32 @@ detect_nonwear_clipping = function(data = c(), windowsizes = c(5, 900, 3600), sf } for (jj in 1:3) { # Clipping - CW[h,jj] = length(which(abs(as.numeric(data[(1 + cliphoc1):cliphoc2,jj])) > clipthres)) - if (length(which(abs(as.numeric(data[(1 + cliphoc1):cliphoc2,jj])) > clipthres*1.5)) > 0) { - CW[h,jj] = window2 # If there is a a value that is more than 150% the dynamic range then ignore entire block. + CW[h, jj] = length(which(abs(data[(1 + cliphoc1):cliphoc2, jj]) > clipthres)) + if (length(which(abs(data[(1 + cliphoc1):cliphoc2,jj]) > clipthres * 1.5)) > 0) { + CW[h, jj] = window2 # If there is a a value that is more than 150% the dynamic range then ignore entire block. } # Non-wear #hoc1 & hoc2 = edges of windows #window is bigger& window2 is smaller one - sdwacc = sd(as.numeric(data[(1 + hoc1):hoc2,jj]), na.rm = TRUE) - maxwacc = max(as.numeric(data[(1 + hoc1):hoc2,jj]), na.rm = TRUE) - minwacc = min(as.numeric(data[(1 + hoc1):hoc2,jj]), na.rm = TRUE) + # (1 + hoc1):hoc2 + indices = seq((1 + hoc1), hoc2, by = ceiling(sf / 5)) + maxwacc = max(data[indices, jj], na.rm = TRUE) + minwacc = min(data[indices, jj], na.rm = TRUE) absrange = abs(maxwacc - minwacc) - if (is.numeric(absrange) == TRUE & is.numeric(sdwacc) == TRUE & is.na(sdwacc) == FALSE) { + if (absrange < racriter) { + sdwacc = sd(data[indices,jj], na.rm = TRUE) if (sdwacc < sdcriter) { - if (absrange < racriter) { - NW[NWflag,jj] = 1 - } + NW[NWflag,jj] = 1 } - } else { - NW[NWflag,jj] = 1 # if sdwacc, maxwacc, or minwacc could not be derived then label as non-wear } } - CW = CW / (window2) #changed 30-1-2012, was window*sf + CW = CW / (window2) if (length(params_rawdata[["rmc.col.wear"]]) == 0) { NWav[h] = (NW[h,1] + NW[h,2] + NW[h,3]) #indicator of non-wear } - CWav[h] = max(c(CW[h,1],CW[h,2],CW[h,3])) #indicator of clipping + CWav[h] = max(c(CW[h, 1], CW[h, 2], CW[h, 3])) #indicator of clipping } } - # In NWav: single 1's surrounded by 2's or 3's --> 2 (so it is considered nonwear) ones = which(NWav == 1) if (length(ones) > 0) { diff --git a/tests/testthat/test_detect_nonwear.R b/tests/testthat/test_detect_nonwear.R index db0678201..f6dd93c4f 100644 --- a/tests/testthat/test_detect_nonwear.R +++ b/tests/testthat/test_detect_nonwear.R @@ -3,27 +3,27 @@ context("detect_nonwear_clipping") test_that("detects non wear time", { skip_on_cran() + # This will produce a 2-day long acc file with a 2-hour block of nonwear # starting at the 5th minute of every day Ndays = 2 - create_test_acc_csv(Nmin = Ndays*1440) + sf = 3 + create_test_acc_csv(Nmin = Ndays * 1440, sf = sf) data = as.matrix(read.csv("123A_testaccfile.csv", skip = 10)) # 2013 algorithm ------ # clipthres to 1.7 to test the clipping detection NWCW = detect_nonwear_clipping(data = data, nonwear_approach = "2013", - sf = 3, clipthres = 1.7) + sf = sf, clipthres = 1.7) CW = NWCW$CWav; NW = NWCW$NWav NW_rle_2013 = rle(NW) CW = sum(NWCW$CWav > 0) - + # 2023 algorithm ------ - NWCW = detect_nonwear_clipping(data = data, nonwear_approach = "2023", sf = 3) + NWCW = detect_nonwear_clipping(data = data, nonwear_approach = "2023", sf = sf) NW = NWCW$NWav NW_rle_2023 = rle(NW) - - # tests ---------------- # Does it find the 2 periods of nonwear? expect_equal(sum(NW_rle_2013$values == 3), 2)