From 398894f1443105a2705be665a70231ac4fea56d3 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Sat, 17 Oct 2026 11:54:56 +0200 Subject: [PATCH 01/14] Fix readAxivity single-block reads at end of file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Use the final block’s recorded sampling frequency when no preceding block is available to estimate timing, and initialise interpolation state so its samples are returned. Add AX3 and AX6 regression tests for both interpolation modes and corrupt starting blocks that leave only the final valid block. --- NEWS.md | 4 +++ R/readAxivity.R | 24 ++++++++++++---- tests/testthat/test_readAxivity.R | 48 +++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index 83637cd..cc1999a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# Development version + +- Fix `readAxivity()` failing when a requested chunk contains only the final data block. Use the block's recorded sampling frequency when no preceding block is available to estimate timing. + # Changes in version 1.0.11 (release date:18-09-2026) - Added GENEActiv light calibration needed for device model GENEActiv 1.2 and above. Contributed by ActivInsights Ltd, the manufacturer of GENEActiv devices. diff --git a/R/readAxivity.R b/R/readAxivity.R index acf83f4..7ec523c 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -526,7 +526,15 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire # Process the last block in the file, if necessary # Calculate pseudo time for the "next" block - newTimes = (prevRaw$start - prevStart) / prevLength * prevRaw$length + prevRaw$start + if (rawPos == 1) { + # Only the final block was requested, so there is no preceding block + # from which to estimate its sample spacing. Use its recorded frequency. + prevStart = prevRaw$start + newTimes = prevStart + prevRaw$length / prevRaw$frequency + raw = prevRaw + } else { + newTimes = (prevRaw$start - prevStart) / prevLength * prevRaw$length + prevRaw$start + } prevLength = prevRaw$length # fill vector rawTime and matrix rawAccel for resampling @@ -574,15 +582,19 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire # fill vector rawTime and matrix rawAccel for resampling rawLast = prevLength + 1 rawTime[2:(rawLast+1)] = seq(prevStart, raw$start, length.out = rawLast) # rawTime[rawLast+1] will be ignored by resampling alg + } - if (rawPos == 1) { - rawAccel[1,] = (prevRaw$data[1,]) - rawTime[1] = prevStart - 0.00001 - rawPos = 2 - } + if (rawPos == 1) { + rawAccel[1,] = prevRaw$data[1,] + rawTime[1] = prevStart - 0.00001 + rawPos = 2 } frequency_observed = rawLast / (rawTime[rawLast] - rawTime[1]) + if (ii == numDBlocks && start == numDBlocks - 1) { + # No inter-block timing observation is available for this fallback. + frequency_observed = prevRaw$frequency + } #------------------------------------------------------------ # Check block integrity: # The following code checks whether any of the following conditions are met: diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index 8137dea..6f4db2f 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -1,5 +1,53 @@ library(GGIRread) context("reading Axivity .cwa data") + +test_that("readAxivity returns the final block when it is the only requested block", { + for (device in c("ax3", "ax6")) { + cwafile = system.file(paste0("testfiles/", device, "_testfile.cwa"), + package = "GGIRread") + n = (file.size(cwafile) - 1024) / 512 + header = readAxivity(cwafile, desiredtz = "UTC")$header + for (interpolationType in c(1, 2)) { + result = readAxivity(cwafile, start = n - 1, end = n, + desiredtz = "UTC", header = header, + interpolationType = interpolationType, + loadbattery = TRUE, frequency_tol = 0) + expect_equal(nrow(result$data), header$blockLength) + expect_equal(ncol(result$data), if (device == "ax3") 7L else 10L) + expect_equal(diff(result$data$time), + rep(1 / header$frequency, header$blockLength - 1), + tolerance = 1e-6) + expect_true(all(result$data$battery > 0)) + expect_null(result$QClog) + # First and last accelerometer samples from the fixture's final packet. + expected = if (device == "ax3") { + rbind(c(0.59375, 0.328125, 0.59375), c(-0.0625, -0.84375, 0.265625)) + } else { + rbind(c(104, 2007, 15), c(98, 2010, 23)) / 2048 + } + expect_equal(unname(as.matrix(result$data[c(1, header$blockLength), + c("x", "y", "z")])), expected, + tolerance = 1e-6) + } + + # The same path is reached when a corrupt start packet is skipped, + # leaving just the final valid packet in the requested chunk. + bytes = readBin(cwafile, "raw", n = file.size(cwafile)) + offset = 1024 + (n - 2) * 512 + 1 + bytes[offset] = as.raw(bitwXor(as.integer(bytes[offset]), 1L)) + corruptfile = tempfile(fileext = ".cwa") + writeBin(bytes, corruptfile) + expect_warning(result <- readAxivity(corruptfile, start = n - 2, end = n, + desiredtz = "UTC"), + "Skipping corrupt start block") + reference = readAxivity(cwafile, start = n - 1, end = n, desiredtz = "UTC") + expect_equal(result$data, reference$data) + expect_equal(nrow(result$QClog), 1L) + expect_false(result$QClog$checksum_pass) + unlink(corruptfile) + } +}) + test_that("readAxivity reads data from AX3 file correctly", { cwafile = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread")[1] AX3 = readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 1, end = 4) From e26df7b974fb065858acc322fee7774981507759 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Sat, 17 Oct 2026 12:38:40 +0200 Subject: [PATCH 02/14] Speed up readAxivity with buffered C++ processing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Move packet decoding, timestamp caching, integrity checks, imputation, and resampling into C++, retaining the R interface and timezone conversion. Add regression tests and a benchmark script excluded from package builds. Real-recording validation produced identical outputs with approximately 4.5× faster reading. --- .Rbuildignore | 3 +- NEWS.md | 2 + R/RcppExports.R | 8 + R/readAxivity.R | 445 ++------------------- src/RcppExports.cpp | 36 ++ src/decodeAxivityBlock.cpp | 140 +++++++ src/init.c | 5 + src/processAxivityBlocks.cpp | 232 +++++++++++ tests/testthat/test_decodeAxivityBlock.R | 159 ++++++++ tests/testthat/test_processAxivityBlocks.R | 110 +++++ tools/benchmark-readAxivity.R | 81 ++++ 11 files changed, 815 insertions(+), 406 deletions(-) create mode 100644 src/decodeAxivityBlock.cpp create mode 100644 src/processAxivityBlocks.cpp create mode 100644 tests/testthat/test_decodeAxivityBlock.R create mode 100644 tests/testthat/test_processAxivityBlocks.R create mode 100644 tools/benchmark-readAxivity.R diff --git a/.Rbuildignore b/.Rbuildignore index 315ca92..bc3532e 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -3,4 +3,5 @@ ^\.github prepareNewRelease.R LICENSE -codecov.yml \ No newline at end of file +codecov.yml +^tools/benchmark-readAxivity\.R$ diff --git a/NEWS.md b/NEWS.md index cc1999a..a93f23d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # Development version +- Speed up `readAxivity()` with C++ packet decoding and a buffered native processing loop for timestamp caching, integrity checks, imputation and resampling. Retain the R interface and timezone conversion. + - Fix `readAxivity()` failing when a requested chunk contains only the final data block. Use the block's recorded sampling frequency when no preceding block is available to estimate timing. # Changes in version 1.0.11 (release date:18-09-2026) diff --git a/R/RcppExports.R b/R/RcppExports.R index a789d39..e5300dd 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,10 +9,18 @@ GENEActivReader <- function(filename, start = 0L, end = 0L, progress_bar = FALSE .Call(`_GGIRread_GENEActivReader`, filename, start, end, progress_bar) } +decodeAxivityBlock <- function(block, complete, parameters, loadbattery) { + .Call(`_GGIRread_decodeAxivityBlock`, block, complete, parameters, loadbattery) +} + 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) +} + resample <- function(raw, rawTime, time, stop, type = 1L) { .Call(`_GGIRread_resample`, raw, rawTime, time, stop, type) } diff --git a/R/readAxivity.R b/R/readAxivity.R index 7ec523c..84bfcc8 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -52,213 +52,34 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire invisible(list(start = start, struc = struc)) } - readDataBlock = function(fid, complete = TRUE, struc = list(0,0L,0), parameters = NULL){ - # Read one block of data and return list with following elements - # frequency is frequency recorded in this block - # start is start time in numeric form. To create string representation - # it is necessary to use - # as.POSIXct(start, origin = "1970-01-01", tz=desiredtz) - # temperature is temperature for the block - # battery is battery charge for the block - # light is light sensor measurement for the block - # length is number of observations in the block - # data is matrix with three columns "x", "y", and "z" - # matrix data is presented if complete == TRUE only. - # - if (!is.null(parameters)) { - accelScaleCode = parameters$accelScaleCode - accelScale = parameters$accelScale - Naxes = parameters$Naxes - frequency_data = parameters$frequency_data - format = parameters$format - } - block = readBin(fid, raw(), n=blockBytes) - if (length(block) < blockBytes) { - return(NULL) - } - - # sampling rate in one of file format U8 at offset 24 - samplerate_dynrange = readBin(block[25], integer(), size = 1, signed = FALSE) - - if (samplerate_dynrange != 0) { # Very old files that have zero at offset 24 don't have a checksum - checksum = sum(readBin(block, n = 256, - integer(), - size = 2, - signed = FALSE, - endian = "little")) - checksum = checksum %% 65536 # 65536 = 2^16; the checksum is calculated as a 16-bit integer - if (checksum != 0) { - # Checksum doesn't match. This means some bits in this block got corrupted. - # We don't know which, so we can't trust this block. We skip it, and impute it later. - rawdata_list = list( - struc = struc, - parameters = parameters, - checksum_pass = FALSE - ) - return(invisible(rawdata_list)) - } - } - - idstr = readChar(block, 2, useBytes = TRUE) - if (idstr != "AX") { - stop("Packet header is incorrect. First two characters must be AX.") - } - - packetLength = readBin(block[3:4], integer(), size = 2, signed = FALSE, endian = "little") - if (packetLength != 508) { - stop("Packet length is incorrect, should always be 508.") - } - - # offset 4: if the top bit set, this contains a 15-bit fraction of a second for the timestamp - tsOffset = readBin(block[5:6], integer(), size = 2, signed = FALSE, endian = "little") - - # offset 10: sequence ID - blockID = readBin(block[11:14], integer(), size = 4, endian = "little") - - # read data for timestamp u32 at offset 14 - timeStamp = readBin(block[15:18], integer(), size = 4, endian = "little") # the "signed" flag of readBin only works when reading 1 or 2 bytes - - # Get light u16 at offset 18 - offset18 = readBin(block[19:20], integer(), size = 2, signed = FALSE, endian = "little") - light = bitwAnd(offset18, 0x03ffL) - - # Read and recalculate temperature, lower 10 bits of u16 at offset 20. - # Formula for the temperature is specified at - # https://github.com/digitalinteraction/openmovement/blob/545564d3bf45fc19914de1ad1523ed86538cfe5e/Docs/ax3/cwa.h#L102 - # Also see the following discussion: - # https://github.com/digitalinteraction/openmovement/issues/11#issuecomment-1622278513 - temperature = bitwAnd(readBin(block[21:22], integer(), size = 2, signed = FALSE, endian = "little"), 0x03ffL) * 75.0 / 256.0 - 50; - if (loadbattery == TRUE) { - # Read and recalculate battery charge u8 in offset 23 - # https://github.com/digitalinteraction/openmovement/blob/master/Docs/ax3/ax3-auxiliary.md#battery-voltage - # Battery is sampled as a 10-bit ADC value, but only the middle 8 bits are stored (the lowest bit is lost, and the highest bit is always 1). - # So to restore the ADC value, double the packed value and add 512. - # Then voltage = ADC_value * 6 / 1024 - battery = 3.0 * (readBin(block[24], integer(), size = 1, signed = FALSE) / 256.0 + 1.0); - } else { - battery = 0 - } - - # offset 25, per documentation: - # "top nibble: number of axes, 3=Axyz, 6=Gxyz/Axyz, 9=Gxyz/Axyz/Mxyz; - # bottom nibble: packing format" (2 means unpacked, 0 packed). - offset25 = readBin(block[26], integer(), size = 1, signed = FALSE) - packed = (bitwAnd(offset25,15) == 0) - - # offset 26 has a int16 (not uint16) value. - # It's the "relative sample index from the start of the buffer where the whole-second timestamp is valid" - offset26 = readBin(block[27:28], integer(), size = 2, endian = "little") - - # number of observations in block U16 at offset 28 - # blockLength is expected to be 40 for AX6, 80 or 120 for AX3. - # Note that if AX6 is configured to only collect accelerometer data - # this will look as if it is a AX3 - blockLength = readBin(block[29:30], integer(), size = 2, signed = FALSE, endian = "little") - - if (is.null(parameters)) { - accelScaleCode = bitwShiftR(offset18, 13) - accelScale = 1 / (2^(8 + accelScaleCode)) - # top nibble of offset25 is the number of axes - Naxes = bitwShiftR(offset25, 4) - } - - # auxiliary variables - shift = 0 - fractional = 0 - - if (samplerate_dynrange == 0) { - # Very old files have zero at offset 24 and frequency at offset 26 - frequency_data = offset26 - } else { - # value at offset 26 is index of measurement with whole number of seconds - shift = offset26 - if (is.null(parameters)) { - frequency_data = round( 3200 / bitwShiftL(1, 15 - bitwAnd(samplerate_dynrange, 15))) - # If the top bit of tsOffset is set, then timestamp offset was artificially - # modified for backwards-compatibility ... therefore undo this... - if (bitwAnd(tsOffset, 0x8000L) != 0) { - format = 1 - } else { - format = 2 - } - } - if (format == 1) { - # Need to undo backwards-compatible shim: - # Take into account how many whole samples the fractional part - # of timestamp accounts for: - # relativeOffset = fifoLength - # - (short)(((unsigned long)timeFractional * AccelFrequency()) >> 16); - # nearest whole sample - # whole-sec | /fifo-pos@time - # | |/ - # [0][1][2][3][4][5][6][7][8][9] - # use 15-bits as 16-bit fractional time - fractional = bitwShiftL(bitwAnd(tsOffset, 0x7fffL), 1); - - # frequency is truncated to int in firmware - shift = shift + bitwShiftR((fractional * frequency_data), 16); - } - } - - # Read data if necessary - if (complete) { - - if (packed) { - # Read 4 bytes for three measurements - packedData = readBin(block[31:510], integer(), size = 4, n = blockLength, endian = "little") - # Unpack data - data = AxivityNumUnpack(packedData) - } else { - # Read unpacked data - xyz = readBin(block[31:510], integer(), size = 2, n = blockLength * Naxes, endian = "little") - data = matrix(xyz, ncol = Naxes, byrow = T) - } - - # Set names and Normalize accelerations - if (Naxes == 3) { - colnames(data) = c("x", "y", "z") - data[,c("x", "y", "z")] = data[,c("x", "y", "z")] * accelScale #/ 256 - } else { - gyroRangeCode = floor(offset18 / 1024) %% 8 - gyroRange = 8000 / (2^gyroRangeCode) - colnames(data) = c("gx", "gy", "gz", "x", "y", "z") - data[,c("gx", "gy", "gz")] = (data[,c("gx", "gy", "gz")] / 2^15) * gyroRange - data[,c("x", "y", "z")] = data[,c("x", "y", "z")] * accelScale - } - } - if (is.null(parameters)) { - parameters = list(accelScaleCode = accelScaleCode, - accelScale = accelScale, - Naxes = Naxes, - frequency_data = frequency_data, - format = format) - } - - if (blockID == 0) { - # force timestampDecoder to extract full timestamp - # this could for example happen when curious participant plugs - # AX device into computer + readDataBlock = function(fid, complete = TRUE, struc = list(0,0L,0), parameters = NULL) { + # Decode one packet in C++; keep timezone conversion and timestamp state in R. + block = readBin(fid, raw(), n = blockBytes) + if (length(block) < blockBytes) return(NULL) + decoded = decodeAxivityBlock(block, complete, parameters, loadbattery) + if (!decoded$checksum_pass) { + return(invisible(list(struc = struc, parameters = parameters, + checksum_pass = FALSE))) + } + if (decoded$blockID == 0) { + # Force full calendar conversion after a device sequence reset. struc[[1]] = 0 } - tsDeco = timestampDecoder(timeStamp, fractional, -shift / frequency_data, struc, configtz) - start = tsDeco$start - struc = tsDeco$struc + tsDeco = timestampDecoder(decoded$timeStamp, decoded$fractional, + decoded$shift, struc, configtz) rawdata_list = list( - frequency = frequency_data, - start = start, - temperature = temperature, - battery = battery, - light = light, - length = blockLength, - struc = struc, - parameters = parameters, + frequency = decoded$frequency, + start = tsDeco$start, + temperature = decoded$temperature, + battery = decoded$battery, + light = decoded$light, + length = decoded$length, + struc = tsDeco$struc, + parameters = decoded$parameters, checksum_pass = TRUE, - blockID = blockID + blockID = decoded$blockID ) - if (complete) { - rawdata_list$data = data - } - + if (complete) rawdata_list$data = decoded$data return(invisible(rawdata_list)) } @@ -496,212 +317,26 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire timeRes = seq(startTimestamp, endTimestamp, step) nr = length(timeRes) - 1 timeRes = as.vector(timeRes[1:nr]) - temp = vector(mode = "double", nr) - battery = vector(mode = "double", nr) - light = vector(mode = "double", nr) - - # Create progress bar if it is necessary + # 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) + calendar = function(coded) { + timestampDecoder(coded, 0, 0, list(0, 0L, 0), configtz)$start + } + progress = NULL if (progressBar) { pb = txtProgressBar(1, nr, style = 3) + progress = function(pos) setTxtProgressBar(pb, pos) } - - # Allocate enough space for the expected number of samples in a block, plus an extra ones needed for resampling. - # Don't rely on the type of device to determine the dimensionality of the data - # because AX6 can be configured to only collect accelerometer data. - if (prevRaw$parameters$Naxes == 3) { # AX3, or AX6 configured to only collect accelerometer data - accelRes = matrix(0, nrow = nr, ncol = 3, dimnames = list(NULL, c("x", "y", "z"))) - rawAccel = matrix(0, nrow = blockLength + 1, ncol = 3) - } else { # AX6 configured to collect gyroscope data - accelRes = matrix(0, nrow = nr, ncol = 6, dimnames = list(NULL, c("gx", "gy", "gz", "x", "y", "z"))) - rawAccel = matrix(0, nrow = blockLength + 1, ncol = 6) - } - rawTime = vector(mode = "numeric", blockLength + 2) - rawPos = 1 - - consecutiveCorrupt = 0 - - # Read the data - for (ii in (start+1):end) { - if (ii == numDBlocks) { - # Process the last block in the file, if necessary - - # Calculate pseudo time for the "next" block - if (rawPos == 1) { - # Only the final block was requested, so there is no preceding block - # from which to estimate its sample spacing. Use its recorded frequency. - prevStart = prevRaw$start - newTimes = prevStart + prevRaw$length / prevRaw$frequency - raw = prevRaw - } else { - newTimes = (prevRaw$start - prevStart) / prevLength * prevRaw$length + prevRaw$start - } - prevLength = prevRaw$length - - # fill vector rawTime and matrix rawAccel for resampling - rawLast = prevLength + 1 - rawTime[2:(rawLast+1)] = seq(prevStart, newTimes, length.out = rawLast) # rawTime[rawLast+1] will be ignored by resampling alg - } - else { # read a new block - raw = readDataBlock(fid, struc = struc, parameters = prevRaw$parameters) - if (!raw$checksum_pass) { - # If the checksum doesn't match, we can't trust any of this block's data, - # so we have to completely skip the block. - # Depending on the nature of the faulty block, the data for the time period it represented - # will probably get imputed later, once we encounter a block with a valid checksum. - QClog = rbind(QClog, data.frame(checksum_pass = FALSE, - blockID_current = ii, # the actual block ID wasn't read b/c the block is corrupt - blockID_next = ii + 1, - start = 0, - end = 0, - blockLengthSeconds = 0, - frequency_blockheader = 0, - frequency_observed = 0, - imputed = FALSE)) - warning("Skipping corrupt block #", ii) - - consecutiveCorrupt = consecutiveCorrupt + 1 - if (consecutiveCorrupt > maxAllowedCorruptBlocks) { - stop("Error reading file. ", maxAllowedCorruptBlocks+1, " consecutive blocks are corrupt.") - } - - next - } - - consecutiveCorrupt = 0 - - if (is.null(raw)) { - # this shouldn't happen - stop(paste0("\nreadAxivity encountered unexpected empty block at #", ii)) - } - - # Save start and length of the previous block - prevStart = prevRaw$start - prevLength = prevRaw$length - struc = raw$struc - - # fill vector rawTime and matrix rawAccel for resampling - rawLast = prevLength + 1 - rawTime[2:(rawLast+1)] = seq(prevStart, raw$start, length.out = rawLast) # rawTime[rawLast+1] will be ignored by resampling alg - } - - if (rawPos == 1) { - rawAccel[1,] = prevRaw$data[1,] - rawTime[1] = prevStart - 0.00001 - rawPos = 2 - } - - frequency_observed = rawLast / (rawTime[rawLast] - rawTime[1]) - if (ii == numDBlocks && start == numDBlocks - 1) { - # No inter-block timing observation is available for this fallback. - frequency_observed = prevRaw$frequency - } - #------------------------------------------------------------ - # Check block integrity: - # The following code checks whether any of the following conditions are met: - # - blockID is not zero and not consecutive from previous blockID - # - checksum_pass is FALSE - # - observed and expected sampling frequency differ by a fraction larger - # than frequency_tol - # If yes, then we consider the block faulty - # and impute the acceleration and if applicable gyroscope values - # We log this event in output object QClog, which will allow the user to - # decide on alternative imputation strategies. - - impute = FALSE - doQClog = FALSE - frequency_bias = abs(frequency_observed - prevRaw$frequency) / prevRaw$frequency - if ((ii < numDBlocks && - raw$blockID != 0 && - raw$blockID - prevRaw$blockID != 1) || - frequency_bias > frequency_tol) { - # Log and impute this event - doQClog = TRUE - impute = TRUE - - # Prepare imputation with last recorded value, from the last non-faulty block, - # normalized to vector 1 g - imputedValues = rawAccel[1, 1:3] - VectorG = sqrt(sum(imputedValues^2)) - if (VectorG > 0.8 & VectorG < 1.2) { - # only trust vector as proxy for orientation if it is between 0.8 and 1.2 - imputedValues = imputedValues / VectorG - } else { - imputedValues = c(0, 0, 1) - } - } else { - # integrity check passes, so use data - rawAccel[2:rawLast,] = prevRaw$data - # If frequency bias is larger than 0.05 then still log it even though it is - # not imputed - if (frequency_bias > 0.05) { - doQClog = TRUE - } - } - if (doQClog == TRUE) { - # Note: This is always a description of the previous block - QClog = rbind(QClog, data.frame(checksum_pass = TRUE, - blockID_current = prevRaw$blockID, - blockID_next = raw$blockID, - start = prevRaw$start, - end = raw$start, - blockLengthSeconds = raw$start - prevRaw$start, - frequency_blockheader = prevRaw$frequency, - frequency_observed = frequency_observed, - imputed = impute)) - } - - ########################################################################### - # resampling of measurements - last = pos + 200; - if (last > nr) last = nr - if (rawTime[rawLast] > timeRes[last]) { - # there has been a time jump - # so, time jump needs to be adjusted for in last index - timejump = rawTime[rawLast] - timeRes[last] - positions2add = floor(timejump * prevRaw$frequency) - last = last + positions2add - if (last > nr) last = nr - } - - if (impute == FALSE) { - tmp = resample(rawAccel, rawTime, timeRes[pos:last], rawLast, type = interpolationType) - } else { - # Impute the data because integrity check did not pass - if (last - pos > prevRaw$frequency * 259200) { # 3600 * 24 * 5 = 259200 - # Error if time gap is very large to avoid filling up memory - stop(paste0("\nreadAxivity encountered a time gap in the file of ", - round((last - pos) / (3600 * 24) / prevRaw$frequency, digits = 2), " days")) - } - - # Figure out the number of points to impute, numImp. - # numImp should be such that timeRes[numImp + pos + 1] is the last point in timeRes[] < rawTime[rawLast] - numImp = length(which(timeRes[pos:last]= pos) { - accelRes[pos:last,] = tmp - light[pos:last] = prevRaw$light - temp[pos:last] = prevRaw$temperature - battery[pos:last] = prevRaw$battery - } - # Remove all rawdata except for the last - rawTime[1] = timeRes[last] - rawAccel[1,] = accelRes[last, ] - rawPos = 2 - # Now current becomes previous - prevRaw = raw - pos = last + 1 - # Refresh progress bar if it is necessary - if (progressBar) { - setTxtProgressBar(pb, pos) - } - } + result = processAxivityBlocks(readBlocks, calendar, progress, prevRaw, timeRes, + start, end, numDBlocks, frequency_tol, + maxAllowedCorruptBlocks, loadbattery, interpolationType) + accelRes = result$accel + temp = result$temp + battery = result$battery + light = result$light + last = result$last + QClog = rbind(QClog, result$QClog) #=============================================================================== # If the user asked for more data than the length of the recording, # there will be 0s at the end of the result lists; get rid of them. diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 531c163..32f44b1 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -35,6 +35,20 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// decodeAxivityBlock +Rcpp::List decodeAxivityBlock(Rcpp::RawVector block, bool complete, Rcpp::Nullable parameters, bool loadbattery); +RcppExport SEXP _GGIRread_decodeAxivityBlock(SEXP blockSEXP, SEXP completeSEXP, SEXP parametersSEXP, SEXP loadbatterySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::RawVector >::type block(blockSEXP); + Rcpp::traits::input_parameter< bool >::type complete(completeSEXP); + Rcpp::traits::input_parameter< Rcpp::Nullable >::type parameters(parametersSEXP); + Rcpp::traits::input_parameter< bool >::type loadbattery(loadbatterySEXP); + rcpp_result_gen = Rcpp::wrap(decodeAxivityBlock(block, complete, parameters, loadbattery)); + return rcpp_result_gen; +END_RCPP +} // find_matrix_packet_start Rcpp::IntegerVector find_matrix_packet_start(std::string file_path, Rcpp::RawVector pattern, int chunk_size); RcppExport SEXP _GGIRread_find_matrix_packet_start(SEXP file_pathSEXP, SEXP patternSEXP, SEXP chunk_sizeSEXP) { @@ -48,6 +62,28 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// processAxivityBlocks +Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calendar, Rcpp::Nullable 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) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::Function >::type readBlocks(readBlocksSEXP); + Rcpp::traits::input_parameter< Rcpp::Function >::type calendar(calendarSEXP); + Rcpp::traits::input_parameter< Rcpp::Nullable >::type progress(progressSEXP); + Rcpp::traits::input_parameter< Rcpp::List >::type first(firstSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type time(timeSEXP); + Rcpp::traits::input_parameter< double >::type start(startSEXP); + Rcpp::traits::input_parameter< double >::type end(endSEXP); + Rcpp::traits::input_parameter< double >::type numDBlocks(numDBlocksSEXP); + Rcpp::traits::input_parameter< double >::type frequency_tol(frequency_tolSEXP); + 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)); + return rcpp_result_gen; +END_RCPP +} // resample NumericMatrix resample(NumericMatrix raw, NumericVector rawTime, NumericVector time, int stop, int type); RcppExport SEXP _GGIRread_resample(SEXP rawSEXP, SEXP rawTimeSEXP, SEXP timeSEXP, SEXP stopSEXP, SEXP typeSEXP) { diff --git a/src/decodeAxivityBlock.cpp b/src/decodeAxivityBlock.cpp new file mode 100644 index 0000000..2693f89 --- /dev/null +++ b/src/decodeAxivityBlock.cpp @@ -0,0 +1,140 @@ +#include +#include +#include + +namespace { +// Decode explicitly: CWA packets are little endian and need not be aligned. +unsigned int u16(const Rcpp::RawVector& b, int offset) { + return static_cast(b[offset]) | + (static_cast(b[offset + 1]) << 8); +} + +int s16(const Rcpp::RawVector& b, int offset) { + const unsigned int value = u16(b, offset); + return value >= 0x8000U ? static_cast(value) - 0x10000 : value; +} + +std::uint32_t u32(const Rcpp::RawVector& b, int offset) { + return static_cast(u16(b, offset)) | + (static_cast(u16(b, offset + 2)) << 16); +} + +int signed10(std::uint32_t value) { + const int sample = static_cast(value & 0x3ffU); + return sample >= 0x200 ? sample - 0x400 : sample; +} +} + +// Packet decoding only. Calendar/timezone conversion and the inter-packet +// integrity, imputation and interpolation rules remain in readAxivity(). +// CWA layout: https://github.com/digitalinteraction/openmovement/blob/master/Docs/ax3/cwa.h +// [[Rcpp::export]] +Rcpp::List decodeAxivityBlock(Rcpp::RawVector block, bool complete, + Rcpp::Nullable parameters, + bool loadbattery) { + using namespace Rcpp; + if (block.size() != 512) stop("Axivity data blocks must contain 512 bytes."); + + const int samplerate = block[24]; + if (samplerate != 0) { + unsigned int checksum = 0; + for (int offset = 0; offset < 512; offset += 2) checksum += u16(block, offset); + if ((checksum & 0xffffU) != 0) return List::create(_["checksum_pass"] = false); + } + if (block[0] != 'A' || block[1] != 'X') { + stop("Packet header is incorrect. First two characters must be AX."); + } + if (u16(block, 2) != 508) { + stop("Packet length is incorrect, should always be 508."); + } + + const unsigned int tsOffset = u16(block, 4); + const unsigned int lightScale = u16(block, 18); + const int axesPacking = block[25]; + const bool packed = (axesPacking & 15) == 0; + const int offset26 = s16(block, 26); + const int count = u16(block, 28); + const bool reuse = parameters.isNotNull(); + List params; + int axes, scaleCode, format; + double accelScale, frequency; + if (reuse) { + params = List(parameters.get()); + axes = as(params["Naxes"]); + scaleCode = as(params["accelScaleCode"]); + accelScale = as(params["accelScale"]); + frequency = as(params["frequency_data"]); + format = as(params["format"]); + } else { + axes = axesPacking >> 4; + scaleCode = lightScale >> 13; + accelScale = 1.0 / std::pow(2.0, 8 + scaleCode); + // nearbyint, like R's round(), rounds ties to even. + frequency = samplerate == 0 ? offset26 : + std::nearbyint(3200.0 / (1U << (15 - (samplerate & 15)))); + format = samplerate == 0 ? 0 : ((tsOffset & 0x8000U) != 0 ? 1 : 2); + params = List::create(_["accelScaleCode"] = scaleCode, + _["accelScale"] = accelScale, + _["Naxes"] = axes, + _["frequency_data"] = frequency, + _["format"] = static_cast(format)); + } + // Old-format packets carry their frequency in each packet, even when the + // other decoding parameters are reused from the preceding packet. + if (samplerate == 0) frequency = offset26; + if (axes != 3 && axes != 6) stop("Unsupported number of Axivity axes."); + if (packed && axes != 3) stop("Packed Axivity data must have three axes."); + const int capacity = packed ? 120 : 480 / (2 * axes); + if (count < 1 || count > capacity) stop("Invalid Axivity sample count for packet format."); + if (!std::isfinite(frequency) || frequency <= 0) stop("Invalid Axivity sampling frequency."); + + int fractional = 0; + double shift = 0; + if (samplerate != 0) { + shift = offset26; + if (format == 1) { + // Undo the firmware's backwards-compatible timestamp shim: the lower + // 15 bits hold a fraction of a second, and offset26 includes its samples. + fractional = (tsOffset & 0x7fffU) << 1; + shift += std::floor(fractional * frequency / 65536.0); + } + } + // Temperature uses the low ten ADC bits. Battery stores only the middle + // eight ADC bits; restoring the implicit high bit gives 3 * (byte/256 + 1). + List result = List::create( + _["frequency"] = frequency, + _["temperature"] = (u16(block, 20) & 0x3ffU) * 75.0 / 256.0 - 50.0, + _["battery"] = loadbattery ? 3.0 * (block[23] / 256.0 + 1.0) : 0.0, + _["light"] = static_cast(lightScale & 0x3ffU), + _["length"] = count, + _["parameters"] = params, + _["checksum_pass"] = true, + _["blockID"] = static_cast(u32(block, 10)), + _["timeStamp"] = static_cast(u32(block, 14)), + _["fractional"] = fractional, + _["shift"] = -shift / frequency); + if (complete) { + NumericMatrix data(count, axes); + const double gyroRange = 8000.0 / (1U << ((lightScale >> 10) & 7U)); + for (int row = 0; row < count; ++row) { + if (packed) { + const std::uint32_t value = u32(block, 30 + 4 * row); + const int multiplier = 1 << (value >> 30); + for (int axis = 0; axis < 3; ++axis) { + // Multiplication avoids left-shifting a negative signed integer. + data(row, axis) = signed10(value >> (10 * axis)) * multiplier * accelScale; + } + } else { + for (int axis = 0; axis < axes; ++axis) { + const int sample = s16(block, 30 + 2 * (row * axes + axis)); + data(row, axis) = axes == 6 && axis < 3 ? + (sample / 32768.0) * gyroRange : sample * accelScale; + } + } + } + colnames(data) = axes == 3 ? CharacterVector::create("x", "y", "z") : + CharacterVector::create("gx", "gy", "gz", "x", "y", "z"); + result["data"] = data; + } + return result; +} diff --git a/src/init.c b/src/init.c index febf426..5b34752 100644 --- a/src/init.c +++ b/src/init.c @@ -4,12 +4,17 @@ #include 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); 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_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}, diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp new file mode 100644 index 0000000..605d5c9 --- /dev/null +++ b/src/processAxivityBlocks.cpp @@ -0,0 +1,232 @@ +#include +#include +#include +#include +#include + +Rcpp::List decodeAxivityBlock(Rcpp::RawVector, bool, + Rcpp::Nullable, bool); + +namespace { +struct Block { + double start, frequency, temperature, battery, light; + int length, id; + Rcpp::NumericMatrix data; + Block(const Rcpp::List& x, double timestamp) : + start(timestamp), frequency(Rcpp::as(x["frequency"])), + temperature(Rcpp::as(x["temperature"])), + battery(Rcpp::as(x["battery"])), light(Rcpp::as(x["light"])), + length(Rcpp::as(x["length"])), id(Rcpp::as(x["blockID"])), + data(x["data"]) {} +}; + +struct QC { + bool checksum, imputed; + int current, next; + double start, end, duration, expected, observed; +}; + +SEXP qualityLog(const std::vector& rows) { + if (rows.empty()) return R_NilValue; + const R_xlen_t n = rows.size(); + Rcpp::LogicalVector checksum(n), imputed(n); + Rcpp::IntegerVector current(n), next(n); + Rcpp::NumericVector start(n), end(n), duration(n), expected(n), observed(n); + bool hasCorrupt = false; + for (R_xlen_t i = 0; i < n; ++i) { + const QC& row = rows[i]; + hasCorrupt = hasCorrupt || !row.checksum; + checksum[i] = row.checksum; imputed[i] = row.imputed; + current[i] = row.current; next[i] = row.next; + start[i] = row.start; end[i] = row.end; duration[i] = row.duration; + expected[i] = row.expected; observed[i] = row.observed; + } + // The R checksum-failure rows use ii + 1 (a double) for blockID_next; + // otherwise both IDs come directly from integer packet headers. + Rcpp::RObject nextIDs = hasCorrupt ? + Rcpp::RObject(Rcpp::as(next)) : Rcpp::RObject(next); + return Rcpp::DataFrame::create( + Rcpp::_["checksum_pass"] = checksum, + Rcpp::_["blockID_current"] = current, Rcpp::_["blockID_next"] = nextIDs, + Rcpp::_["start"] = start, Rcpp::_["end"] = end, + Rcpp::_["blockLengthSeconds"] = duration, + Rcpp::_["frequency_blockheader"] = expected, + Rcpp::_["frequency_observed"] = observed, Rcpp::_["imputed"] = imputed); +} +} + +// Process one requested chunk. R retains the file connection (including Unicode +// paths and connection cleanup), endpoint selection and calendar/timezone rules. +// readBlocks supplies a bounded byte buffer; calendar is called only when the +// cached whole-minute timestamp changes or the device sequence resets. +// [[Rcpp::export]] +Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calendar, + Rcpp::Nullable progress, + Rcpp::List first, Rcpp::NumericVector time, + double start, double end, double numDBlocks, + double frequency_tol, int maxAllowedCorruptBlocks, + bool loadbattery, int interpolationType) { + using namespace Rcpp; + if (time.size() == 0 || time.size() > std::numeric_limits::max()) { + stop("Invalid number of samples in Axivity chunk."); + } + const int nr = time.size(); + Block previous(first, as(first["start"])); + List parameters = first["parameters"]; + const int axes = as(parameters["Naxes"]); + List state = first["struc"]; + double timestamp = as(state[0]); + int oldSeconds = as(state[1]); + unsigned int oldMinute = as(state[2]); + NumericMatrix accel(nr, axes); + colnames(accel) = axes == 3 ? CharacterVector::create("x", "y", "z") : + CharacterVector::create("gx", "gy", "gz", "x", "y", "z"); + NumericVector temperature(nr), battery(nr), light(nr); + std::vector qc; + Function sum = Environment::base_env()["sum"]; + Function power = Environment::base_env()["^"]; + Function squareRoot = Environment::base_env()["sqrt"]; + std::vector anchor(axes, 0.0), rawTime; + double anchorTime = 0, prevStart = 0; + int prevLength = 0, position = 0, consecutiveCorrupt = 0; + bool firstIteration = true; + RawVector buffer(0), packet(512); + R_xlen_t offset = 0; + + for (double ii = start + 1; ii <= end; ++ii) { + if (static_cast(ii - start) % 128 == 0) checkUserInterrupt(); + Block current = previous; + double nextTime; + if (ii == numDBlocks) { + // Preserve the established final-packet extrapolation, including the + // nominal-frequency fallback for a chunk containing only that packet. + if (firstIteration) { + prevStart = previous.start; + nextTime = prevStart + previous.length / previous.frequency; + } else { + nextTime = (previous.start - prevStart) / prevLength * previous.length + previous.start; + } + prevLength = previous.length; + } else { + if (offset + 512 > buffer.size()) { + buffer = as(readBlocks()); + offset = 0; + } + if (offset + 512 > buffer.size()) { + stop("\nreadAxivity encountered unexpected empty block at #%.0f", ii); + } + std::copy(buffer.begin() + offset, buffer.begin() + offset + 512, packet.begin()); + offset += 512; + List decoded = decodeAxivityBlock(packet, true, parameters, loadbattery); + if (!as(decoded["checksum_pass"])) { + qc.push_back(QC{false, false, static_cast(ii), static_cast(ii + 1), + 0, 0, 0, 0, 0}); + warning("Skipping corrupt block #%.0f", ii); + if (++consecutiveCorrupt > maxAllowedCorruptBlocks) { + stop("Error reading file. %d consecutive blocks are corrupt.", maxAllowedCorruptBlocks + 1); + } + continue; + } + consecutiveCorrupt = 0; + const int coded = as(decoded["timeStamp"]); + const unsigned int minute = static_cast(coded) >> 6; + const int seconds = static_cast(coded) & 0x3fU; + if (as(decoded["blockID"]) == 0 || minute != oldMinute || timestamp == 0) { + timestamp = as(calendar(coded)); + } else { + if (seconds < oldSeconds) oldSeconds -= 60; + timestamp += seconds - oldSeconds; + } + oldMinute = minute; + oldSeconds = seconds; + // Match timestampDecoder's arithmetic order, including the fraction. + const double blockStart = timestamp + as(decoded["fractional"]) / 65536.0 + + as(decoded["shift"]); + current = Block(decoded, blockStart); + prevStart = previous.start; + prevLength = previous.length; + nextTime = current.start; + } + + // Equivalent to seq(prevStart, nextTime, length.out = prevLength + 1). + // The endpoint is a time anchor, not a sample supplied to interpolation. + if (!std::isfinite(prevStart)) stop("'from' must be a finite number"); + if (!std::isfinite(nextTime)) stop("'to' must be a finite number"); + const int length = prevLength; + rawTime.resize(length + 2); + const double spacing = (nextTime - prevStart) / length; + rawTime[1] = prevStart; + for (int i = 1; i < length; ++i) rawTime[i + 1] = prevStart + i * spacing; + rawTime[length + 1] = nextTime; + if (firstIteration) { + for (int axis = 0; axis < axes; ++axis) anchor[axis] = previous.data(0, axis); + anchorTime = prevStart - 0.00001; + firstIteration = false; + } + rawTime[0] = anchorTime; + double observed = (length + 1.0) / (rawTime[length] - anchorTime); + if (ii == numDBlocks && start == numDBlocks - 1) observed = previous.frequency; + const double bias = std::abs(observed - previous.frequency) / previous.frequency; + const bool impute = (ii < numDBlocks && current.id != 0 && + static_cast(current.id) - previous.id != 1) || bias > frequency_tol; + if (impute || bias > 0.05) { + qc.push_back(QC{true, impute, previous.id, current.id, previous.start, current.start, + current.start - previous.start, previous.frequency, observed}); + } + int last = std::min(static_cast(position) + 200, nr - 1); + if (rawTime[length] > time[last]) { + const double add = std::floor((rawTime[length] - time[last]) * previous.frequency); + last = static_cast(std::min(static_cast(last) + add, static_cast(nr - 1))); + } + if (impute && last - position > previous.frequency * 259200) { + Function round = Environment::base_env()["round"]; + Function paste0 = Environment::base_env()["paste0"]; + stop(as(paste0("\nreadAxivity encountered a time gap in the file of ", + round((last - position) / (3600.0 * 24) / previous.frequency, 2), " days"))); + } + double imputed[3] = {0, 0, 1}; + if (impute) { + // Retain R's platform-specific summation for the three-value norm. This + // rare callback avoids changing repeated normalisation by one ULP. + NumericVector orientation(3); + for (int axis = 0; axis < 3; ++axis) { + orientation[axis] = anchor[axis]; + } + const double norm = as(squareRoot(sum(power(orientation, 2.0)))); + if (norm > 0.8 && norm < 1.2) { + for (int axis = 0; axis < 3; ++axis) imputed[axis] = anchor[axis] / norm; + } + } + int source = 1; + while (position <= last && (impute ? time[position] < rawTime[length] : + time[position] <= rawTime[length])) { + if (impute) { + for (int axis = 0; axis < axes; ++axis) accel(position, axis) = axis < 3 ? imputed[axis] : 0; + } else { + while (source < length && rawTime[source] < time[position]) ++source; + const double weight = (time[position] - rawTime[source - 1]) / + (rawTime[source] - rawTime[source - 1]); + for (int axis = 0; axis < axes; ++axis) { + const double left = source == 1 ? anchor[axis] : previous.data(source - 2, axis); + const double right = previous.data(source - 1, 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] = previous.temperature; + battery[position] = previous.battery; + light[position] = previous.light; + ++position; + } + if (position == 0) stop("No samples could be resampled from the Axivity block."); + anchorTime = time[position - 1]; + for (int axis = 0; axis < axes; ++axis) anchor[axis] = accel(position - 1, axis); + previous = current; + if (progress.isNotNull()) Function(progress.get())(position + 1); + } + if (position == 0) stop("No samples could be resampled from the Axivity chunk."); + RObject log = qualityLog(qc); + return List::create(_["accel"] = accel, _["temp"] = temperature, + _["battery"] = battery, _["light"] = light, + _["last"] = position, _["QClog"] = log); +} diff --git a/tests/testthat/test_decodeAxivityBlock.R b/tests/testthat/test_decodeAxivityBlock.R new file mode 100644 index 0000000..655c765 --- /dev/null +++ b/tests/testthat/test_decodeAxivityBlock.R @@ -0,0 +1,159 @@ +library(GGIRread) + +axivity_checksum = function(packet) { + packet[511:512] = as.raw(0) + checksum = (-sum(readBin(packet, integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packet[511:512] = writeBin(as.integer(checksum), raw(), size = 2, + endian = "little") + packet +} + +axivity_packet = function(device = "ax3") { + filename = system.file(paste0("testfiles/", device, "_testfile.cwa"), + package = "GGIRread") + readBin(filename, "raw", n = 1536)[1025:1536] +} + +test_that("native decoder matches fixture samples and auxiliary measurements", { + for (device in c("ax3", "ax6")) { + filename = system.file(paste0("testfiles/", device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(filename, "raw", n = file.size(filename)) + n = (length(bytes) - 1024) / 512 + parameters = NULL + for (i in seq_len(n)) { + packet = bytes[1024 + (i - 1) * 512 + seq_len(512)] + decoded = GGIRread:::decodeAxivityBlock(packet, TRUE, parameters, TRUE) + count = readBin(packet[29:30], integer(), size = 2, signed = FALSE, + endian = "little") + lightScale = readBin(packet[19:20], integer(), size = 2, signed = FALSE, + endian = "little") + if (device == "ax3") { + packed = readBin(packet[31:510], integer(), n = count, size = 4, + endian = "little") + expected = GGIRread:::AxivityNumUnpack(packed) / 2^(8 + bitwShiftR(lightScale, 13)) + } else { + expected = matrix(readBin(packet[31:510], integer(), n = count * 6, + size = 2, endian = "little"), ncol = 6, byrow = TRUE) + expected[,1:3] = expected[,1:3] / 32768 * (8000 / 2^(floor(lightScale / 1024) %% 8)) + expected[,4:6] = expected[,4:6] / 2^(8 + bitwShiftR(lightScale, 13)) + } + expect_identical(unname(decoded$data), unname(expected)) + expect_equal(decoded$light, bitwAnd(lightScale, 1023L)) + expect_equal(decoded$temperature, + bitwAnd(readBin(packet[21:22], integer(), size = 2, + signed = FALSE, endian = "little"), 1023L) * 75 / 256 - 50) + expect_equal(decoded$battery, 3 * (as.integer(packet[24]) / 256 + 1)) + parameters = decoded$parameters + } + metadata = GGIRread:::decodeAxivityBlock(packet, FALSE, parameters, FALSE) + expect_false("data" %in% names(metadata)) + expect_equal(metadata$battery, 0) + expect_identical(metadata$parameters, parameters) + } +}) + +test_that("packed decoding handles signed ten-bit values and all exponents", { + packet = axivity_packet() + # x = -512, y = -1, z = 511. Build words as bytes to include high-bit values. + for (exponent in 0:3) { + word = 512 + 1023 * 2^10 + 511 * 2^20 + exponent * 2^30 + packet[31:34] = as.raw(floor(word / 256^(0:3)) %% 256) + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE) + expect_equal(unname(decoded$data[1,]), + c(-512, -1, 511) * 2^exponent * decoded$parameters$accelScale) + } +}) + +test_that("unpacked decoding supports accelerometer-only packets and signed limits", { + packet = axivity_packet("ax6") + packet[26] = as.raw(0x32) # Three axes, signed 16-bit unpacked samples. + packet[29:30] = writeBin(80L, raw(), size = 2, endian = "little") + values = rep(c(-32768L, -1L, 32767L), 80) + packet[31:510] = writeBin(values, raw(), size = 2, endian = "little") + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE) + expect_equal(unname(decoded$data), matrix(values, ncol = 3, byrow = TRUE) * + decoded$parameters$accelScale) + expect_identical(colnames(decoded$data), c("x", "y", "z")) +}) + +test_that("native decoding preserves timestamp offsets and parameter reuse", { + packet = axivity_packet() + packet[5:6] = writeBin(as.integer(0x8000 + 12345), raw(), size = 2, endian = "little") + packet[27:28] = writeBin(-23L, raw(), size = 2, endian = "little") + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, FALSE, NULL, FALSE) + expect_equal(decoded$fractional, 24690) + expect_equal(decoded$shift, + -(-23 + floor(24690 * decoded$frequency / 65536)) / decoded$frequency) + expect_identical(GGIRread:::decodeAxivityBlock(packet, FALSE, decoded$parameters, FALSE), decoded) + + packet[5:6] = as.raw(0) + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, FALSE, NULL, FALSE) + expect_equal(decoded$fractional, 0) + expect_equal(decoded$shift, 23 / decoded$frequency) + + # Older packets encode frequency directly, without checksum or time shift. + packet[25] = as.raw(0) + packet[27:28] = writeBin(100L, raw(), size = 2, endian = "little") + decoded = GGIRread:::decodeAxivityBlock(packet, FALSE, NULL, FALSE) + expect_true(decoded$checksum_pass) + expect_equal(decoded$frequency, 100) + expect_equal(decoded$shift, 0) + expect_identical(GGIRread:::decodeAxivityBlock(packet, FALSE, decoded$parameters, FALSE), decoded) + packet[27:28] = writeBin(200L, raw(), size = 2, endian = "little") + expect_equal(GGIRread:::decodeAxivityBlock(packet, FALSE, decoded$parameters, FALSE)$frequency, 200) +}) + +test_that("native decoding applies all accelerometer and gyroscope scale codes", { + for (device in c("ax3", "ax6")) { + original = axivity_packet(device) + lightScale = readBin(original[19:20], integer(), size = 2, + signed = FALSE, endian = "little") + base = GGIRread:::decodeAxivityBlock(original, TRUE, NULL, FALSE) + originalGyroCode = floor(lightScale / 1024) %% 8 + for (code in 0:7) { + packet = original + replacement = lightScale %% 1024 + code * 1024 + code * 8192 + packet[19:20] = writeBin(as.integer(replacement), raw(), size = 2, endian = "little") + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE) + expect_equal(decoded$data[,c("x", "y", "z")], + base$data[,c("x", "y", "z")] * + (1 / 2^(8 + code)) / base$parameters$accelScale) + if (device == "ax6") { + expect_equal(decoded$data[,c("gx", "gy", "gz")], + base$data[,c("gx", "gy", "gz")] * 2^(originalGyroCode - code)) + } + # Accelerometer scale is cached, whereas gyro range is packet-specific. + reused = GGIRread:::decodeAxivityBlock(packet, TRUE, base$parameters, FALSE) + expect_identical(reused$data[,c("x", "y", "z")], base$data[,c("x", "y", "z")]) + } + } +}) + +test_that("native decoder validates packets before accessing sample payloads", { + packet = axivity_packet() + expect_error(GGIRread:::decodeAxivityBlock(packet[-1], TRUE, NULL, FALSE), "512 bytes") + packet[1] = as.raw(0) + expect_identical(GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE), + list(checksum_pass = FALSE)) + expect_error(GGIRread:::decodeAxivityBlock(axivity_checksum(packet), TRUE, NULL, FALSE), + "Packet header is incorrect") + packet = axivity_packet() + packet[3:4] = as.raw(0) + expect_error(GGIRread:::decodeAxivityBlock(axivity_checksum(packet), TRUE, NULL, FALSE), + "Packet length is incorrect") + packet = axivity_packet() + packet[29:30] = writeBin(121L, raw(), size = 2, endian = "little") + expect_error(GGIRread:::decodeAxivityBlock(axivity_checksum(packet), TRUE, NULL, FALSE), + "Invalid Axivity sample count") + packet = axivity_packet("ax6") + packet[29:30] = writeBin(41L, raw(), size = 2, endian = "little") + expect_error(GGIRread:::decodeAxivityBlock(axivity_checksum(packet), TRUE, NULL, FALSE), + "Invalid Axivity sample count") +}) diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R new file mode 100644 index 0000000..578d0fb --- /dev/null +++ b/tests/testthat/test_processAxivityBlocks.R @@ -0,0 +1,110 @@ +library(GGIRread) + +test_that("native loop crosses input buffers and calendar days at a low sample rate", { + # Two days at 1 Hz exercise date changes without allocating a multi-day + # 100 Hz output. Payloads are constant; timestamps and checksums are genuine. + fixture = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") + bytes = readBin(fixture, "raw", n = 1536) + packet = bytes[1025:1536] + packet[25] = as.raw(3) # round(3200 / 2^(15 - 3)) = 1 Hz + packet[5:6] = as.raw(0) + packet[27:28] = as.raw(0) + packet[31:510] = rep(writeBin(256L, raw(), size = 4, endian = "little"), 120) + blocks = 1442L + origin = as.POSIXct("2020-01-01", tz = "UTC") + times = as.POSIXlt(origin + (seq_len(blocks) - 1) * 120, tz = "UTC") + coded = bitwShiftL(as.integer(times$year - 100), 26) + + bitwShiftL(as.integer(times$mon + 1), 22) + bitwShiftL(as.integer(times$mday), 17) + + bitwShiftL(as.integer(times$hour), 12) + bitwShiftL(as.integer(times$min), 6) + packets = matrix(rep(packet, blocks), nrow = 512) + packets[11:14,] = writeBin(seq_len(blocks) - 1L, raw(), size = 4, endian = "little") + packets[15:18,] = writeBin(as.integer(coded), raw(), size = 4, endian = "little") + packets[511:512,] = as.raw(0) + for (i in seq_len(blocks)) { + checksum = (-sum(readBin(packets[,i], integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packets[511:512,i] = writeBin(as.integer(checksum), raw(), size = 2, endian = "little") + } + bytes[37] = as.raw(3) + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + writeBin(c(bytes[1:1024], as.vector(packets)), filename) + result = readAxivity(filename, start = 0, end = blocks - 1, configtz = "UTC") + samples = (blocks - 1) * 120 + expect_equal(nrow(result$data), samples) + expect_identical(result$data$time, as.numeric(origin) + seq_len(samples) - 1) + scale = 1 / 2^(8 + bitwShiftR(as.integer(packet[20]), 5)) + expect_true(all(result$data$x == 256 * scale)) + expect_true(all(result$data$y == 0)) + expect_true(all(result$data$z == 0)) + expect_null(result$QClog) +}) + +test_that("native buffering preserves corrupt-packet logs and limits", { + fixture = system.file("testfiles/ax6_testfile.cwa", package = "GGIRread") + bytes = readBin(fixture, "raw", n = file.size(fixture)) + # After the initial packet, each input buffer contains another 128 packets. + # Corrupt packets straddle a refill boundary. + for (block in 128:129) { + offset = 1024 + block * 512 + 1 + bytes[offset] = as.raw(bitwXor(as.integer(bytes[offset]), 1L)) + } + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + writeBin(bytes, filename) + expect_warning(result <- readAxivity(filename, start = 0, end = 132, configtz = "UTC"), + "Skipping corrupt block #128") + failed = result$QClog[!result$QClog$checksum_pass,] + expect_identical(failed$blockID_current, 128:129) + expect_identical(failed$blockID_next, c(129, 130)) + expect_true(any(result$QClog$blockID_current == 127 & result$QClog$imputed)) + expect_warning(expect_error(readAxivity(filename, start = 0, end = 132, + configtz = "UTC", maxAllowedCorruptBlocks = 1), + "2 consecutive blocks are corrupt"), + "Skipping corrupt block") +}) + +test_that("repeated imputation retains the R normalisation arithmetic", { + filename = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") + result = readAxivity(filename, start = 13, end = 16, + configtz = "UTC", frequency_tol = 0) + expected = c(62, 12, 11) / 64 + for (block in 1:3) { + expected = expected / sqrt(sum(expected^2)) + first = which(result$data$time >= result$QClog$start[block])[1] + expect_identical(as.numeric(result$data[first, c("x", "y", "z")]), expected) + } +}) + +test_that("a corrupt tail without a usable endpoint fails explicitly", { + filename = system.file("testfiles/ax3_testfile_corrupt_blocks_0_13_14_142_143_144.cwa", + package = "GGIRread") + # The former R loop left 'last' undefined after skipping every lookahead + # packet. Do not return a successful but empty data frame in this case. + expect_warning(expect_error(readAxivity(filename, start = 141, end = 145, + configtz = "UTC"), + "No samples could be resampled from the Axivity chunk"), + "Skipping corrupt") +}) + +test_that("a device sequence reset refreshes the native timestamp cache", { + fixture = system.file("testfiles/ax6_testfile.cwa", package = "GGIRread") + bytes = readBin(fixture, "raw", n = file.size(fixture)) + blocks = (length(bytes) - 1024) / 512 + for (block in 130:(blocks - 1)) { + positions = 1024 + block * 512 + seq_len(512) + packet = bytes[positions] + packet[11:14] = writeBin(as.integer(block - 130), raw(), size = 4, endian = "little") + packet[511:512] = as.raw(0) + checksum = (-sum(readBin(packet, integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packet[511:512] = writeBin(as.integer(checksum), raw(), size = 2, endian = "little") + bytes[positions] = packet + } + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + writeBin(bytes, filename) + original = readAxivity(fixture, start = 0, end = 200, configtz = "Europe/Berlin") + reset = readAxivity(filename, start = 0, end = 200, configtz = "Europe/Berlin") + expect_identical(reset, original) +}) diff --git a/tools/benchmark-readAxivity.R b/tools/benchmark-readAxivity.R new file mode 100644 index 0000000..271a179 --- /dev/null +++ b/tools/benchmark-readAxivity.R @@ -0,0 +1,81 @@ +# Source this file after loading the development package, then run: +# benchmark_readAxivity("/path/to/previous/readAxivity.R") +# Named paths can compare several stages in the same run: +# benchmark_readAxivity(c(original = "original.R", decoder = "decoder.R")) +# The reference is a saved copy of R/readAxivity.R from before the optimisation. +# Synthetic data repeat sensor payloads, but have continuous timestamps, +# consecutive sequence IDs and valid checksums. They measure throughput, not +# compatibility with the full variation of real recordings. + +make_axivity_benchmark_file = function(device, blocks = 10000L) { + stopifnot(device %in% c("ax3", "ax6"), blocks > 0) + fixture = system.file(paste0("testfiles/", device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(fixture, "raw", n = 1536) + packet = bytes[1025:1536] + samples = readBin(packet[29:30], integer(), size = 2, signed = FALSE, endian = "little") + rate = 3200 / 2^(15 - bitwAnd(as.integer(packet[25]), 15L)) + times = as.POSIXlt(as.POSIXct("2020-01-01", tz = "UTC") + + (seq_len(blocks) - 1) * samples / rate, tz = "UTC") + coded = bitwShiftL(as.integer(times$year - 100), 26) + + bitwShiftL(as.integer(times$mon + 1), 22) + bitwShiftL(as.integer(times$mday), 17) + + bitwShiftL(as.integer(times$hour), 12) + bitwShiftL(as.integer(times$min), 6) + as.integer(times$sec) + packets = matrix(rep(packet, blocks), nrow = 512) + packets[5:6,] = as.raw(0) + packets[11:14,] = writeBin(as.integer(seq_len(blocks) - 1), raw(), size = 4, endian = "little") + packets[15:18,] = writeBin(as.integer(coded), raw(), size = 4, endian = "little") + offsets = -as.integer(round((times$sec - floor(times$sec)) * rate)) + packets[27:28,] = writeBin(offsets, raw(), size = 2, endian = "little") + packets[511:512,] = as.raw(0) + for (i in seq_len(blocks)) { + checksum = (-sum(readBin(packets[,i], integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packets[511:512,i] = writeBin(as.integer(checksum), raw(), size = 2, endian = "little") + } + filename = tempfile(fileext = ".cwa") + writeBin(c(bytes[1:1024], as.vector(packets)), filename) + filename +} + +benchmark_readAxivity = function(reference_file, blocks = c(1000L, 10000L), repetitions = 5L) { + readers = lapply(reference_file, function(filename) { + reference = new.env(parent = asNamespace("GGIRread")) + sys.source(filename, reference) + reference$readAxivity + }) + if (is.null(names(readers))) { + names(readers) = if (length(readers) == 1) "previous" else paste0("reference", seq_along(readers)) + } + readers$native = GGIRread::readAxivity + results = list() + for (device in c("ax3", "ax6")) { + for (n in blocks) { + filename = make_axivity_benchmark_file(device, n) + tryCatch({ + # Warm up both paths, and compare the entire output before timing. + outputs = lapply(readers, function(reader) { + reader(filename, start = 0, end = n, configtz = "UTC") + }) + stopifnot(all(vapply(outputs, identical, logical(1), outputs$native))) + rm(outputs) + timings = matrix(NA_real_, repetitions, length(readers), dimnames = list(NULL, names(readers))) + for (i in seq_len(repetitions)) { + # Alternate order to reduce systematic cache/order effects. + order = if (i %% 2) seq_along(readers) else rev(seq_along(readers)) + for (j in order) { + gc() + timings[i,j] = system.time(invisible(readers[[j]]( + filename, start = 0, end = n, configtz = "UTC")))[["elapsed"]] + } + } + elapsed = apply(timings, 2, median) + row = data.frame(device = device, blocks = n) + for (name in names(elapsed)) row[[paste0(name, "_seconds")]] = elapsed[[name]] + row$speedup = unname(elapsed[1] / elapsed["native"]) + print(row) + results[[length(results) + 1L]] = row + }, finally = unlink(filename)) + } + } + invisible(do.call(rbind, results)) +} From d3b8881317d0ef2c387c42000cd1417d546922bc Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Sat, 17 Oct 2026 12:46:34 +0200 Subject: [PATCH 03/14] Update NEWS.md --- NEWS.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index a93f23d..a10868f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,8 @@ -# Development version +# Changes in version 1.0.12 (release date:??-??-2026) - Speed up `readAxivity()` with C++ packet decoding and a buffered native processing loop for timestamp caching, integrity checks, imputation and resampling. Retain the R interface and timezone conversion. -- Fix `readAxivity()` failing when a requested chunk contains only the final data block. Use the block's recorded sampling frequency when no preceding block is available to estimate timing. +- Fix `readAxivity()` failing when a requested chunk contains only the final data block. Use the block's recorded sampling frequency when no preceding block is available to estimate timing. #98 # Changes in version 1.0.11 (release date:18-09-2026) From 25706708d35ef1ba1b510f18df52624ec80f5ad5 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Oct 2026 17:27:41 +0200 Subject: [PATCH 04/14] Address PR100 feedback: Fix Axivity compatibility, final-block timing and memory retention MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Use explicit Rcpp conversions for C++11/14 compatibility - Correct the final packet’s time axis in reads ending at EOF - Release native buffers before assembling trimmed output - Add AX3/AX6 regression tests and update reference expectations - Update NEWS and benchmark parity checks for corrected EOF behavior Validation: full test suite passes; compilation checked under C++11/14/17. --- NEWS.md | 4 ++++ R/readAxivity.R | 3 +++ src/processAxivityBlocks.cpp | 10 +++++---- tests/testthat/test_readAxivity.R | 37 +++++++++++++++++++++++++++---- tools/benchmark-readAxivity.R | 5 +++-- 5 files changed, 49 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index a10868f..28902b2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,10 @@ - Fix `readAxivity()` failing when a requested chunk contains only the final data block. Use the block's recorded sampling frequency when no preceding block is available to estimate timing. #98 +- Fix the final block's time axis in `readAxivity()` reads containing multiple blocks up to the end of a file. This corrects an existing timing error and can change the final samples and output row count. + +- Keep the native Axivity processing loop compatible with C++11 and C++14, and release its full output buffers before assembling trimmed results to reduce memory retention. + # Changes in version 1.0.11 (release date:18-09-2026) - Added GENEActiv light calibration needed for device model GENEActiv 1.2 and above. Contributed by ActivInsights Ltd, the manufacturer of GENEActiv devices. diff --git a/R/readAxivity.R b/R/readAxivity.R index 84bfcc8..8e9db25 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -337,6 +337,9 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire light = result$light last = result$last QClog = rbind(QClog, result$QClog) + # Release the references to the full native buffers before trimming them. + # Otherwise result keeps those buffers alive while the output is assembled. + rm(result) #=============================================================================== # If the user asked for more data than the length of the recording, # there will be 0s at the end of the result lists; get rid of them. diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index 605d5c9..d0caf68 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -17,7 +17,7 @@ struct Block { temperature(Rcpp::as(x["temperature"])), battery(Rcpp::as(x["battery"])), light(Rcpp::as(x["light"])), length(Rcpp::as(x["length"])), id(Rcpp::as(x["blockID"])), - data(x["data"]) {} + data(Rcpp::as(x["data"])) {} }; struct QC { @@ -44,7 +44,8 @@ SEXP qualityLog(const std::vector& rows) { // The R checksum-failure rows use ii + 1 (a double) for blockID_next; // otherwise both IDs come directly from integer packet headers. Rcpp::RObject nextIDs = hasCorrupt ? - Rcpp::RObject(Rcpp::as(next)) : Rcpp::RObject(next); + Rcpp::RObject(Rcpp::wrap(Rcpp::as(next))) : + Rcpp::RObject(Rcpp::wrap(next)); return Rcpp::DataFrame::create( Rcpp::_["checksum_pass"] = checksum, Rcpp::_["blockID_current"] = current, Rcpp::_["blockID_next"] = nextIDs, @@ -98,14 +99,15 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend Block current = previous; double nextTime; if (ii == numDBlocks) { - // Preserve the established final-packet extrapolation, including the - // nominal-frequency fallback for a chunk containing only that packet. + // Estimate the final packet's duration from the preceding packet, or + // use its nominal frequency when it is the only requested packet. if (firstIteration) { prevStart = previous.start; nextTime = prevStart + previous.length / previous.frequency; } else { nextTime = (previous.start - prevStart) / prevLength * previous.length + previous.start; } + prevStart = previous.start; prevLength = previous.length; } else { if (offset + 512 > buffer.size()) { diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index 6f4db2f..19ccb7d 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -48,6 +48,35 @@ test_that("readAxivity returns the final block when it is the only requested blo } }) +test_that("multi-block reads place the final packet on its own time axis", { + for (device in c("ax3", "ax6")) { + cwafile = system.file(paste0("testfiles/", device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(cwafile, "raw", n = file.size(cwafile)) + n = (length(bytes) - 1024) / 512 + packet = bytes[1024 + (n - 1) * 512 + seq_len(512)] + decoded = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE) + single = readAxivity(cwafile, start = n - 1, end = n, desiredtz = "UTC") + preceding = readAxivity(cwafile, start = n - 2, end = n - 1, desiredtz = "UTC") + finalStart = single$data$time[1] + # Multi-block reads estimate spacing from the preceding packet, whereas + # the single-block path must use the nominal frequency. + spacing = (finalStart - preceding$data$time[1]) / preceding$header$blockLength + rawTime = finalStart + (seq_len(decoded$length) - 1) * spacing + axes = colnames(decoded$data) + for (start in c(0, n - 2)) { + result = readAxivity(cwafile, start = start, end = n, desiredtz = "UTC") + tail = result$data[result$data$time >= finalStart,] + expect_gt(nrow(tail), decoded$length / 2) + expect_lte(max(tail$time), max(rawTime) + 1e-6) + for (axis in axes) { + expected = approx(rawTime, decoded$data[,axis], xout = tail$time)$y + expect_equal(tail[[axis]], expected, tolerance = 1e-6, scale = 1) + } + } + } +}) + test_that("readAxivity reads data from AX3 file correctly", { cwafile = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread")[1] AX3 = readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 1, end = 4) @@ -64,11 +93,11 @@ test_that("readAxivity reads data from AX3 file correctly", { # ask for more data then there is in the file AX3b = readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 0, end = 1000) expect_equal(AX3b$header$device, "Axivity") - expect_equal(nrow(AX3b$data), 17599) + expect_equal(nrow(AX3b$data), 17600) expect_equal(ncol(AX3b$data), 7) expect_equal(AX3b$data$time[5], 1551174906.040, tolerance = .001, scale = 1) expect_equal(AX3b$data$temp[3], 25.5859, tolerance = 0.0001, scale = 1) - expect_equal(floor(sum(abs(AX3b$data[,2:4]))), 25367) + expect_equal(floor(sum(abs(AX3b$data[,2:4]))), 25405) }) @@ -329,11 +358,11 @@ test_that("readAxivity reads data from AX6 file correctly", { # ask for more data then there is in the file AX6b = readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 0, end = 1000) expect_equal(AX6b$header$device, "Axivity") - expect_equal(nrow(AX6b$data), 11428) + expect_equal(nrow(AX6b$data), 11429) expect_equal(ncol(AX6b$data), 10) expect_equal(AX6b$data$time[5], 1577131446.740, tolerance = .001, scale = 1) expect_equal(AX6b$data$temp[3], 27.34375, tolerance = 0.0001) - expect_equal(floor(sum(abs(AX6b$data[,2:4]))), 960046) + expect_equal(floor(sum(abs(AX6b$data[,2:4]))), 960062) expect_true(is.null(AX6b$QClog)) }) diff --git a/tools/benchmark-readAxivity.R b/tools/benchmark-readAxivity.R index 271a179..992a6ca 100644 --- a/tools/benchmark-readAxivity.R +++ b/tools/benchmark-readAxivity.R @@ -52,9 +52,10 @@ benchmark_readAxivity = function(reference_file, blocks = c(1000L, 10000L), repe for (n in blocks) { filename = make_axivity_benchmark_file(device, n) tryCatch({ - # Warm up both paths, and compare the entire output before timing. + # Compare away from EOF: older readers place the final packet on an + # incorrect time axis. Full-file timings below include the fixed tail. outputs = lapply(readers, function(reader) { - reader(filename, start = 0, end = n, configtz = "UTC") + reader(filename, start = 0, end = n - 1, configtz = "UTC") }) stopifnot(all(vapply(outputs, identical, logical(1), outputs$native))) rm(outputs) From 3dd07204596afd4132d61d0c45367aa662190cf5 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Oct 2026 18:15:59 +0200 Subject: [PATCH 05/14] Fix Axivity timestamp rounding on macOS ARM MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Preserve R’s rounding when constructing native timestamps Resolve interpolation test failures without relaxing tolerances Update checkout actions to v5 and document the fix in NEWS --- .github/workflows/r_codecov.yml | 2 +- .github/workflows/r_multiOS.yml | 4 ++-- .github/workflows/r_ubuntuOnly.yml | 4 ++-- NEWS.md | 2 ++ src/processAxivityBlocks.cpp | 12 ++++++++++-- 5 files changed, 17 insertions(+), 7 deletions(-) diff --git a/.github/workflows/r_codecov.yml b/.github/workflows/r_codecov.yml index eb448e9..855de02 100644 --- a/.github/workflows/r_codecov.yml +++ b/.github/workflows/r_codecov.yml @@ -17,7 +17,7 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: r-lib/actions/setup-r@v2 with: diff --git a/.github/workflows/r_multiOS.yml b/.github/workflows/r_multiOS.yml index 55219db..73acdb7 100644 --- a/.github/workflows/r_multiOS.yml +++ b/.github/workflows/r_multiOS.yml @@ -30,7 +30,7 @@ jobs: R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: r-lib/actions/setup-pandoc@v2 @@ -48,4 +48,4 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: upload-snapshots: true - build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")' \ No newline at end of file + build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")' diff --git a/.github/workflows/r_ubuntuOnly.yml b/.github/workflows/r_ubuntuOnly.yml index e49d56c..7d8aa8e 100644 --- a/.github/workflows/r_ubuntuOnly.yml +++ b/.github/workflows/r_ubuntuOnly.yml @@ -26,7 +26,7 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: r-lib/actions/setup-tinytex@v2 - run: tlmgr --version - uses: r-lib/actions/setup-r@v2 @@ -45,4 +45,4 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: - upload-snapshots: true \ No newline at end of file + upload-snapshots: true diff --git a/NEWS.md b/NEWS.md index 28902b2..fdc0be3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,8 @@ - Keep the native Axivity processing loop compatible with C++11 and C++14, and release its full output buffers before assembling trimmed results to reduce memory retention. +- Preserve R's timestamp rounding in the native Axivity loop on platforms using fused multiply-add, fixing small interpolation differences on macOS ARM. + # Changes in version 1.0.11 (release date:18-09-2026) - Added GENEActiv light calibration needed for device model GENEActiv 1.2 and above. Contributed by ActivInsights Ltd, the manufacturer of GENEActiv devices. diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index d0caf68..fde7d33 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -105,7 +105,9 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend prevStart = previous.start; nextTime = prevStart + previous.length / previous.frequency; } else { - nextTime = (previous.start - prevStart) / prevLength * previous.length + previous.start; + // As in R, round the duration before adding the epoch timestamp. + volatile double duration = (previous.start - prevStart) / prevLength * previous.length; + nextTime = duration + previous.start; } prevStart = previous.start; prevLength = previous.length; @@ -158,7 +160,13 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend rawTime.resize(length + 2); const double spacing = (nextTime - prevStart) / length; rawTime[1] = prevStart; - for (int i = 1; i < length; ++i) rawTime[i + 1] = prevStart + i * spacing; + for (int i = 1; i < length; ++i) { + // R's seq() rounds the multiplication before adding the epoch timestamp. + // Prevent fused multiply-add (notably on ARM) from moving a time anchor + // by one ULP, which is amplified by interpolation between nearby samples. + volatile double elapsed = i * spacing; + rawTime[i + 1] = prevStart + elapsed; + } rawTime[length + 1] = nextTime; if (firstIteration) { for (int axis = 0; axis < axes; ++axis) anchor[axis] = previous.data(0, axis); From b6540997636861c0f2b84fda3003c4b8955e4046 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 20 Oct 2026 08:47:59 +0200 Subject: [PATCH 06/14] Fix AX6 channel selection during gap imputation Use accelerometer channels for orientation normalization and the 1 g fallback, and fill gyroscope channels with zero. Add regression tests for both cases with linear and nearest-neighbour interpolation, and update NEWS.md. --- NEWS.md | 2 + src/processAxivityBlocks.cpp | 10 +++-- tests/testthat/test_processAxivityBlocks.R | 50 ++++++++++++++++++++++ 3 files changed, 59 insertions(+), 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index fdc0be3..7f8f79d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # Changes in version 1.0.12 (release date:??-??-2026) +- Fix AX6 gap imputation in `readAxivity()` to use the accelerometer channels for the normalised 1 g orientation or fallback, and fill gyroscope channels with zero. + - Speed up `readAxivity()` with C++ packet decoding and a buffered native processing loop for timestamp caching, integrity checks, imputation and resampling. Retain the R interface and timezone conversion. - Fix `readAxivity()` failing when a requested chunk contains only the final data block. Use the block's recorded sampling frequency when no preceding block is available to estimate timing. #98 diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index fde7d33..08b3b8e 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -75,6 +75,8 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend Block previous(first, as(first["start"])); List parameters = first["parameters"]; const int axes = as(parameters["Naxes"]); + // Acceleration follows the gyroscope channels in six-axis packets. + const int accelOffset = axes == 6 ? 3 : 0; List state = first["struc"]; double timestamp = as(state[0]); int oldSeconds = as(state[1]); @@ -200,18 +202,20 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend // rare callback avoids changing repeated normalisation by one ULP. NumericVector orientation(3); for (int axis = 0; axis < 3; ++axis) { - orientation[axis] = anchor[axis]; + orientation[axis] = anchor[accelOffset + axis]; } const double norm = as(squareRoot(sum(power(orientation, 2.0)))); if (norm > 0.8 && norm < 1.2) { - for (int axis = 0; axis < 3; ++axis) imputed[axis] = anchor[axis] / norm; + for (int axis = 0; axis < 3; ++axis) imputed[axis] = anchor[accelOffset + axis] / norm; } } int source = 1; while (position <= last && (impute ? time[position] < rawTime[length] : time[position] <= rawTime[length])) { if (impute) { - for (int axis = 0; axis < axes; ++axis) accel(position, axis) = axis < 3 ? imputed[axis] : 0; + for (int axis = 0; axis < axes; ++axis) { + accel(position, axis) = axis < accelOffset ? 0 : imputed[axis - accelOffset]; + } } else { while (source < length && rawTime[source] < time[position]) ++source; const double weight = (time[position] - rawTime[source - 1]) / diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R index 578d0fb..bdfc0cb 100644 --- a/tests/testthat/test_processAxivityBlocks.R +++ b/tests/testthat/test_processAxivityBlocks.R @@ -76,6 +76,56 @@ test_that("repeated imputation retains the R normalisation arithmetic", { } }) +test_that("AX6 gaps impute acceleration orientation and zero angular velocity", { + fixture = system.file("testfiles/ax6_testfile.cwa", package = "GGIRread") + original = readBin(fixture, "raw", n = 1024 + 6 * 512) + first = GGIRread:::decodeAxivityBlock(original[1025:1536], TRUE, NULL, FALSE) + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + + # Exercise both normalisation of a plausible orientation and the 1 g fallback. + # Non-zero gyro values must never be used as the acceleration orientation. + for (orientation in list(c(0.6, 0.8, 0.2), c(2, 0, 0))) { + bytes = original + rawAccel = round(orientation / first$parameters$accelScale) + acceleration = rawAccel * first$parameters$accelScale + norm = sqrt(sum(acceleration^2)) + expected = if (norm > 0.8 && norm < 1.2) acceleration / norm else c(0, 0, 1) + for (block in 0:5) { + positions = 1024 + block * 512 + seq_len(512) + packet = bytes[positions] + count = readBin(packet[29:30], integer(), size = 2, signed = FALSE, + endian = "little") + packet[30 + seq_len(count * 12)] = writeBin( + as.integer(rep(c(1000, 2000, 3000, rawAccel), count)), raw(), + size = 2, endian = "little") + packet[511:512] = as.raw(0) + checksum = (-sum(readBin(packet, integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packet[511:512] = writeBin(as.integer(checksum), raw(), size = 2, + endian = "little") + bytes[positions] = packet + } + # Corrupt block 2, leaving recorded data on both sides of the gap. + bytes[1024 + 2 * 512 + 1] = as.raw(0) + writeBin(bytes, filename) + for (interpolation in 1:2) { + expect_warning(result <- readAxivity(filename, start = 0, end = 5, + configtz = "UTC", + interpolationType = interpolation), + "Skipping corrupt block #2") + gap = result$QClog[result$QClog$imputed, ] + expect_equal(nrow(gap), 1) + # Stay inside the imputed span, away from interpolation at its boundary. + rows = result$data$time >= gap$start & result$data$time < gap$end - 0.1 + expect_true(any(rows)) + expect_true(all(as.matrix(result$data[rows, c("gx", "gy", "gz")]) == 0)) + expect_equal(unname(as.matrix(result$data[rows, c("x", "y", "z")])), + matrix(rep(expected, each = sum(rows)), ncol = 3)) + } + } +}) + test_that("a corrupt tail without a usable endpoint fails explicitly", { filename = system.file("testfiles/ax3_testfile_corrupt_blocks_0_13_14_142_143_144.cwa", package = "GGIRread") From 927bc59f5bdcf013caf89b7b7307caae6d4b2609 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 20 Oct 2026 08:50:24 +0200 Subject: [PATCH 07/14] Fix loss of final valid Axivity packet before corrupt tail Treat the last valid packet as the effective end of the recording so it is processed without lookahead, including single-packet reads. Add AX3 and AX6 regression coverage, update existing tests and NEWS.md. --- NEWS.md | 2 + R/readAxivity.R | 4 ++ tests/testthat/test_processAxivityBlocks.R | 43 +++++++++++++++++----- tests/testthat/test_readAxivity.R | 28 +++++++------- 4 files changed, 53 insertions(+), 24 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7f8f79d..9ab183c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # Changes in version 1.0.12 (release date:??-??-2026) +- Fix `readAxivity()` dropping the last valid data block before a corrupt file tail, including reads that request only that valid block and the corrupt tail. + - Fix AX6 gap imputation in `readAxivity()` to use the accelerometer channels for the normalised 1 g orientation or fallback, and fill gyroscope channels with zero. - Speed up `readAxivity()` with C++ packet decoding and a buffered native processing loop for timestamp caching, integrity checks, imputation and resampling. Retain the R interface and timezone conversion. diff --git a/R/readAxivity.R b/R/readAxivity.R index 8e9db25..8d0f556 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -267,6 +267,10 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire } stop("Error reading file. The last ", maxAllowedCorruptBlocks+1, " blocks are corrupt.") } + # Treat the last valid packet as the physical end for native processing. + # Otherwise the loop expects a lookahead packet and never resamples it. + # The header retains the original file's block count. + numDBlocks = end # the end timestamp should fall right after the actual very last timestamp of the file endTimestamp = as.numeric(lastBlock$start) + blockLength * step # now pad it generously in case there are gaps in the last block diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R index bdfc0cb..af3fbb5 100644 --- a/tests/testthat/test_processAxivityBlocks.R +++ b/tests/testthat/test_processAxivityBlocks.R @@ -126,15 +126,40 @@ test_that("AX6 gaps impute acceleration orientation and zero angular velocity", } }) -test_that("a corrupt tail without a usable endpoint fails explicitly", { - filename = system.file("testfiles/ax3_testfile_corrupt_blocks_0_13_14_142_143_144.cwa", - package = "GGIRread") - # The former R loop left 'last' undefined after skipping every lookahead - # packet. Do not return a successful but empty data frame in this case. - expect_warning(expect_error(readAxivity(filename, start = 141, end = 145, - configtz = "UTC"), - "No samples could be resampled from the Axivity chunk"), - "Skipping corrupt") +test_that("corrupt tails retain the last valid packet as an effective final packet", { + filename = tempfile(fileext = ".cwa") + trimmed = tempfile(fileext = ".cwa") + on.exit(unlink(c(filename, trimmed))) + for (device in c("ax3", "ax6")) { + fixture = system.file("testfiles", paste0(device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(fixture, "raw", n = file.size(fixture)) + blocks = (length(bytes) - 1024) / 512 + validBlocks = blocks - 3 + writeBin(bytes[seq_len(1024 + validBlocks * 512)], trimmed) + for (block in validBlocks:(blocks - 1)) { + bytes[1024 + block * 512 + 1] = as.raw(0) + } + writeBin(bytes, filename) + for (interpolation in 1:2) { + # Cover whole-file reads, two valid packets and a single valid packet. + for (start in c(0, validBlocks - 2, validBlocks - 1)) { + expected = readAxivity(trimmed, start = start, end = validBlocks, + configtz = "UTC", interpolationType = interpolation) + expect_gt(nrow(expected$data), 0) + # Include requests whose endpoint search reaches the corrupt tail. + for (end in validBlocks:blocks) { + expect_warning(actual <- readAxivity(filename, start = start, end = end, + configtz = "UTC", + interpolationType = interpolation), + "Skipping corrupt end block") + expect_identical(actual$data, expected$data) + expect_identical(actual$QClog, expected$QClog) + expect_equal(actual$header$blocks, blocks) + } + } + } + } }) test_that("a device sequence reset refreshes the native timestamp cache", { diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index 19ccb7d..2b2461f 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -297,10 +297,14 @@ test_that("readAxivity handles failed checksums correctly", { "All requested blocks are corrupt"), "Skipping corrupt end block") - # If there are corrupt blocks at the very end of the file, then the last non-corrupt block is - # treated as the end block and no imputation is done. - # So in our case, since the end blocks 142-144 are corrupt, they are ignored and the file is - # treated as if it had only 142 blocks (0-141) + # Ignore the corrupt tail and include the final valid block, just as if + # the file ended after block 141. + trimmed = tempfile(fileext = ".cwa") + on.exit(unlink(trimmed)) + writeBin(readBin(cwafile, "raw", n = 1024 + 142 * 512), trimmed) + expect_warning(expectedTail <- readAxivity(trimmed, desiredtz = "Europe/Berlin", + start = 140, end = 142), + "Skipping corrupt block #0") expect_warning(AX3_140_141 <- readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 140, end = 141), "Skipping corrupt block #0") expect_warning(AX3_140_142 <- readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 140, end = 142), @@ -312,17 +316,11 @@ test_that("readAxivity handles failed checksums correctly", { expect_warning(AX3_140_145 <- readAxivity(filename = cwafile, desiredtz = "Europe/Berlin", start = 140, end = 145), "Skipping corrupt end block") - expect_equal(nrow(AX3_140_141$data), nrow(AX3_140_142$data)) - expect_equal(floor(sum(abs(AX3_140_141$data[,2:4]))), floor(sum(abs(AX3_140_142$data[,2:4])))) - - expect_equal(nrow(AX3_140_141$data), nrow(AX3_140_143$data)) - expect_equal(floor(sum(abs(AX3_140_141$data[,2:4]))), floor(sum(abs(AX3_140_143$data[,2:4])))) - - expect_equal(nrow(AX3_140_141$data), nrow(AX3_140_144$data)) - expect_equal(floor(sum(abs(AX3_140_141$data[,2:4]))), floor(sum(abs(AX3_140_144$data[,2:4])))) - - expect_equal(nrow(AX3_140_141$data), nrow(AX3_140_145$data)) - expect_equal(floor(sum(abs(AX3_140_141$data[,2:4]))), floor(sum(abs(AX3_140_145$data[,2:4])))) + expect_gt(nrow(expectedTail$data), nrow(AX3_140_141$data)) + for (actual in list(AX3_140_142, AX3_140_143, AX3_140_144, AX3_140_145)) { + expect_identical(actual$data, expectedTail$data) + expect_identical(actual$QClog, expectedTail$QClog) + } # If the file has more consecutive corrupt blocks than allowed by maxAllowedCorruptBlocks, # we should get an error. From 13c87a328c27094c646fdfc3473035ef2c07e647 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 20 Oct 2026 08:55:49 +0200 Subject: [PATCH 08/14] Preserve fractional nominal Axivity sampling rates Remove integer rounding from header and packet rate decoding so rates such as 12.5 Hz are retained in timestamp corrections and the uniform output grid. Add regression coverage and update NEWS.md. --- NEWS.md | 2 ++ R/readAxivity.R | 2 +- src/decodeAxivityBlock.cpp | 3 +- tests/testthat/test_decodeAxivityBlock.R | 22 ++++++++++++ tests/testthat/test_processAxivityBlocks.R | 16 +++++---- tests/testthat/test_readAxivity.R | 39 ++++++++++++++++++++++ 6 files changed, 75 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index 9ab183c..d8d60a1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # Changes in version 1.0.12 (release date:??-??-2026) +- Preserve fractional Axivity sampling rates (such as 12.5, 6.25 and 3.125 Hz) in header and packet decoding, timestamp corrections and the resampling grid instead of rounding them to integers. + - Fix `readAxivity()` dropping the last valid data block before a corrupt file tail, including reads that request only that valid block and the corrupt tail. - Fix AX6 gap imputation in `readAxivity()` to use the accelerometer channels for the normalised 1 g orientation or fallback, and fill gyroscope channels with zero. diff --git a/R/readAxivity.R b/R/readAxivity.R index 8d0f556..77c9c2f 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -127,7 +127,7 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire uniqueSerialCode = bitwOr(bitwShiftL(upperDeviceId, 16), lowerDeviceId) # sample rate and dynamic range accelerometer samplerate_dynrange = readBin(block[37], integer(), size = 1, signed = FALSE) #offset 36 - frequency_header = round( 3200 / bitwShiftL(1, 15 - bitwAnd(samplerate_dynrange, 15))) + frequency_header = 3200 / bitwShiftL(1, 15 - bitwAnd(samplerate_dynrange, 15)) accrange = bitwShiftR(16, (bitwShiftR(samplerate_dynrange, 6))) version = readBin(block[42], integer(), size = 1, signed = FALSE) #offset 41 diff --git a/src/decodeAxivityBlock.cpp b/src/decodeAxivityBlock.cpp index 2693f89..b0b476a 100644 --- a/src/decodeAxivityBlock.cpp +++ b/src/decodeAxivityBlock.cpp @@ -69,9 +69,8 @@ Rcpp::List decodeAxivityBlock(Rcpp::RawVector block, bool complete, axes = axesPacking >> 4; scaleCode = lightScale >> 13; accelScale = 1.0 / std::pow(2.0, 8 + scaleCode); - // nearbyint, like R's round(), rounds ties to even. frequency = samplerate == 0 ? offset26 : - std::nearbyint(3200.0 / (1U << (15 - (samplerate & 15)))); + 3200.0 / (1U << (15 - (samplerate & 15))); format = samplerate == 0 ? 0 : ((tsOffset & 0x8000U) != 0 ? 1 : 2); params = List::create(_["accelScaleCode"] = scaleCode, _["accelScale"] = accelScale, diff --git a/tests/testthat/test_decodeAxivityBlock.R b/tests/testthat/test_decodeAxivityBlock.R index 655c765..a99a278 100644 --- a/tests/testthat/test_decodeAxivityBlock.R +++ b/tests/testthat/test_decodeAxivityBlock.R @@ -109,6 +109,28 @@ test_that("native decoding preserves timestamp offsets and parameter reuse", { expect_equal(GGIRread:::decodeAxivityBlock(packet, FALSE, decoded$parameters, FALSE)$frequency, 200) }) +test_that("native decoding preserves fractional rates in timestamp corrections", { + for (code in 0:15) { + rate = 3200 / 2^(15 - code) + for (fractional in c(FALSE, TRUE)) { + packet = axivity_packet() + # Keep the rate byte nonzero even for code 0 (zero denotes legacy format). + packet[25] = as.raw(0x40 + code) + packet[5:6] = writeBin(if (fractional) 0xc000L else 0L, raw(), + size = 2, endian = "little") + packet[27:28] = writeBin(-5L, raw(), size = 2, endian = "little") + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, FALSE, NULL, FALSE) + expect_equal(decoded$frequency, rate) + expect_equal(decoded$parameters$frequency_data, rate) + expect_equal(decoded$fractional, if (fractional) 32768 else 0) + expect_equal(decoded$shift, (5 - if (fractional) floor(rate / 2) else 0) / rate) + expect_identical(GGIRread:::decodeAxivityBlock(packet, FALSE, + decoded$parameters, FALSE), decoded) + } + } +}) + test_that("native decoding applies all accelerometer and gyroscope scale codes", { for (device in c("ax3", "ax6")) { original = axivity_packet(device) diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R index af3fbb5..fa51662 100644 --- a/tests/testthat/test_processAxivityBlocks.R +++ b/tests/testthat/test_processAxivityBlocks.R @@ -1,21 +1,25 @@ library(GGIRread) test_that("native loop crosses input buffers and calendar days at a low sample rate", { - # Two days at 1 Hz exercise date changes without allocating a multi-day + # Two days at 0.78125 Hz exercise date changes without allocating a multi-day # 100 Hz output. Payloads are constant; timestamps and checksums are genuine. fixture = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") bytes = readBin(fixture, "raw", n = 1536) packet = bytes[1025:1536] - packet[25] = as.raw(3) # round(3200 / 2^(15 - 3)) = 1 Hz + packet[25] = as.raw(3) # 3200 / 2^(15 - 3) = 0.78125 Hz + rate = 0.78125 + count = 100L # 128 seconds per packet, keeping calendar timestamps integral. + packet[29:30] = writeBin(count, raw(), size = 2, endian = "little") packet[5:6] = as.raw(0) packet[27:28] = as.raw(0) packet[31:510] = rep(writeBin(256L, raw(), size = 4, endian = "little"), 120) blocks = 1442L origin = as.POSIXct("2020-01-01", tz = "UTC") - times = as.POSIXlt(origin + (seq_len(blocks) - 1) * 120, tz = "UTC") + times = as.POSIXlt(origin + (seq_len(blocks) - 1) * count / rate, tz = "UTC") coded = bitwShiftL(as.integer(times$year - 100), 26) + bitwShiftL(as.integer(times$mon + 1), 22) + bitwShiftL(as.integer(times$mday), 17) + - bitwShiftL(as.integer(times$hour), 12) + bitwShiftL(as.integer(times$min), 6) + bitwShiftL(as.integer(times$hour), 12) + bitwShiftL(as.integer(times$min), 6) + + as.integer(times$sec) packets = matrix(rep(packet, blocks), nrow = 512) packets[11:14,] = writeBin(seq_len(blocks) - 1L, raw(), size = 4, endian = "little") packets[15:18,] = writeBin(as.integer(coded), raw(), size = 4, endian = "little") @@ -30,9 +34,9 @@ test_that("native loop crosses input buffers and calendar days at a low sample r on.exit(unlink(filename)) writeBin(c(bytes[1:1024], as.vector(packets)), filename) result = readAxivity(filename, start = 0, end = blocks - 1, configtz = "UTC") - samples = (blocks - 1) * 120 + samples = (blocks - 1) * count expect_equal(nrow(result$data), samples) - expect_identical(result$data$time, as.numeric(origin) + seq_len(samples) - 1) + expect_identical(result$data$time, as.numeric(origin) + (seq_len(samples) - 1) / rate) scale = 1 / 2^(8 + bitwShiftR(as.integer(packet[20]), 5)) expect_true(all(result$data$x == 256 * scale)) expect_true(all(result$data$y == 0)) diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index 2b2461f..d8a902c 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -1,6 +1,45 @@ library(GGIRread) context("reading Axivity .cwa data") +test_that("fractional sampling rates are preserved in headers and output grids", { + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + for (device in c("ax3", "ax6")) { + fixture = system.file("testfiles", paste0(device, "_testfile.cwa"), + package = "GGIRread") + original = readBin(fixture, "raw", n = 1536) + for (code in 0:7) { + rate = 3200 / 2^(15 - code) + bytes = original + bytes[37] = as.raw(0x40 + code) + packet = bytes[1025:1536] + packet[25] = as.raw(0x40 + code) + packet[5:6] = writeBin(0xc000L, raw(), size = 2, endian = "little") + packet[15:18] = writeBin(as.integer(20 * 2^26 + 2^22 + 2^17), raw(), + size = 4, endian = "little") # 2020-01-01 UTC + packet[27:28] = writeBin(-5L, raw(), size = 2, endian = "little") + packet[511:512] = as.raw(0) + checksum = (-sum(readBin(packet, integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packet[511:512] = writeBin(as.integer(checksum), raw(), size = 2, + endian = "little") + bytes[1025:1536] = packet + writeBin(bytes, filename) + result = readAxivity(filename, start = 0, end = 1, configtz = "UTC") + expectedStart = as.numeric(as.POSIXct("2020-01-01", tz = "UTC")) + + 0.5 + (5 - floor(rate / 2)) / rate + expect_equal(result$header$frequency, rate) + expect_equal(as.numeric(result$header$start), expectedStart, tolerance = 1e-6, + scale = 1) + expect_equal(nrow(result$data), result$header$blockLength) + expect_equal(result$data$time, + expectedStart + (seq_len(result$header$blockLength) - 1) / rate, + tolerance = 1e-6, scale = 1) + expect_null(result$QClog) + } + } +}) + test_that("readAxivity returns the final block when it is the only requested block", { for (device in c("ax3", "ax6")) { cwafile = system.file(paste0("testfiles/", device, "_testfile.cwa"), From 8b66687c7bc8ee1cb41f7bc66b2bfa3b04caa713 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 20 Oct 2026 09:07:10 +0200 Subject: [PATCH 09/14] Fix Axivity truncated-packet handling and final-sample timing Count only complete packets, warn about ignored trailing bytes, and report an error when no complete data packets remain. Use nominal spacing for single final blocks to prevent floating-point rounding from dropping the last sample at fractional sampling rates. Add regression tests and update NEWS.md. Full suite: 2,780 tests pass. --- NEWS.md | 4 ++ R/readAxivity.R | 16 ++++++- src/processAxivityBlocks.cpp | 7 ++- tests/testthat/test_processAxivityBlocks.R | 20 +++++++++ tests/testthat/test_readAxivity.R | 52 ++++++++++++++++++++++ 5 files changed, 96 insertions(+), 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index d8d60a1..4af52ea 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # Changes in version 1.0.12 (release date:??-??-2026) +- Fix a floating-point timing error that could drop the final sample when reading a single final Axivity block at a fractional sampling rate. + +- Fix `readAxivity()` handling of truncated final packets: count only complete data blocks and warn about ignored trailing bytes. Report an explicit error when no complete data blocks remain. + - Preserve fractional Axivity sampling rates (such as 12.5, 6.25 and 3.125 Hz) in header and packet decoding, timestamp corrections and the resampling grid instead of rounding them to integers. - Fix `readAxivity()` dropping the last valid data block before a corrupt file tail, including reads that request only that valid block and the corrupt tail. diff --git a/R/readAxivity.R b/R/readAxivity.R index 77c9c2f..96a7620 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -181,14 +181,26 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire if (nargin < 1) { stop("At least file must be specified") } - # Get file size in data blocks - numDBlocks = round(file.size(filename) / blockBytes) - 2 + # Only complete packets can be decoded; a partial final packet is excluded. + fileBytes = file.size(filename) + numDBlocks = floor(fileBytes / blockBytes) - 2 # Open file fid = file(filename,"rb") on.exit({ close(fid) }) + if (fileBytes >= headerBytes) { + trailingBytes = (fileBytes - headerBytes) %% blockBytes + if (trailingBytes > 0) { + warning("Ignoring ", trailingBytes, + " trailing bytes from an incomplete Axivity data block.") + } + if (numDBlocks == 0) { + stop("Error reading file. No complete Axivity data blocks.") + } + } + QClog = NULL # Initialise log of data quality issues ############################################################################# diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index 08b3b8e..4ed5453 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -100,6 +100,7 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend if (static_cast(ii - start) % 128 == 0) checkUserInterrupt(); Block current = previous; double nextTime; + const bool nominalFinal = ii == numDBlocks && firstIteration; if (ii == numDBlocks) { // Estimate the final packet's duration from the preceding packet, or // use its nominal frequency when it is the only requested packet. @@ -160,7 +161,11 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend if (!std::isfinite(nextTime)) stop("'to' must be a finite number"); const int length = prevLength; rawTime.resize(length + 2); - const double spacing = (nextTime - prevStart) / length; + // A single final packet has a known nominal spacing. Subtracting rounded + // epoch timestamps to recover it can move its last sample one ULP before + // the output grid and incorrectly exclude that sample. + const double spacing = nominalFinal ? 1.0 / previous.frequency : + (nextTime - prevStart) / length; rawTime[1] = prevStart; for (int i = 1; i < length; ++i) { // R's seq() rounds the multiplication before adding the epoch timestamp. diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R index fa51662..36d9ca0 100644 --- a/tests/testthat/test_processAxivityBlocks.R +++ b/tests/testthat/test_processAxivityBlocks.R @@ -1,5 +1,25 @@ library(GGIRread) +test_that("single final packets retain their last sample despite epoch rounding", { + for (length in c(40L, 120L)) { + rate = 1.5625 + start = 1577836800.5 + 5 / rate + data = cbind(x = seq_len(length), y = 0, z = 0) + first = list(start = start, frequency = rate, temperature = 20, battery = 0, + light = 0, length = length, blockID = 0L, data = data, + parameters = list(Naxes = 3L), struc = list(0, 0L, 0)) + time = start + (seq_len(length) - 1) * (1 / rate) + for (interpolation in 1:2) { + result = GGIRread:::processAxivityBlocks( + function() stop("No lookahead is needed"), function(x) 0, NULL, + first, time, 0, 1, 1, 0.1, 20L, FALSE, interpolation) + expect_equal(result$last, length) + expect_equal(unname(result$accel), unname(data)) + expect_null(result$QClog) + } + } +}) + test_that("native loop crosses input buffers and calendar days at a low sample rate", { # Two days at 0.78125 Hz exercise date changes without allocating a multi-day # 100 Hz output. Payloads are constant; timestamps and checksums are genuine. diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index d8a902c..a9cbfd0 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -1,6 +1,58 @@ library(GGIRread) context("reading Axivity .cwa data") +test_that("incomplete final packets are reported and excluded consistently", { + filename = tempfile(fileext = ".cwa") + complete = tempfile(fileext = ".cwa") + on.exit(unlink(c(filename, complete))) + for (device in c("ax3", "ax6")) { + fixture = system.file("testfiles", paste0(device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(fixture, "raw", n = file.size(fixture)) + # The reference ends immediately before the packet that will be truncated. + writeBin(head(bytes, -512), complete) + for (interpolation in 1:2) { + expected = readAxivity(complete, start = 0, end = Inf, configtz = "UTC", + interpolationType = interpolation) + # Cover both sides of the former rounding boundary, including the tie. + for (remaining in c(1, 255, 256, 257, 511)) { + writeBin(head(bytes, -(512 - remaining)), filename) + message = paste0("Ignoring ", remaining, + " trailing bytes from an incomplete Axivity data block") + expect_warning(actual <- readAxivity(filename, start = 0, end = Inf, + configtz = "UTC", + interpolationType = interpolation), + message) + expect_identical(actual, expected) + # Also cover cached-header reads requesting only the final complete packet. + start = expected$header$blocks - 1 + reference = readAxivity(complete, start = start, end = Inf, + configtz = "UTC", header = expected$header, + interpolationType = interpolation) + expect_warning(actual <- readAxivity(filename, start = start, end = Inf, + configtz = "UTC", header = expected$header, + interpolationType = interpolation), + message) + expect_identical(actual, reference) + } + } + } +}) + +test_that("a truncated first data packet gives an explicit error", { + fixture = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") + bytes = readBin(fixture, "raw", n = 1536) + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + for (remaining in c(1, 255, 256, 257, 511)) { + writeBin(bytes[seq_len(1024 + remaining)], filename) + expect_warning( + expect_error(readAxivity(filename, start = 0, end = Inf, configtz = "UTC"), + "No complete Axivity data blocks"), + paste0("Ignoring ", remaining, " trailing bytes")) + } +}) + test_that("fractional sampling rates are preserved in headers and output grids", { filename = tempfile(fileext = ".cwa") on.exit(unlink(filename)) From 820e113d4367da251a22f02918f991825c837874 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 20 Oct 2026 09:13:35 +0200 Subject: [PATCH 10/14] Fix Axivity interpolation anchors at packet boundaries Use the preceding raw sample and timestamp instead of the last resampled output for linear and nearest-neighbour interpolation. Preserve the existing pre-gap imputation policy. Add regression coverage and update NEWS.md. Full suite: 2,802 tests pass. --- NEWS.md | 2 + src/processAxivityBlocks.cpp | 23 +++++++-- tests/testthat/test_processAxivityBlocks.R | 60 ++++++++++++++++++++++ tests/testthat/test_readAxivity.R | 3 +- 4 files changed, 83 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index 4af52ea..603ddcf 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # Changes in version 1.0.12 (release date:??-??-2026) +- Fix interpolation between valid Axivity packets to use the preceding raw sample and timestamp instead of the last resampled output, for both linear and nearest-neighbour interpolation. + - Fix a floating-point timing error that could drop the final sample when reading a single final Axivity block at a fractional sampling rate. - Fix `readAxivity()` handling of truncated final packets: count only complete data blocks and warn about ignored trailing bytes. Report an explicit error when no complete data blocks remain. diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index 4ed5453..33709e7 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -89,8 +89,10 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend Function sum = Environment::base_env()["sum"]; Function power = Environment::base_env()["^"]; Function squareRoot = Environment::base_env()["sqrt"]; - std::vector anchor(axes, 0.0), rawTime; - double anchorTime = 0, prevStart = 0; + // Keep the output anchor for quality checks and imputation, independently + // of the raw-sample anchor used for interpolation between valid packets. + std::vector anchor(axes, 0.0), interpolationAnchor(axes, 0.0), rawTime; + double anchorTime = 0, interpolationTime = 0, prevStart = 0; int prevLength = 0, position = 0, consecutiveCorrupt = 0; bool firstIteration = true; RawVector buffer(0), packet(512); @@ -178,10 +180,12 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend if (firstIteration) { for (int axis = 0; axis < axes; ++axis) anchor[axis] = previous.data(0, axis); anchorTime = prevStart - 0.00001; + interpolationAnchor = anchor; + interpolationTime = anchorTime; firstIteration = false; } - rawTime[0] = anchorTime; double observed = (length + 1.0) / (rawTime[length] - anchorTime); + rawTime[0] = interpolationTime; if (ii == numDBlocks && start == numDBlocks - 1) observed = previous.frequency; const double bias = std::abs(observed - previous.frequency) / previous.frequency; const bool impute = (ii < numDBlocks && current.id != 0 && @@ -226,7 +230,7 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend const double weight = (time[position] - rawTime[source - 1]) / (rawTime[source] - rawTime[source - 1]); for (int axis = 0; axis < axes; ++axis) { - const double left = source == 1 ? anchor[axis] : previous.data(source - 2, axis); + const double left = source == 1 ? interpolationAnchor[axis] : previous.data(source - 2, axis); const double right = previous.data(source - 1, 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; @@ -240,6 +244,17 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend if (position == 0) stop("No samples could be resampled from the Axivity block."); anchorTime = time[position - 1]; for (int axis = 0; axis < axes; ++axis) anchor[axis] = accel(position - 1, axis); + if (impute) { + // Continue from the imputed output after a rejected interval; do not + // reintroduce its raw samples at the boundary after a gap. + interpolationTime = anchorTime; + interpolationAnchor = anchor; + } else { + interpolationTime = rawTime[length]; + for (int axis = 0; axis < axes; ++axis) { + interpolationAnchor[axis] = previous.data(length - 1, axis); + } + } previous = current; if (progress.isNotNull()) Function(progress.get())(position + 1); } diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R index 36d9ca0..6f654f7 100644 --- a/tests/testthat/test_processAxivityBlocks.R +++ b/tests/testthat/test_processAxivityBlocks.R @@ -1,5 +1,65 @@ library(GGIRread) +test_that("packet boundaries interpolate through raw samples on every axis", { + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + origin = as.numeric(as.POSIXct("2020-01-01", tz = "UTC")) + for (device in c("ax3", "ax6")) { + fixture = system.file("testfiles", paste0(device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(fixture, "raw", n = 1536) + axes = if (device == "ax3") 3L else 6L + count = 32L + packets = list() + values = list() + for (block in 0:3) { + packet = bytes[1025:1536] + elapsed = block * 0.375 + fraction = elapsed %% 1 + packet[5:6] = writeBin(as.integer(0x8000 + fraction * 32768), raw(), + size = 2, endian = "little") + packet[11:14] = writeBin(block, raw(), size = 4, endian = "little") + packet[15:18] = writeBin(as.integer(20 * 2^26 + 2^22 + 2^17 + floor(elapsed)), + raw(), size = 4, endian = "little") + packet[25] = as.raw(10) # Nominal output rate: 100 Hz. + packet[26] = as.raw(axes * 16 + 2) # Unpacked signed 16-bit samples. + packet[27:28] = writeBin(as.integer(-floor(fraction * 100)), raw(), + size = 2, endian = "little") + packet[29:30] = writeBin(count, raw(), size = 2, endian = "little") + samples = matrix((seq_len(count * axes) * 137 + block * 997) %% 4000 - 2000, + ncol = axes) + packet[30 + seq_len(count * axes * 2)] = writeBin( + as.integer(t(samples)), raw(), size = 2, endian = "little") + packet[511:512] = as.raw(0) + checksum = (-sum(readBin(packet, integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packet[511:512] = writeBin(as.integer(checksum), raw(), size = 2, + endian = "little") + packets[[block + 1]] = packet + values[[block + 1]] = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE)$data + } + 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]) + for (interpolation in 1:2) { + result = readAxivity(filename, start = 0, end = 3, configtz = "UTC", + interpolationType = interpolation, frequency_tol = 0.5) + expect_false(any(result$QClog$imputed)) + for (axis in colnames(rawValues)) { + if (interpolation == 1) { + expected = approx(rawTime, rawValues[,axis], xout = result$data$time)$y + } else { + nearest = vapply(result$data$time, function(t) which.min(abs(rawTime - t)), 1L) + expected = rawValues[nearest, axis] + } + expect_equal(result$data[[axis]], expected, tolerance = 1e-9) + } + } + } +}) + test_that("single final packets retain their last sample despite epoch rounding", { for (length in c(40L, 120L)) { rate = 1.5625 diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index a9cbfd0..a8ca49f 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -451,7 +451,8 @@ test_that("readAxivity reads data from AX6 file correctly", { expect_equal(ncol(AX6b$data), 10) expect_equal(AX6b$data$time[5], 1577131446.740, tolerance = .001, scale = 1) expect_equal(AX6b$data$temp[3], 27.34375, tolerance = 0.0001) - expect_equal(floor(sum(abs(AX6b$data[,2:4]))), 960062) + # Includes interpolation through raw samples at packet boundaries. + expect_equal(floor(sum(abs(AX6b$data[,2:4]))), 960636) expect_true(is.null(AX6b$QClog)) }) From df708c80dc5e2c347cdf599dc6f9d925297ee241 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Fri, 23 Oct 2026 12:03:54 +0200 Subject: [PATCH 11/14] Fix Axivity fractional timestamp shim and repeated chunk warnings Use the firmware's truncated integer frequency to undo the timestamp shim while preserving fractional rates for sample timing. Warn about trailing bytes only when reading the header. Add regression tests and update NEWS. --- NEWS.md | 4 +++- R/readAxivity.R | 2 +- src/decodeAxivityBlock.cpp | 4 +++- tests/testthat/test_decodeAxivityBlock.R | 18 ++++++++++++++++++ tests/testthat/test_readAxivity.R | 6 +++--- 5 files changed, 28 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index 603ddcf..dba74df 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,10 +4,12 @@ - Fix a floating-point timing error that could drop the final sample when reading a single final Axivity block at a fractional sampling rate. -- Fix `readAxivity()` handling of truncated final packets: count only complete data blocks and warn about ignored trailing bytes. Report an explicit error when no complete data blocks remain. +- Fix `readAxivity()` handling of truncated final packets: count only complete data blocks and warn about ignored trailing bytes only when reading the header, avoiding repeated warnings for subsequent chunks. Report an explicit error when no complete data blocks remain. - Preserve fractional Axivity sampling rates (such as 12.5, 6.25 and 3.125 Hz) in header and packet decoding, timestamp corrections and the resampling grid instead of rounding them to integers. +- Undo the Axivity firmware's fractional timestamp shim using its truncated integer sampling frequency while retaining the exact nominal rate for sample timing. + - Fix `readAxivity()` dropping the last valid data block before a corrupt file tail, including reads that request only that valid block and the corrupt tail. - Fix AX6 gap imputation in `readAxivity()` to use the accelerometer channels for the normalised 1 g orientation or fallback, and fill gyroscope channels with zero. diff --git a/R/readAxivity.R b/R/readAxivity.R index 96a7620..8bec4aa 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -192,7 +192,7 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire if (fileBytes >= headerBytes) { trailingBytes = (fileBytes - headerBytes) %% blockBytes - if (trailingBytes > 0) { + if (trailingBytes > 0 && is.null(header)) { warning("Ignoring ", trailingBytes, " trailing bytes from an incomplete Axivity data block.") } diff --git a/src/decodeAxivityBlock.cpp b/src/decodeAxivityBlock.cpp index b0b476a..bb677e8 100644 --- a/src/decodeAxivityBlock.cpp +++ b/src/decodeAxivityBlock.cpp @@ -94,8 +94,10 @@ Rcpp::List decodeAxivityBlock(Rcpp::RawVector block, bool complete, if (format == 1) { // Undo the firmware's backwards-compatible timestamp shim: the lower // 15 bits hold a fraction of a second, and offset26 includes its samples. + // The firmware truncates the frequency to an integer here (e.g. 12 for + // 12.5 Hz), so undo the shim with the same truncated value. fractional = (tsOffset & 0x7fffU) << 1; - shift += std::floor(fractional * frequency / 65536.0); + shift += std::floor(fractional * std::floor(frequency) / 65536.0); } } // Temperature uses the low ten ADC bits. Battery stores only the middle diff --git a/tests/testthat/test_decodeAxivityBlock.R b/tests/testthat/test_decodeAxivityBlock.R index a99a278..01411da 100644 --- a/tests/testthat/test_decodeAxivityBlock.R +++ b/tests/testthat/test_decodeAxivityBlock.R @@ -158,6 +158,24 @@ test_that("native decoding applies all accelerometer and gyroscope scale codes", } }) +test_that("the fractional timestamp shim uses the firmware's integer frequency", { + for (code in 5:7) { + packet = axivity_packet() + rate = 3200 / 2^(15 - code) # 3.125, 6.25 and 12.5 Hz. + packet[25] = as.raw(0x40 + code) + # Fraction 0.96875 s: at 12.5 Hz the firmware shifts 11 samples, + # whereas using the exact rate would incorrectly shift 12 samples. + packet[5:6] = writeBin(0xfc00L, raw(), size = 2, endian = "little") + packet[27:28] = writeBin(-5L, raw(), size = 2, endian = "little") + packet = axivity_checksum(packet) + decoded = GGIRread:::decodeAxivityBlock(packet, FALSE, NULL, FALSE) + expect_equal(decoded$frequency, rate) + expect_equal(decoded$shift, -(-5 + floor(0.96875 * floor(rate))) / rate) + expect_identical(GGIRread:::decodeAxivityBlock(packet, FALSE, + decoded$parameters, FALSE), decoded) + } +}) + test_that("native decoder validates packets before accessing sample payloads", { packet = axivity_packet() expect_error(GGIRread:::decodeAxivityBlock(packet[-1], TRUE, NULL, FALSE), "512 bytes") diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index a8ca49f..14e62a7 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -29,10 +29,10 @@ test_that("incomplete final packets are reported and excluded consistently", { reference = readAxivity(complete, start = start, end = Inf, configtz = "UTC", header = expected$header, interpolationType = interpolation) - expect_warning(actual <- readAxivity(filename, start = start, end = Inf, + # Warn only when reading the header, not on subsequent cached-header chunks. + expect_no_warning(actual <- readAxivity(filename, start = start, end = Inf, configtz = "UTC", header = expected$header, - interpolationType = interpolation), - message) + interpolationType = interpolation)) expect_identical(actual, reference) } } From 2dd9ce920de0221d0ed8533bce1379fce0fbd135 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Fri, 23 Oct 2026 12:12:11 +0200 Subject: [PATCH 12/14] Fix Axivity buffer leak when warning handlers exit Emit corrupt-block warnings through Rcpp's protected call to base R so output buffers are released when warning handlers exit or warn = 2. Require Rcpp >= 1.0.10, add memory regression tests and update NEWS. --- DESCRIPTION | 4 +-- NEWS.md | 2 ++ src/processAxivityBlocks.cpp | 6 +++- tests/testthat/test_processAxivityBlocks.R | 41 ++++++++++++++++++++++ 4 files changed, 50 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fa833b0..5b5cc46 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,8 +24,8 @@ URL: https://github.com/wadpac/GGIRread/ BugReports: https://github.com/wadpac/GGIRread/issues License: Apache License (== 2.0) Suggests: testthat -Imports: matlab, bitops, Rcpp (>= 0.12.10), data.table, readxl, jsonlite, digest +Imports: matlab, bitops, Rcpp (>= 1.0.10), data.table, readxl, jsonlite, digest Depends: stats, utils, R (>= 3.5.0) NeedsCompilation: yes -LinkingTo: Rcpp +LinkingTo: Rcpp (>= 1.0.10) ByteCompile: yes diff --git a/NEWS.md b/NEWS.md index dba74df..35c5ae1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # Changes in version 1.0.12 (release date:??-??-2026) +- Release native Axivity output buffers when a corrupt-block warning exits through a warning handler or `options(warn = 2)`. Require Rcpp >= 1.0.10 for unwind protection. + - Fix interpolation between valid Axivity packets to use the preceding raw sample and timestamp instead of the last resampled output, for both linear and nearest-neighbour interpolation. - Fix a floating-point timing error that could drop the final sample when reading a single final Axivity block at a fractional sampling rate. diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index 33709e7..c48668a 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -130,7 +130,11 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend if (!as(decoded["checksum_pass"])) { qc.push_back(QC{false, false, static_cast(ii), static_cast(ii + 1), 0, 0, 0, 0, 0}); - warning("Skipping corrupt block #%.0f", ii); + // Rcpp::warning() can longjmp past C++ destructors when a handler + // exits or warn = 2. Function calls use Rcpp's unwind protection so + // the preserved output buffers are released on those paths too. + Function(Environment::base_env()["warning"])( + tfm::format("Skipping corrupt block #%.0f", ii), _["call."] = false); if (++consecutiveCorrupt > maxAllowedCorruptBlocks) { stop("Error reading file. %d consecutive blocks are corrupt.", maxAllowedCorruptBlocks + 1); } diff --git a/tests/testthat/test_processAxivityBlocks.R b/tests/testthat/test_processAxivityBlocks.R index 6f654f7..135a68d 100644 --- a/tests/testthat/test_processAxivityBlocks.R +++ b/tests/testthat/test_processAxivityBlocks.R @@ -1,5 +1,46 @@ library(GGIRread) +test_that("exiting warning handlers release native output buffers", { + fixture = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") + packet = readBin(fixture, "raw", n = 1536)[1025:1536] + first = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE) + first$start = 0 + first$struc = list(0, 0L, 0) + packet[1] = as.raw(0) # A failed checksum triggers the native warning path. + samples = 100000L + time = (seq_len(samples) - 1) / first$frequency + read = function() GGIRread:::processAxivityBlocks( + function() packet, function(x) 0, NULL, first, time, + 0, 2, 2, 0.1, 20L, FALSE, 1L) + old = options(warn = 1) + on.exit(options(old)) + + condition = tryCatch(read(), warning = identity) + expect_s3_class(condition, "warning") + expect_identical(conditionMessage(condition), "Skipping corrupt block #1") + expect_null(conditionCall(condition)) + options(warn = 2) + expect_error(read(), "(converted from warning) Skipping corrupt block #1", fixed = TRUE) + + for (warn in c(1, 2)) { + options(warn = warn) + # Warm up condition handling before measuring live R vector cells. RSS + # need not fall after gc(), because the allocator can retain freed pages. + interrupt = if (warn == 1) { + function() tryCatch(read(), warning = function(w) NULL) + } else { + function() tryCatch(read(), error = function(e) NULL) + } + expect_null(interrupt()) + before = gc()["Vcells", "used"] + for (i in 1:3) expect_null(interrupt()) + retained = gc()["Vcells", "used"] - before + # Allow a whole output allocation of headroom; three leaked allocations + # retain at least 18 vector cells per sample (six numeric columns each). + expect_lt(retained, samples * 6) + } +}) + test_that("packet boundaries interpolate through raw samples on every axis", { filename = tempfile(fileext = ".cwa") on.exit(unlink(filename)) From d2e62b2545b456d493146b6a49499d929fc4a2c5 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Fri, 23 Oct 2026 12:27:52 +0200 Subject: [PATCH 13/14] Fix Axivity edge cases, documentation and Clang compatibility Preserve sensor columns in single-sample reads, enforce corrupt-block limits during endpoint searches and reject unsupported packing formats. Use a named warning callback to fix macOS compilation while retaining unwind protection. Update documentation, NEWS and regression tests. --- NEWS.md | 6 +- R/readAxivity.R | 17 +++-- man/readAxivity.Rd | 85 ++++++++++++++++-------- src/decodeAxivityBlock.cpp | 4 +- src/processAxivityBlocks.cpp | 4 +- tests/testthat/test_decodeAxivityBlock.R | 5 ++ tests/testthat/test_readAxivity.R | 50 ++++++++++++++ 7 files changed, 129 insertions(+), 42 deletions(-) diff --git a/NEWS.md b/NEWS.md index 35c5ae1..de426b4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,10 @@ # Changes in version 1.0.12 (release date:??-??-2026) -- Release native Axivity output buffers when a corrupt-block warning exits through a warning handler or `options(warn = 2)`. Require Rcpp >= 1.0.10 for unwind protection. +- Preserve sensor columns for Axivity reads returning a single sample, and report an explicit error when a chunk's endpoint search exceeds the consecutive corrupt-block limit. + +- Reject unsupported Axivity sample packing formats rather than interpreting them as 16-bit samples. Update the reader documentation to describe chunk indices, output units, header fields and quality logs. + +- Release native Axivity output buffers when a corrupt-block warning exits through a warning handler or `options(warn = 2)`. Require Rcpp >= 1.0.10 for unwind protection, and use a named warning callback for compatibility with Clang on macOS. - Fix interpolation between valid Axivity packets to use the preceding raw sample and timestamp instead of the last resampled output, for both linear and nearest-neighbour interpolation. diff --git a/R/readAxivity.R b/R/readAxivity.R index 8bec4aa..92cff46 100755 --- a/R/readAxivity.R +++ b/R/readAxivity.R @@ -96,7 +96,7 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire # as.POSIXct(start, origin = "1970-01-01", tz=desiredtz) # device is "Axivity" # firmwareVersion is version of firmware - # blocks is number of datablocks with 80 or 120 observations in each + # blocks is the number of complete data blocks (typically 40, 80 or 120 samples) # Unfortunately frequency of measurement is varied in this device. # @@ -246,8 +246,9 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire warning("Skipping corrupt end block #", ii) end = end + 1 } - if (is_corrupt && ii == end+maxAllowedCorruptBlocks) { - stop("Error reading file. The last ", maxAllowedCorruptBlocks+1, " blocks are corrupt.") + if (is_corrupt && end < numDBlocks) { + stop("Error reading file. ", maxAllowedCorruptBlocks+1, + " consecutive blocks are corrupt.") } endTimestamp = as.numeric(endBlock$start) } @@ -255,11 +256,11 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire if (end == numDBlocks) { # end == numDBlocks, meaning we'll be reading all the remaining blocks. # There is no block #numDBlocks, so we can't get the ending timestamp from the start of that block. - # Instead read the very last block of the file (if the last block is corrupt, fing the last non-corrupt one), + # Instead read the very last block of the file (if corrupt, find the last non-corrupt one), # then project what the ending timestamp should be. # Skip any corrupt blocks, up to maxAllowedCorruptBlocks in number. - # SInce we are at the last block, we have to go backwards. + # Since we are at the last block, we have to go backwards. is_corrupt = TRUE for (ii in (end-1) : max(start, end-maxAllowedCorruptBlocks-1)) { seek(fid, headerBytes + blockBytes * ii, origin = 'start') @@ -292,7 +293,6 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire # Read the start block. # Reinitiate file and skip header as well as the initial start-1 blocks seek(fid, headerBytes + blockBytes * start, origin = 'start') - pos = 1 # position of the first element to complete in data # Skip any corrupt blocks, up to maxAllowedCorruptBlocks in number. is_corrupt = TRUE @@ -326,7 +326,6 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire if (is.null(prevRaw)) { return(invisible(list(header = header, data = NULL))) } - struc = prevRaw$struc startTimestamp = as.numeric(prevRaw$start) # allocate memory for results @@ -360,8 +359,8 @@ readAxivity = function(filename, start = 0, end = 0, progressBar = FALSE, desire # If the user asked for more data than the length of the recording, # there will be 0s at the end of the result lists; get rid of them. if (last < nrow(accelRes)) { - cut = c(last+1:nrow(accelRes)) - accelRes = accelRes[-cut,] + cut = seq.int(last + 1L, nrow(accelRes)) + accelRes = accelRes[-cut, , drop = FALSE] battery = battery[-cut] light = light[-cut] temp = temp[-cut] diff --git a/man/readAxivity.Rd b/man/readAxivity.Rd index 1d5c500..9940677 100644 --- a/man/readAxivity.Rd +++ b/man/readAxivity.Rd @@ -16,12 +16,14 @@ filename (required) is name of cwa file to read } \item{start}{ - non-negative integer which is a cwa file block number. + Zero-based, non-negative integer index of the first data block to read. + The metadata header is not counted as a data block. } \item{end}{ - non-negative integer which is a cwa file block number. - End must be not be less than start. - If end is less or equal to start, then there is no data read. + Exclusive data block index. Use \code{Inf} to read through the last complete + data block. Values beyond the file's block count are limited to that count. + If \code{end <= start}, only the header is returned, with \code{data = NULL}. + The default \code{start = end = 0} therefore reads only the header. } \item{progressBar}{ Is trigger to switch on/off the text progress bar. If progressBar @@ -43,7 +45,9 @@ Boolean to indicate whether battery voltage should be loaded } \item{header}{ - Header to be reused if it was extracted earlier + Header returned by an earlier call for the same file, to avoid reading it + again when processing chunks. Use the same timezone settings as in that call. + Warnings about incomplete trailing bytes are issued only when \code{header = NULL}. } \item{frequency_tol}{ Numeric value representing the tolerance for frequency bias expressed @@ -65,42 +69,65 @@ \item \code{data} dataframe with data from start inclusive till end exclusive. If start == end then data section of final structure is empty. \itemize{ - \item timestamp in numeric format - \item gx gyroscope x-axis (only available in AX6) - \item gy gyroscope y-axis (only available in AX6) - \item gz gyroscope z-axis (only available in AX6) - \item x accelerometer x-axis - \item y accelerometer y-axis - \item z accelerometer z-axis - \item temperature in Celsius (only recorded once per block) - \item battery one value per block (only recorded once per block) - \item light in LUX on (only recorded once per block) + \item time numeric timestamp in seconds since the Unix epoch + \item gx gyroscope x-axis in degrees per second (AX6 with gyroscope enabled) + \item gy gyroscope y-axis in degrees per second (AX6 with gyroscope enabled) + \item gz gyroscope z-axis in degrees per second (AX6 with gyroscope enabled) + \item x accelerometer x-axis in g + \item y accelerometer y-axis in g + \item z accelerometer z-axis in g + \item temp temperature in Celsius (recorded once per block) + \item battery voltage in volts (recorded once per block), or zero when \code{loadbattery = FALSE} + \item light raw light sensor value (recorded once per block), without conversion to lux } \item \code{header} file header. list of header information: \itemize{ - \item uniqueSerialCode is unque serial code of used device - \item frequency is measurement frequency. All data will be resampled for this - frequency - \item start is timestamp in numeric form. To get text representation it is - enough to use as.POSIXct(start, origin = "1970-01-01", tz=desiredtz) + \item uniqueSerialCode is the device serial code + \item frequency is the nominal measurement frequency in Hz, including + fractional rates such as 12.5 Hz. All data are resampled at this frequency + \item start is a \code{POSIXct} timestamp for the first valid data block, + displayed in \code{desiredtz} \item device is "Axivity" \item firmwareVersion - \item blocks is number of datablocks with 80 or 120 raw observations in each. + \item blocks is the number of complete data blocks, including corrupt blocks + \item accrange is the accelerometer range in g (for example, 8 denotes +/-8 g) + \item hardwareType is \code{"AX3"} or \code{"AX6"} + \item blockLength is the sample count of the first valid data block, + typically 40, 80 or 120 depending on the sensor configuration } - \item \code{QClog} dataframe with log of integrity checks on each block including: + \item \code{QClog} dataframe logging failed checksums, imputed intervals and + frequency deviations above 5\%; \code{NULL} when there are no such events. + It is omitted for header-only reads. Columns are: \itemize{ \item checksum_pass Boolean to indicate whether checksum was zero. - \item blockID_previous ID for previous block - \item blockID_current ID for current block - \item start_previous numeric start time of previous block - \item start_current numeric start time of current block - \item blockLengthSeconds length of previous block + \item blockID_current ID of the block whose interval is being checked + \item blockID_next ID of the following block + \item start numeric start time of the checked interval, in seconds since the Unix epoch + \item end numeric end time of the checked interval, in seconds since the Unix epoch + \item blockLengthSeconds duration of the interval in seconds \item frequency_blockheader frequency derived from block header \item frequency_observed frequency as observed - \item imputed Boolean to indicate whether this block was inputed + \item imputed Boolean indicating whether the interval was imputed } } } +\details{ + Only complete 512-byte data packets are decoded. An incomplete final packet + is excluded with a warning when reading the header. An error is raised when + no complete data packets remain. Corrupt packets at the start or end of a + requested interval are skipped up to \code{maxAllowedCorruptBlocks}; internal + corruption can cause the intervening interval to be imputed. + + Acceleration is interpolated using the selected method. Rejected intervals + are filled with the preceding acceleration orientation normalised to 1 g + 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. + + 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 + any imputation between the surrounding valid packets. +} \examples{ \donttest{ cwafile = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread")[1] @@ -111,4 +138,4 @@ Evgeny Mirkes Lena Kushleyeva Vincent van Hees -} \ No newline at end of file +} diff --git a/src/decodeAxivityBlock.cpp b/src/decodeAxivityBlock.cpp index bb677e8..089af23 100644 --- a/src/decodeAxivityBlock.cpp +++ b/src/decodeAxivityBlock.cpp @@ -51,7 +51,9 @@ Rcpp::List decodeAxivityBlock(Rcpp::RawVector block, bool complete, const unsigned int tsOffset = u16(block, 4); const unsigned int lightScale = u16(block, 18); const int axesPacking = block[25]; - const bool packed = (axesPacking & 15) == 0; + const int packing = axesPacking & 15; + if (packing != 0 && packing != 2) stop("Unsupported Axivity sample packing format."); + const bool packed = packing == 0; const int offset26 = s16(block, 26); const int count = u16(block, 28); const bool reuse = parameters.isNotNull(); diff --git a/src/processAxivityBlocks.cpp b/src/processAxivityBlocks.cpp index c48668a..0193cce 100644 --- a/src/processAxivityBlocks.cpp +++ b/src/processAxivityBlocks.cpp @@ -133,8 +133,8 @@ Rcpp::List processAxivityBlocks(Rcpp::Function readBlocks, Rcpp::Function calend // Rcpp::warning() can longjmp past C++ destructors when a handler // exits or warn = 2. Function calls use Rcpp's unwind protection so // the preserved output buffers are released on those paths too. - Function(Environment::base_env()["warning"])( - tfm::format("Skipping corrupt block #%.0f", ii), _["call."] = false); + Function warn = Environment::base_env()["warning"]; + warn(tfm::format("Skipping corrupt block #%.0f", ii), _["call."] = false); if (++consecutiveCorrupt > maxAllowedCorruptBlocks) { stop("Error reading file. %d consecutive blocks are corrupt.", maxAllowedCorruptBlocks + 1); } diff --git a/tests/testthat/test_decodeAxivityBlock.R b/tests/testthat/test_decodeAxivityBlock.R index 01411da..a0adfc4 100644 --- a/tests/testthat/test_decodeAxivityBlock.R +++ b/tests/testthat/test_decodeAxivityBlock.R @@ -179,6 +179,11 @@ test_that("the fractional timestamp shim uses the firmware's integer frequency", test_that("native decoder validates packets before accessing sample payloads", { packet = axivity_packet() expect_error(GGIRread:::decodeAxivityBlock(packet[-1], TRUE, NULL, FALSE), "512 bytes") + packet[26] = as.raw(0x31) # Unknown packing must not be treated as signed 16-bit data. + packet[29:30] = writeBin(1L, raw(), size = 2, endian = "little") + expect_error(GGIRread:::decodeAxivityBlock(axivity_checksum(packet), TRUE, NULL, FALSE), + "Unsupported Axivity sample packing format", fixed = TRUE) + packet = axivity_packet() packet[1] = as.raw(0) expect_identical(GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, FALSE), list(checksum_pass = FALSE)) diff --git a/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index 14e62a7..17cd09b 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -1,6 +1,56 @@ library(GGIRread) context("reading Axivity .cwa data") +test_that("one-sample recordings retain sensor columns and one output row", { + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + for (device in c("ax3", "ax6")) { + fixture = system.file("testfiles", paste0(device, "_testfile.cwa"), + package = "GGIRread") + bytes = readBin(fixture, "raw", n = 1536) + packet = bytes[1025:1536] + packet[29:30] = writeBin(1L, raw(), size = 2, endian = "little") + packet[511:512] = as.raw(0) + checksum = (-sum(readBin(packet, integer(), n = 256, size = 2, + signed = FALSE, endian = "little"))) %% 65536 + packet[511:512] = writeBin(as.integer(checksum), raw(), size = 2, endian = "little") + writeBin(c(bytes[1:1024], packet), filename) + expected = GGIRread:::decodeAxivityBlock(packet, TRUE, NULL, TRUE) + axes = colnames(expected$data) + for (interpolation in 1:2) { + result = readAxivity(filename, start = 0, end = Inf, configtz = "UTC", + loadbattery = TRUE, interpolationType = interpolation) + expect_equal(nrow(result$data), 1L) + expect_identical(names(result$data), c("time", axes, "temp", "battery", "light")) + expect_equal(unname(as.matrix(result$data[, axes])), unname(expected$data)) + expect_equal(result$data$temp, expected$temperature) + expect_equal(result$data$battery, expected$battery) + expect_equal(result$data$light, expected$light) + } + } +}) + +test_that("corrupt lookahead packets respect the consecutive-block limit", { + fixture = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") + bytes = readBin(fixture, "raw", n = file.size(fixture)) + for (block in 5:6) bytes[1024 + block * 512 + 1] = as.raw(0) + filename = tempfile(fileext = ".cwa") + on.exit(unlink(filename)) + writeBin(bytes, filename) + expect_warning( + expect_error(readAxivity(filename, start = 0, end = 5, configtz = "UTC", + maxAllowedCorruptBlocks = 1), + "2 consecutive blocks are corrupt", fixed = TRUE), + "Skipping corrupt end block #5", fixed = TRUE) + expect_warning(actual <- readAxivity(filename, start = 0, end = 5, configtz = "UTC", + maxAllowedCorruptBlocks = 2), + "Skipping corrupt end block #5", fixed = TRUE) + expect_warning(expected <- readAxivity(filename, start = 0, end = 7, configtz = "UTC", + maxAllowedCorruptBlocks = 2), + "Skipping corrupt block #5", fixed = TRUE) + expect_identical(actual, expected) +}) + test_that("incomplete final packets are reported and excluded consistently", { filename = tempfile(fileext = ".cwa") complete = tempfile(fileext = ".cwa") From 4ceed99184c0d6c58c81fc734bae3133cfde97cc Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Fri, 23 Oct 2026 19:12:41 +0200 Subject: [PATCH 14/14] Document AX6 gyro scale codes and code-zero behavior --- src/decodeAxivityBlock.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/decodeAxivityBlock.cpp b/src/decodeAxivityBlock.cpp index 089af23..5f5e52f 100644 --- a/src/decodeAxivityBlock.cpp +++ b/src/decodeAxivityBlock.cpp @@ -118,6 +118,10 @@ Rcpp::List decodeAxivityBlock(Rcpp::RawVector block, bool complete, _["shift"] = -shift / frequency); if (complete) { NumericMatrix data(count, axes); + // Per Daniel Jackson's AX6 clarification, public firmware uses gyro codes + // 2-6 (2000-125 dps); code 0 means disabled gyro in three-axis packets. + // Keep the formula for unused codes as an extension of the range pattern. + // omconvert's code-0 fallback to 2000 dps concerns pre-release prototypes. const double gyroRange = 8000.0 / (1U << ((lightScale >> 10) & 7U)); for (int row = 0; row < count; ++row) { if (packed) {