Issue 101 readAxivity: keep one time grid for the whole file so chunked reads join exactly

NEWS.md

@@ -1,3 +1,7 @@
+# Changes in version 1.2.1 (release date:??-11-2026)
+
+- `readAxivity()`: output samples now lie on one time grid for the whole file, anchored at the start time in the file header, and each call returns the grid points up to the start of block `end`. Reading a file in consecutive chunks (as GGIR does) now gives the same samples as a single read. Previously each call started its own grid at its first block and stopped at the last raw sample of the chunk, so the grid of the next chunk was shifted by a fraction of a sample and one sample could be lost at each border. This changes the output slightly: timestamps can shift by less than one sample period, values are interpolated at those times, and the number of rows per call can change by one.
+
 # Changes in version 1.2.0 (release date:28-10-2026)
 
 All changes below concern `readAxivity()`.

R/RcppExports.R

@@ -17,8 +17,8 @@ find_matrix_packet_start <- function(file_path, pattern, chunk_size = 1048576L)
     .Call(`_GGIRread_find_matrix_packet_start`, file_path, pattern, chunk_size)
 }
 
-processAxivityBlocks <- function(readBlocks, calendar, progress, first, time, start, end, numDBlocks, frequency_tol, maxAllowedCorruptBlocks, loadbattery, interpolationType) {
-    .Call(`_GGIRread_processAxivityBlocks`, readBlocks, calendar, progress, first, time, start, end, numDBlocks, frequency_tol, maxAllowedCorruptBlocks, loadbattery, interpolationType)
+processAxivityBlocks <- function(readBlocks, calendar, progress, first, time, start, end, numDBlocks, frequency_tol, maxAllowedCorruptBlocks, loadbattery, interpolationType, fillToNext = FALSE, gridTolerance = 0) {
+    .Call(`_GGIRread_processAxivityBlocks`, readBlocks, calendar, progress, first, time, start, end, numDBlocks, frequency_tol, maxAllowedCorruptBlocks, loadbattery, interpolationType, fillToNext, gridTolerance)
 }
 
 resample <- function(raw, rawTime, time, stop, type = 1L) {

R/readAxivity.R

@@ -329,9 +329,20 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire
   startTimestamp = as.numeric(prevRaw$start)
 
   # allocate memory for results
-  timeRes = seq(startTimestamp, endTimestamp, step)
-  nr = length(timeRes) - 1
-  timeRes = as.vector(timeRes[1:nr])
+  # Output samples lie on one time grid for the whole file, anchored at the
+  # start time in the file header: header$start + k * step. The first sample is
+  # the first grid point at or after the first sample of the start block, and
+  # the last one is the last grid point before the start of the end block
+  # (filled in processAxivityBlocks). Calls for consecutive block ranges
+  # therefore join without lost or repeated samples, and give the same result
+  # as reading those blocks in one call.
+  gridStart = as.numeric(header$start)
+  # tolerance (s) for deciding on which side of a block border a grid point is
+  gridTolerance = 0.001 * step
+  k0 = ceiling((startTimestamp - gridTolerance - gridStart) / step)
+  k1 = ceiling((endTimestamp - gridStart) / step) + 1
+  timeRes = gridStart + (k0:k1) * step
+  nr = length(timeRes)
   # Keep the R connection for platform-independent filenames and cleanup. The
   # native loop requests at most 128 packets at a time, independent of file size.
   readBlocks = function() readBin(fid, raw(), n = blockBytes * 128)
@@ -345,7 +356,8 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire
   }
   result = processAxivityBlocks(readBlocks, calendar, progress, prevRaw, timeRes,
                                 start, end, numDBlocks, frequency_tol,
-                                maxAllowedCorruptBlocks, loadbattery, interpolationType)
+                                maxAllowedCorruptBlocks, loadbattery, interpolationType,
+                                end < numDBlocks, gridTolerance)
   accelRes = result$accel
   temp = result$temp
   battery = result$battery

man/readAxivity.Rd

@@ -123,6 +123,16 @@
   when its magnitude lies between 0.8 and 1.2 g, otherwise with \code{c(0, 0, 1)}.
   Gyroscope channels are filled with zero during imputation.
 
+  Output samples lie on one time grid for the whole file, starting at the
+  start time in the file header and spaced 1 / frequency. A call returns the
+  grid points from the first sample of block \code{start} up to, but not
+  including, the start of block \code{end}. Reading consecutive block ranges
+  in separate calls therefore gives the same samples as reading them in one
+  call. Next to corrupt blocks, a chunk that starts or ends next to a corrupt
+  interval, or that contains only the last valid block of the file, can still
+  differ slightly from a single read, because the frequency check and the
+  duration of the last block are estimated from fewer blocks.
+
   For failed-checksum rows in \code{QClog}, the packet's contents cannot be
   trusted: block IDs are inferred from file positions and timing and frequency
   fields are zero. Such rows have \code{imputed = FALSE}; a separate row records
@@ -132,6 +142,10 @@
   \donttest{
     cwafile  = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread")[1]
     AX3 = readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 1, end = 2)
+    # read in two chunks that join seamlessly
+    header = readAxivity(cwafile, start = 0, end = 0)$header
+    part1 = readAxivity(cwafile, start = 0, end = 50, header = header)
+    part2 = readAxivity(cwafile, start = 50, end = 100, header = header)
   }
 }
 \author{

src/RcppExports.cpp

@@ -63,8 +63,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // processAxivityBlocks
-Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calendar, Rcpp::Nullable<Rcpp::Function> progress, Rcpp::List first, Rcpp::NumericVector time, double start, double end, double numDBlocks, double frequency_tol, int maxAllowedCorruptBlocks, bool loadbattery, int interpolationType);
-RcppExport SEXP _GGIRread_processAxivityBlocks(SEXP readBlocksSEXP, SEXP calendarSEXP, SEXP progressSEXP, SEXP firstSEXP, SEXP timeSEXP, SEXP startSEXP, SEXP endSEXP, SEXP numDBlocksSEXP, SEXP frequency_tolSEXP, SEXP maxAllowedCorruptBlocksSEXP, SEXP loadbatterySEXP, SEXP interpolationTypeSEXP) {
+Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calendar, Rcpp::Nullable<Rcpp::Function> progress, Rcpp::List first, Rcpp::NumericVector time, double start, double end, double numDBlocks, double frequency_tol, int maxAllowedCorruptBlocks, bool loadbattery, int interpolationType, bool fillToNext, double gridTolerance);
+RcppExport SEXP _GGIRread_processAxivityBlocks(SEXP readBlocksSEXP, SEXP calendarSEXP, SEXP progressSEXP, SEXP firstSEXP, SEXP timeSEXP, SEXP startSEXP, SEXP endSEXP, SEXP numDBlocksSEXP, SEXP frequency_tolSEXP, SEXP maxAllowedCorruptBlocksSEXP, SEXP loadbatterySEXP, SEXP interpolationTypeSEXP, SEXP fillToNextSEXP, SEXP gridToleranceSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -80,7 +80,9 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type maxAllowedCorruptBlocks(maxAllowedCorruptBlocksSEXP);
     Rcpp::traits::input_parameter< bool >::type loadbattery(loadbatterySEXP);
     Rcpp::traits::input_parameter< int >::type interpolationType(interpolationTypeSEXP);
-    rcpp_result_gen = Rcpp::wrap(processAxivityBlocks(readBlocks, calendar, progress, first, time, start, end, numDBlocks, frequency_tol, maxAllowedCorruptBlocks, loadbattery, interpolationType));
+    Rcpp::traits::input_parameter< bool >::type fillToNext(fillToNextSEXP);
+    Rcpp::traits::input_parameter< double >::type gridTolerance(gridToleranceSEXP);
+    rcpp_result_gen = Rcpp::wrap(processAxivityBlocks(readBlocks, calendar, progress, first, time, start, end, numDBlocks, frequency_tol, maxAllowedCorruptBlocks, loadbattery, interpolationType, fillToNext, gridTolerance));
     return rcpp_result_gen;
 END_RCPP
 }

src/init.c

@@ -6,15 +6,16 @@
 extern SEXP _GGIRread_AxivityNumUnpack(SEXP);
 extern SEXP _GGIRread_decodeAxivityBlock(SEXP, SEXP, SEXP, SEXP);
 extern SEXP _GGIRread_processAxivityBlocks(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP,
-                                         SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+                                         SEXP, SEXP, SEXP, SEXP, SEXP, SEXP,
+                                         SEXP, SEXP);
 extern SEXP _GGIRread_resample(SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _GGIRread_GENEActivReader(SEXP, SEXP, SEXP, SEXP);
 extern SEXP _GGIRread_find_matrix_packet_start(SEXP, SEXP, SEXP);
 
 static const R_CallMethodDef CallEntries[] = {
     {"_GGIRread_AxivityNumUnpack", (DL_FUNC) &_GGIRread_AxivityNumUnpack, 1},
     {"_GGIRread_decodeAxivityBlock", (DL_FUNC) &_GGIRread_decodeAxivityBlock, 4},
-    {"_GGIRread_processAxivityBlocks", (DL_FUNC) &_GGIRread_processAxivityBlocks, 12},
+    {"_GGIRread_processAxivityBlocks", (DL_FUNC) &_GGIRread_processAxivityBlocks, 14},
     {"_GGIRread_GENEActivReader", (DL_FUNC) &_GGIRread_GENEActivReader, 4},
     {"_GGIRread_resample",  (DL_FUNC) &_GGIRread_resample,  5},
     {"_GGIRread_find_matrix_packet_start",  (DL_FUNC) &_GGIRread_find_matrix_packet_start,  3},

src/processAxivityBlocks.cpp

@@ -66,7 +66,8 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend
                                Rcpp::List first, Rcpp::NumericVector time,
                                double start, double end, double numDBlocks,
                                double frequency_tol, int maxAllowedCorruptBlocks,
-                               bool loadbattery, int interpolationType) {
+                               bool loadbattery, int interpolationType,
+                               bool fillToNext = false, double gridTolerance = 0) {
   using namespace Rcpp;
   if (time.size() == 0 || time.size() > std::numeric_limits<int>::max()) {
     stop("Invalid number of samples in Axivity chunk.");
@@ -259,6 +260,29 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend
         interpolationAnchor[axis] = previous.data(length - 1, axis);
       }
     }
+    if (fillToNext && ii == end && ii < numDBlocks) {
+      // Last iteration of a chunk that ends before the end of the file. The
+      // lookahead packet (current) is not part of the output, but the grid
+      // points before its start belong to this chunk. Fill them as the next
+      // iteration would when the whole file is read: interpolate between the
+      // interpolation anchor and the first sample of the lookahead packet,
+      // with the auxiliary values of that packet. The next chunk then starts
+      // at the first grid point at or after current.start.
+      while (position < nr && time[position] < current.start - gridTolerance) {
+        const double weight = (time[position] - interpolationTime) /
+          (current.start - interpolationTime);
+        for (int axis = 0; axis < axes; ++axis) {
+          const double left = interpolationAnchor[axis];
+          const double right = current.data(0, axis);
+          if (interpolationType == 1) accel(position, axis) = weight * (right - left) + left;
+          else if (interpolationType == 2) accel(position, axis) = std::ceil(weight - 0.5) * (right - left) + left;
+        }
+        temperature[position] = current.temperature;
+        battery[position] = current.battery;
+        light[position] = current.light;
+        ++position;
+      }
+    }
     previous = current;
     if (progress.isNotNull()) Function(progress.get())(position + 1);
   }

tests/testthat/test_processAxivityBlocks.R

@@ -81,9 +81,12 @@ test_that("packet boundaries interpolate through raw samples on every axis", {
     }
     bytes[37] = as.raw(10)
     writeBin(c(bytes[1:1024], unlist(packets)), filename)
-    rawTime = origin + rep((0:2) * 0.375, each = count) +
-      rep((0:(count - 1)) * (0.375 / count), 3)
-    rawValues = do.call(rbind, values[1:3])
+    # Blocks 0-2 are read; the output continues up to the start of block 3,
+    # interpolating towards its first sample.
+    rawTime = c(origin + rep((0:2) * 0.375, each = count) +
+                  rep((0:(count - 1)) * (0.375 / count), 3),
+                origin + 3 * 0.375)
+    rawValues = rbind(do.call(rbind, values[1:3]), values[[4]][1, ])
     for (interpolation in 1:2) {
       result = readAxivity(filename, start = 0, end = 3, configtz = "UTC",
                              interpolationType = interpolation, frequency_tol = 0.5)