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/.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/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 83637cd..de426b4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,35 @@ +# Changes in version 1.0.12 (release date:??-??-2026) + +- 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. + +- 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 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. + +- 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 + +- 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. + +- 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/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 acf83f4..92cff46 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)) } @@ -275,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. # @@ -306,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 @@ -360,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 && is.null(header)) { + 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 ############################################################################# @@ -413,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) } @@ -422,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') @@ -446,6 +280,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 @@ -455,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 @@ -489,213 +326,41 @@ 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 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 - 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]) - #------------------------------------------------------------ - # 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) + # 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. 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/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..5f5e52f --- /dev/null +++ b/src/decodeAxivityBlock.cpp @@ -0,0 +1,147 @@ +#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 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(); + 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); + frequency = samplerate == 0 ? offset26 : + 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. + // 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 * std::floor(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); + // 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) { + 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..0193cce --- /dev/null +++ b/src/processAxivityBlocks.cpp @@ -0,0 +1,270 @@ +#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(Rcpp::as(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::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, + 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"]); + // 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]); + 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"]; + // 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); + 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; + 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. + if (firstIteration) { + prevStart = previous.start; + nextTime = prevStart + previous.length / previous.frequency; + } else { + // 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; + } 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}); + // 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 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); + } + 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); + // 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. + // 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); + anchorTime = prevStart - 0.00001; + interpolationAnchor = anchor; + interpolationTime = anchorTime; + firstIteration = false; + } + 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 && + 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[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[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 < accelOffset ? 0 : imputed[axis - accelOffset]; + } + } 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 ? 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; + } + } + 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); + 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); + } + 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..a0adfc4 --- /dev/null +++ b/tests/testthat/test_decodeAxivityBlock.R @@ -0,0 +1,204 @@ +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 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) + 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("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") + 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)) + 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..135a68d --- /dev/null +++ b/tests/testthat/test_processAxivityBlocks.R @@ -0,0 +1,310 @@ +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)) + 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 + 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. + fixture = system.file("testfiles/ax3_testfile.cwa", package = "GGIRread") + bytes = readBin(fixture, "raw", n = 1536) + packet = bytes[1025:1536] + 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) * 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) + + 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") + 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) * count + expect_equal(nrow(result$data), samples) + 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)) + 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("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("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", { + 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/tests/testthat/test_readAxivity.R b/tests/testthat/test_readAxivity.R index 8137dea..17cd09b 100644 --- a/tests/testthat/test_readAxivity.R +++ b/tests/testthat/test_readAxivity.R @@ -1,5 +1,223 @@ 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") + 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) + # 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)) + 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)) + 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"), + 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("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) @@ -16,11 +234,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) }) @@ -220,10 +438,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), @@ -235,17 +457,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. @@ -281,11 +497,12 @@ 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) + # Includes interpolation through raw samples at packet boundaries. + expect_equal(floor(sum(abs(AX6b$data[,2:4]))), 960636) expect_true(is.null(AX6b$QClog)) }) diff --git a/tools/benchmark-readAxivity.R b/tools/benchmark-readAxivity.R new file mode 100644 index 0000000..992a6ca --- /dev/null +++ b/tools/benchmark-readAxivity.R @@ -0,0 +1,82 @@ +# 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({ + # 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 - 1, 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)) +}