speeding up detect_nonwear_clipping for #945

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) {

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)