Repository navigation
Speed up readAxivity with C++ processing and fix final-block reads - #100
vincentvanhees wants to merge 14 commits into
Conversation
Use the final block’s recorded sampling frequency when no preceding block is available to estimate timing, and initialise interpolation state so its samples are returned. Add AX3 and AX6 regression tests for both interpolation modes and corrupt starting blocks that leave only the final valid block.
Move packet decoding, timestamp caching, integrity checks, imputation, and resampling into C++, retaining the R interface and timezone conversion. Add regression tests and a benchmark script excluded from package builds. Real-recording validation produced identical outputs with approximately 4.5× faster reading.
|
@j262byuu @muschellij2 @jhmigueles By the way - I am not planning to suddenly migrate all R code to C++, but in this case I think it is worth it. |
|
I checked it with some data and got this out. Independent check of #100:
|
| Scenario | Files OK in both | Identical | Notes |
|---|---|---|---|
| header only | 15 | 15/15 | |
| full read | 15 | 15/15 | includes 60.5M, 60.1M, 54.7M and 50.5M-row reads |
| first block / last two / beyond EOF | 15 each | 45/45 | |
| interp=2 + battery + tz | 15 | 15/15 | |
frequency_tol = 0.001 (many QClog rows) |
15 | 15/15 | |
| GGIR-style chunked | 15 | 15/15 | |
| last block only (#98) | 0 in base | n/a | base: object 'prevStart' not found; PR: returns data |
| synthetic corrupt/gap/reset/truncated | 18 | 18/18 | same QClog, imputation flags and warning text |
| 21 / 25 consecutive corrupt blocks | both error | same message | 21 consecutive blocks are corrupt |
| corrupt tail, last block only | both error | same message | All requested blocks are corrupt |
Warnings (Skipping corrupt block #…, Skipping corrupt start/end block #…) match exactly in text and count. The QClog$blockID_next column type (integer vs. double when checksum-failure rows exist) also matches.
2. Speed
These are medians of 2 runs per version, run sequentially and interleaved on an otherwise idle machine. "Chunked" means GGIR-style 3000-block chunks with header reused.
| File | Size | Rows | Full read: base → PR | Speedup | Chunked: base → PR | Speedup |
|---|---|---|---|---|---|---|
| 114890_0000000000 | 2.3 MB | 0.55M | 1.34 s → 0.57 s | 2.3× | 1.30 s → 0.51 s | 2.6× |
| tiny-sample | 4.4 MB | 1.07M | 2.51 s → 1.03 s | 2.4× | 2.43 s → 0.80 s | 3.0× |
| example_90001_0_0 | 125 MB | 54.7M | 71.3 s → 29.4 s | 2.4× | 63.3 s → 23.8 s | 2.7× |
| sample | 219 MB | 50.5M | 121.5 s → 42.4 s | 2.9× | 114.7 s → 36.7 s | 3.1× |
| longitudinal_data | 257 MB | 60.1M | 137.1 s → 49.9 s | 2.7× | 129.7 s → 43.5 s | 3.0× |
On the first 5000 blocks with frequency_tol = 0.001 the gain is larger: 1.8–5.8 s → 0.14–1.1 s, or 4–12×. The base reader's repeated rbind() of QClog was the expensive part there.
Rprof of the PR build on a 54.7M-row full read:
processAxivityBlockstakes about 60% of self time.- The remaining R code in
readAxivity(building the outputdata.frame) takes about 15%. strptime/as.POSIXcttake about 15%.
Most of the remaining C++ time is likely decodeAxivityBlock() building and returning an Rcpp::List (plus a named NumericMatrix) for every packet, followed by as<>() look-ups by name. Decoding straight into a C++ struct inside the loop would probably give a further speedup. That is a possible follow-up, not a blocker.
3. Memory
/usr/bin/time -l, full read of example_90001_0_0.cwa (54.7M rows):
| Max RSS | |
|---|---|
| base | 7.57 GB |
| PR | 9.17 GB (+21%) |
This matters mainly for full-file reads of multi-day recordings. GGIR reads in chunks, so it is much less affected.
4. Argument names and behaviour
- Signature: identical:
filename, start = 0, end = 0, progressBar = FALSE, desiredtz = "", configtz = c(), interpolationType = 1, loadbattery = FALSE, header = NULL, frequency_tol = 0.1, maxAllowedCorruptBlocks = 20. - Return value: same structure, names, column order, classes and types.
progressBar = TRUE: identical console output.- New behaviour:
start = n - 1, end = nnow works (readAxivity()` fails with "object 'prevStart' not found" when only the final block is requested #98). With no preceding block, it uses the block's nominal sampling frequency, and the observed frequency is set to the nominal one, so no QClog row is written.- Stricter packet validation in
decodeAxivityBlock(), with newstop()messages:Unsupported number of Axivity axes.(anything other than 3 or 6; for example, 9-axis Gxyz/Axyz/Mxyz packets)Packed Axivity data must have three axes.Invalid Axivity sample count for packet format.Invalid Axivity sampling frequency.
The base reader would instead have errored somewhere downstream or produced misaligned data. None of the 16 files triggered these.
- For very old files (offset 24 = 0),
parameters$formatis now0. Previouslyformatwas never assigned, so the list picked upbase::format, a function. This is internal only. - The block loop now calls
checkUserInterrupt()every 128 packets, so a long read can still be interrupted with Ctrl-C. - Input is read in 128-packet (64 KB) buffers instead of one
readBin()per packet.
5. ⚠️ Build failure under C++11/C++14
My ~/.R/Makevars forces -std=gnu++11, and with that the package fails to install:
processAxivityBlocks.cpp:20:5: error: call to constructor of 'Rcpp::NumericMatrix' (aka 'Matrix<14>') is ambiguous
data(x["data"]) {}
processAxivityBlocks.cpp:47:5: error: ambiguous conversion for functional-style cast from 'Rcpp::Vector<14>' to 'Rcpp::RObject'
processAxivityBlocks.cpp:47:58: error: ambiguous conversion for functional-style cast from 'Rcpp::IntegerVector' to 'Rcpp::RObject'
-fsyntax-only gives 3 errors under gnu++11 and gnu++14 and 0 under gnu++17. R ≥ 4.3 defaults to C++17, so CRAN is probably fine. But DESCRIPTION declares R (>= 3.5.0), whose default standard was C++11/14, and some users pin a standard. The following change compiles cleanly under gnu++11, gnu++14 and gnu++17:
- data(x["data"]) {}
+ data(Rcpp::as<Rcpp::NumericMatrix>(x["data"])) {}
...
- Rcpp::RObject(Rcpp::as<Rcpp::NumericVector>(next)) : Rcpp::RObject(next);
+ Rcpp::RObject(Rcpp::wrap(Rcpp::as<Rcpp::NumericVector>(next))) : Rcpp::RObject(Rcpp::wrap(next));The alternative is to declare SystemRequirements: C++17 and/or raise the R dependency.
6. 🐛 Final-block timing bug: present before this PR, carried over unchanged
When a read of two or more blocks runs up to EOF, the final block of the file is placed on the wrong time axis. In the ii == numDBlocks branch, prevStart is not advanced to the final block's own start. As a result, the final block's samples are spread over the interval from the previous block's start to the extrapolated end. That is twice the real duration, starting one block too early.
Effect on the output covering the final block's time span (about 1.2 s at 100 Hz):
- The first half of the final block's samples falls before that span and never appears in the output.
- The second half is stretched over the whole span at half speed.
I confirmed this exactly. Modelling the final block's raw samples as seq(prevStart_{n-2}, extrapolated_end, length.out = len + 1) reproduces the last_two output with 0 residual on 5 files. Against a correct time axis the residual is 0.1–2.7 g.
- In base R this sits in
rawTime[2:(rawLast+1)] = seq(prevStart, newTimes, length.out = rawLast). The PR keeps it intentionally ("Preserve the established final-packet extrapolation"). - The new
firstIterationpath for readAxivity()` fails with "object 'prevStart' not found" when only the final block is requested #98 setsprevStart = previous.startand is therefore correct. So the same final block now comes out differently depending on whether the chunk contains one block or more. - The impact is small: about one block (~1–2.5 s) at the very end of a recording, usually off-body. It does still affect every full read, and every GGIR chunk that ends at EOF.
A tested one-line fix in src/processAxivityBlocks.cpp. Compute nextTime first, then advance prevStart:
} else {
nextTime = (previous.start - prevStart) / prevLength * previous.length + previous.start;
}
+ prevStart = previous.start;
prevLength = previous.length;With this change, the final block in a last_two read matches the raw packet samples on the expected time axis (spacing estimated from the previous block) to ≤ 1e-5 g on 6 files. Two existing reference tests then shift, because they encode the old behaviour:
test_readAxivity.R:67,:71(AX3):nrow17599 → 17600;floor(sum(abs(…)))25367 → 25405test_readAxivity.R:332,:336(AX6):nrow11428 → 11429; sum 960046 → 960062
Since #100 is otherwise a pure refactor with exact output parity, this may be cleaner as a separate follow-up PR.
Summary
The C++ port matches the current R reader exactly on everything I could throw at it, fixes #98, and is roughly 2.5–3× faster for typical reads. Before merging, I'd suggest:
- Applying the two-line C++11/14 compatibility fix, or declaring C++17.
- Optionally, a follow-up for the final-block
prevStartbug.
The scripts used for this comparison are available on request.
… and memory retention - Use explicit Rcpp conversions for C++11/14 compatibility - Correct the final packet’s time axis in reads ending at EOF - Release native buffers before assembling trimmed output - Add AX3/AX6 regression tests and update reference expectations - Update NEWS and benchmark parity checks for corrected EOF behavior Validation: full test suite passes; compilation checked under C++11/14/17.
|
Data available at:
GGIRread
|
GGIRread readAxivity |
actipy _read_device |
|
|---|---|---|
| Raw x/y/z per packet | identical: all 25 files, every sample, max |diff| = 0 | identical |
| Output samples | resampled to a uniform grid at the header rate | raw device samples, no resampling |
| Gaps / bad blocks | imputed: last orientation normalised to 1 g, logged in QClog |
left as gaps (no rows); counted in ReadErrors / NumInterrupts |
| Timestamps | double-precision epoch seconds; uses the fractional-second field; timezone-aware (configtz → desiredtz) |
datetime64[ns], truncated to ms, timezone-naive device clock; ignores the fractional-second field and chains blocks at the nominal rate |
| Typical sample-time offset (actipy − GGIR) | — | −78 to +27 ms (file medians), max about ±100 ms |
| Temperature | raw * 75/256 − 50 (cwa.h, OM reference reader) |
(raw * 150 − 20500)/1000: 2–16 °C lower |
| Light | raw 10-bit ADC value | 10^(raw/341) (neither raw nor the cwa.h lux formula) |
| Battery | optional (loadbattery = TRUE) |
not returned |
| Gyro (AX6) | 8000/2^n dps; 8000 dps when range code = 0 | 8000/2^n dps; 2000 dps when code = 0 (same as OM's reader) |
| Value precision | double | float32 (acceleration exact; gyro differs by ≤ 1.5e-5 dps) |
| Device serial | 32-bit (upper and lower words) | lower 16 bits only (e.g. 114890 → 49354) |
| Speed, 60M-sample file | 50 s (PR), 137 s (v1.0.11) | 3.1 s raw; 21 s with resample_hz='uniform' |
With resample_hz='uniform' (section 9) |
linear (default) or nearest-neighbour (interpolationType = 2) |
nearest-neighbour onto a grid starting at its first timestamp. Row counts match GGIR to within ±6, but only 10–43% of rows are identical even against GGIR's nearest-neighbour mode, because of the timestamp offsets. Every actipy NaN row lies inside a GGIR-imputed span. |
1. Raw sample decoding: identical
For every file I decoded each packet with GGIRread's decoder and lined the samples up against actipy's output in order.
- Sample counts are identical on all 25 files: same valid blocks, same samples per block.
- x, y, z: max |diff| = 0 in every file. Acceleration values are multiples of 1/256 g, so float32 stores them exactly.
- Gyro (
ax6_testfile.cwa, range code 5, 250 dps): max |diff| = 1.5e-5 dps. That is float32 rounding; the median ratio is 1.00000006. - Packed and unpacked formats decode the same way. The packed 10-bit + exponent unpacking is equivalent in both readers.
Every difference below therefore comes from timing, unit conversion, gap handling or metadata. None comes from decoding acceleration.
2. Output shape and columns
| GGIRread | actipy | |
|---|---|---|
| Container | list(header, data, QClog); data is a data.frame |
(DataFrame, info dict); time is the index |
| AX3 columns | time, x, y, z, temp, battery, light |
x, y, z, temperature, light |
| AX6 columns | time, gx, gy, gz, x, y, z, temp, battery, light |
x, y, z, gyro_x, gyro_y, gyro_z, temperature, light |
| Rows | uniform grid at the header frequency over the whole span, gaps included | one row per recorded sample |
Row counts differ by design. Two examples:
114890_0000000000: actipy has 530,760 raw samples, GGIR has 545,607 rows. The device really sampled at about 48.6 Hz, and GGIR resamples to the nominal 50 Hz.example_90001_0_0: actipy has 29.2M samples, GGIR has 54.7M rows. The recording has 3.45 days of data spread over a 6.33-day span; GGIR fills the roughly 25M-sample gaps by imputation.
With resample_hz = 'uniform', actipy gives 54,686,811 rows on the same file, one fewer than GGIR. 24,890,906 of them are NaN where GGIR imputes values (see section 9).
3. Timestamps
actipy's rule (AxivityReader.java):
- Block start = calendar time −
timestampOffset/rate. The fractional-second field (offset 4, top bit set) is ignored. - If a block starts less than 1 s after the previous block's nominal end (start + n/rate), it is snapped to that end.
- Samples are spaced evenly from there to the block's own nominal end.
- Times are truncated to whole milliseconds.
I emulated this rule in R and reproduced every actipy timestamp to within 1 ms (the residual is Java float32 arithmetic).
GGIR's rule: start = calendar time + fractional/65536 − (offset + fractional-correction)/rate. Samples are spaced evenly up to the next block's start, which is the actual timing.
The difference between them comes from two parts:
- The ignored fractional-second field. This makes GGIR's start 0 to 1 sample period later; the median is 5 ms at 100 Hz and 10 ms at 50 Hz. It is also visible in the first-sample difference between GGIR's
header$startand actipy, which is 0–10 ms. - Chaining at the nominal rate. actipy places block k in the window [start₍ₖ₋₁₎ + n/rate, startₖ + n/rate). The offset therefore equals the device's clock-rate error over one block.
| File | Actual rate | Median offset, actipy − GGIR |
|---|---|---|
| 114890_0000000000 | 48.6 Hz (nominal 50) | −78.0 ms |
| tiny-sample | 97.5 Hz | −34.5 ms |
| 32610_0000000000 | 99.0 Hz | −17.2 ms |
| ax6_testfile | 99.0 Hz | −9.1 ms |
| longitudinal_data | ≈100 Hz | +4.0 ms |
| sample | — | +14.5 ms |
Across the files, offsets range from about −95 to +40 ms. Larger outliers occur only next to real recording gaps.
Consequence for values. I interpolated actipy's raw samples onto GGIR's time grid. The mean |diff| is 0.003–0.13 g, the 99th percentile 0.03–1.4 g, and the maximum up to 16 g. This is almost entirely the timing offset: with GGIR's timeline applied to the same raw values, the gap closes. Vigorous movement magnifies small time shifts.
Other timestamp differences:
- Timezone: actipy has none; it returns the device's local clock as a naive time. GGIR interprets the device clock in
configtzand returns time indesiredtz, so DST transitions are handled. - Precision: actipy truncates to milliseconds. GGIR keeps full double precision.
- Start and end summary: actipy's
StartTimeandEndTimeare strings truncated to seconds. GGIR'sheader$startis a POSIXct.
4. Unit conversions
Temperature. cwa.h and openmovement's own Python reader (cwa_load.py) both give centigrade = raw * 75/256 − 50. GGIR uses this. actipy uses (raw*150 − 20500)/1000. The two agree only at raw ≈ 206 (≈ 10 °C). At body temperatures actipy reads 7–11 °C lower; across all files the difference was 2–16 °C. Example file means: 31.9 °C (GGIR) vs. 21.4 °C (actipy).
Light. GGIR returns the raw 10-bit value (0–1023), as OM's reader does. actipy returns 10^(raw/341) (range 1–1000), confirmed to 6e-8 relative error. This is not the cwa.h lux formula, 10^((raw+512)*6/1024).
Gyro range. Both use 8000/2ⁿ dps. When the range code is 0, GGIR gives 8000 dps and actipy gives 2000 dps; OM's reader also defaults to 2000. Not exercised here, because the fixture uses code 5.
Battery. actipy does not return it.
Sample rate (from code, not exercised).
- GGIR computes the rate as
round(3200 / 2^(15−code)). Because R'sround()rounds half to even, a 12.5 Hz recording is treated as 12 Hz, and 6.25 / 3.125 Hz become 6 and 3. actipy keeps 12.5 as a float. - actipy's
info['SampleRate']comes from the last valid block. GGIR'sheader$frequencycomes from the file header, with a warning if the first block differs.
5. Metadata
GGIRread header |
actipy info |
|---|---|
uniqueSerialCode (32-bit: upper word at offset 11, lower at 5) |
DeviceID: lower 16 bits only. 114890 → 49354; the AX6 fixture 6011834 → 48058 |
frequency (header) |
SampleRate (last block, float) |
start (POSIXct, desiredtz) |
StartTime, EndTime (strings, seconds; from read_device) |
device = "Axivity", hardwareType (AX3 / AX6), firmwareVersion, accrange, blocks, blockLength |
Device = "Axivity" only (no hardware type, firmware or range) |
| — | Filename, Filesize(MB) |
QClog (per event: block IDs, start/end, observed vs. header frequency, imputed flag) |
ReadOK, ReadErrors, NumTicks, WearTime(days), DataSpan(days), NumInterrupts, Covers24hOK |
6. Corrupt, gapped and unusual files
| Case | GGIRread | actipy |
|---|---|---|
| 1 or 5 corrupt blocks mid-file | skipped, gap imputed, QClog rows, warning per block |
skipped, gap left, ReadErrors = 1 / 5, NumInterrupts +1 |
| 25 consecutive corrupt blocks | error: 21 consecutive blocks are corrupt (maxAllowedCorruptBlocks = 20) |
reads; ReadErrors = 25 |
| 100 blocks removed (time gap) | imputed, QClog imputed = TRUE |
gap left |
Corrupt tail (3 blocks; also GGIRread's own ax3_testfile_corrupt_*) |
skips them and drops the last valid block: ends 1.2 s earlier than actipy (15:19:29.044 vs. 15:19:30.249) | keeps the last valid block |
| Truncated final packet | ignored silently (block count = floor) | stops; ReadErrors = 1, message on stderr |
left_wrist_data.cwa (2011; 254 blank blocks, 16 SI blocks, 240 blank, data, 256 zero blocks, 228 blank, data) |
fails: first 21 blocks are corrupt. With maxAllowedCorruptBlocks = 2000: wrong sign in 'by' argument |
reads 7.70M samples. Non-AX blocks are skipped without counting as errors; quality_control removes 204,799 non-increasing timestamps (ReadErrors = 373); DeviceID = 0 |
Packet without the AX marker but with a valid checksum |
stop("Packet header is incorrect") |
skipped silently |
GGIR's imputation also overwrites the valid block just before each gap. The QC check runs on the previous block, from its start to the next valid block's start, so a gap fails the frequency check and that block's real samples are replaced by imputed ones.
Example, corrupt_one.cwa (block 100 corrupt): the imputed span is 2.46 s and covers valid block 99 plus corrupt block 100. actipy keeps 117 real samples from block 99 in that span. Section 9 shows this costs about 120–245 extra imputed rows per gap.
7. A GGIR resampling detail
GGIR's output is exactly linear interpolation of the raw samples on its own timeline, with one exception. Between a block's last sample and the next block's first sample, GGIR interpolates from its previous output value (the "anchor") instead of from the last raw sample. On 3 test files this changed 0.5% of output points, by up to 1.5 g in active data. All other non-imputed points matched to < 1e-9 g. The behaviour is the same in v1.0.11 and the PR. It is separate from the final-block timing issue reported on #100.
8. Speed and memory
Raw-read columns are medians of 2 runs; the uniform column is a single run. actipy times include JVM start-up, about 0.15 s.
| File | Samples | actipy _read_device |
actipy read_device (minimal) |
actipy resample_hz='uniform' |
GGIR PR | GGIR v1.0.11 |
|---|---|---|---|---|---|---|
| 114890_0000000000 | 0.53M | 0.23 s | 0.18 s | 0.34 s | 0.57 s | 1.34 s |
| tiny-sample | 1.02M | 0.20 s | 0.25 s | 0.46 s | 1.03 s | 2.51 s |
| example_90001_0_0 | 29.2M | 1.52 s | 3.03 s | 14.3 s | 29.4 s | 71.3 s |
| sample | 51.4M | 2.62 s | 5.40 s | 17.4 s | 42.4 s | 121.5 s |
| longitudinal_data | 60.3M | 3.05 s | 5.96 s | 21.0 s | 49.9 s | 137.1 s |
The uniform column is the closest like-for-like comparison, since GGIR always resamples (and also imputes). On that basis actipy is about 2–2.4× faster than the PR build on large files.
Peak RSS on example_90001_0_0:
- actipy, raw: 1.0 GB
- actipy, uniform resampling: 4.5 GB
- GGIR PR: 9.2 GB (double-precision
data.framewith 7 columns)
9. actipy with resample_hz='uniform'
Settings: actipy.read_device(f, lowpass_hz=None, calibrate_gravity=False, detect_nonwear=False, resample_hz='uniform'). GGIR was run with desiredtz="UTC" and both interpolationType = 1 (linear, the default) and interpolationType = 2 (nearest neighbour). Each actipy grid point was paired with the nearest GGIR grid point.
How the two resamplers differ
| GGIRread | actipy resample |
|
|---|---|---|
| Method | linear (default) or nearest neighbour | nearest neighbour only (get_indexer(method='nearest')) |
| Grid start | first valid block's start, including the fractional-second field | first raw timestamp: ms-truncated, fractional ignored |
| Step | 1 / header frequency | 1 / SampleRate (taken from the last block) |
| Grid length | up to the last block's extrapolated end | round((t_last − t_first) * rate) + 1 |
| Gaps | imputed with last orientation (normalised to 1 g, else 0,0,1); temperature and light carried forward | NaN wherever the nearest raw sample is more than 1 s away; temperature and light are NaN there too |
| Values | double | float32 |
| Other | — | skipped if the data already has a perfectly uniform rate; optional start_first_complete_minute |
Results (selected files; all 25 were run)
| File | Rows actipy / GGIR | Grid start, actipy − GGIR | Rows identical vs. GGIR nearest | Mean |Δ| vs. GGIR linear / nearest | Mean ENMO, mg: actipy / GGIR nearest / GGIR linear |
|---|---|---|---|---|---|
| 114890_0000000000 (50 Hz) | 545,605 / 545,607 | −1.9 ms | 13.9% | 0.053 / 0.054 g | 26.44 / 26.45 / 25.44 |
| 31447_0000000000 | 227,313 / 227,314 | −2.7 ms | 21.5% | 0.128 / 0.131 g | 116.53 / 116.61 / 113.34 |
| 32816_0000000000 | 216,322 / 216,320 | −2.5 ms | 43.2% | 0.042 / 0.043 g | 273.03 / 272.80 / 269.82 |
| example_90001_0_0 | 54,686,811 / 54,686,812 | −1.1 ms | 29.0% | 0.025 / 0.025 g | 58.37 / 58.37 / 58.01 |
| longitudinal_data | 60,110,388 / 60,110,386 | 0.0 ms | 76.0% | 0.0025 / 0.0026 g | 20.16 / 20.16 / 20.10 |
| sample | 50,524,347 / 50,524,343 | −9.3 ms | 33.8% | 0.025 / 0.026 g | 44.29 / 44.29 / 43.55 |
| tiny-sample | 1,070,944 / 1,070,945 | −4.3 ms | 16.5% | 0.027 / 0.028 g | 25.93 / 25.91 / 25.40 |
| ax6_testfile (AX6) | 11,430 / 11,428 | −9.8 ms | 10.0% (gyro 10.0%) | 0.106 / 0.108 g | 617.87 / 618.72 / 610.81 |
ENMO and the identical-row percentages are computed on rows where both outputs are valid (actipy not NaN, GGIR not imputed).
Findings
- The grids line up; the samples don't.
- Row counts agree to within ±6 everywhere except the corrupt-tail files. There actipy has 121 more rows, because GGIR drops the last valid block (section 6).
- actipy's grid starts 0–10 ms earlier (the ignored fractional second) and keeps a constant sub-sample phase relative to GGIR's grid.
- Even against GGIR's nearest-neighbour mode, only 10–43% of rows are identical. The exception is 76% for
longitudinal_data, whose clock runs almost exactly at 100 Hz. - The mean |Δ| (0.003–0.13 g, max 16 g) is unchanged from the raw comparison. The cause is actipy's timestamp model (section 3), not the resampler: shifts of tens of ms make nearest-neighbour pick a raw sample 1–8 positions away.
- ENMO depends on the interpolation method, not the reader. actipy (nearest) matches GGIR's nearest-neighbour mode to within 0.25 mg on every file except the AX6 fixture (0.85 mg). GGIR's default linear interpolation gives 0.06–3.3 mg (0.3–8%) lower mean ENMO than either, because it smooths peaks.
- Gaps. Across all 25 files, every actipy NaN row falls inside a GGIR-imputed span, with zero exceptions. GGIR imputes more rows than actipy leaves empty: 240–5,601 extra per file, about 120–245 per gap. These extra rows are mostly the valid block before each gap that GGIR overwrites (section 6), where actipy keeps real samples.
- Summaries over all rows diverge on gappy files. GGIR's imputed rows count as near-zero ENMO, while actipy's NaN rows drop out. For
example_90001_0_0, mean ENMO is 31.60 mg in GGIR (all rows) vs. 58.37 mg in actipy (NaN removed). - Files GGIR cannot read (
left_wrist_data.cwa,corrupt_25mid.cwa) resample fine in actipy: 7,667,354 rows with 31,343 NaN, and 1,070,944 rows with 27,337 NaN respectively. - The temperature, light, serial and timezone differences from sections 3–5 still apply. Resampling does not change them.
Bottom line
Both readers decode exactly the same acceleration and gyro samples. Everything that differs comes after decoding:
- actipy keeps the raw samples and gaps; GGIR resamples and imputes.
- Even with
resample_hz='uniform', the outputs share a grid (to within ±6 rows) but not values: only 10–43% of rows are identical against GGIR's nearest-neighbour mode. GGIR's default linear interpolation also gives 0.3–8% lower ENMO than nearest-neighbour. - actipy's NaN rows always fall inside GGIR's imputed spans. GGIR also imputes over the valid block before each gap.
- Timestamps differ by tens of milliseconds, because actipy ignores the fractional-second field and chains blocks at the nominal rate.
- Temperature and light are in different units, and actipy's temperature formula disagrees with the cwa.h spec.
- actipy truncates the device serial to 16 bits.
- GGIR stops on unusual or heavily corrupted files that actipy reads.
Preserve R’s rounding when constructing native timestamps Resolve interpolation test failures without relaxing tolerances Update checkout actions to v5 and document the fix in NEWS
Use accelerometer channels for orientation normalization and the 1 g fallback, and fill gyroscope channels with zero. Add regression tests for both cases with linear and nearest-neighbour interpolation, and update NEWS.md.
Treat the last valid packet as the effective end of the recording so it is processed without lookahead, including single-packet reads. Add AX3 and AX6 regression coverage, update existing tests and NEWS.md.
Remove integer rounding from header and packet rate decoding so rates such as 12.5 Hz are retained in timestamp corrections and the uniform output grid. Add regression coverage and update NEWS.md.
Count only complete packets, warn about ignored trailing bytes, and report an error when no complete data packets remain. Use nominal spacing for single final blocks to prevent floating-point rounding from dropping the last sample at fractional sampling rates. Add regression tests and update NEWS.md. Full suite: 2,780 tests pass.
Use the preceding raw sample and timestamp instead of the last resampled output for linear and nearest-neighbour interpolation. Preserve the existing pre-gap imputation policy. Add regression coverage and update NEWS.md. Full suite: 2,802 tests pass.
|
Thanks, John, for the detailed comparison. We investigated the findings and made the following changes to Fixes following the comparison
Earlier fixes in this PR
Both complete AX3 and AX6 fixtures matched independent linear interpolation through their decoded raw samples exactly. The full package suite passes locally: 2,802 assertions, zero failures. Temperature conversionGGIRread uses the documented MCP9700 conversion:
Daniel previously confirmed that AX devices use the MCP9700, and that the alternative formula, I verified that actipy’s Java reader at 4eb2bea—the revision used in the comparison and the current This is a difference in sensor conversion, rather than units. No temperature-formula change is needed in GGIRread. Gyro range code zeroThe manufacturer’s implementations disagree:
For acceleration-only packets, the gyro scale code has no effect. The unresolved case is a six-axis @danielgjackson, could you clarify whether this combination can occur in valid recordings, particularly with older firmware, and what scaling a reader should apply? We are unsure whether zero should indicate a missing scale value, a disabled gyro, or something else. We have left GGIRread’s behavior unchanged pending clarification. Findings left unchanged
These fixes should reduce some discrepancies with actipy and OpenMovement by correcting identifiable errors in GGIRread, although matching their output was not the objective. We have not yet repeated the comparison to quantify the remaining differences, and the readers retain different timestamp models and gap-handling policies. |
|
I double-checked that the latest changes address the C++ version issue, and I run all the Change 1: the timestamp shim should use the firmware's integer frequency (bug)Commit Effect: at 12.5 Hz, when the fractional second is above about 0.96, the block start shifts by one sample (80 ms). I confirmed this with a crafted 12.5 Hz packet: the PR gives File: Before: if (format == 1) {
// Undo the firmware's backwards-compatible timestamp shim: the lower
// 15 bits hold a fraction of a second, and offset26 includes its samples.
fractional = (tsOffset & 0x7fffU) << 1;
shift += std::floor(fractional * frequency / 65536.0);
}After: 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);
}For integer rates (25 Hz and up) the result is unchanged, so no existing test changes. Regression test. File test_that("the fractional timestamp shim uses the firmware's integer frequency", {
packet = axivity_packet()
packet[25] = as.raw(0x47) # 12.5 Hz
# Fraction 0.96875 s: the firmware shifts floor(0.96875 * 12) = 11 samples,
# whereas the exact rate would give floor(0.96875 * 12.5) = 12.
packet[5:6] = writeBin(0xfc00L, raw(), size = 2, endian = "little")
packet[27:28] = writeBin(-5L, raw(), size = 2, endian = "little")
decoded = GGIRread:::decodeAxivityBlock(axivity_checksum(packet), FALSE, NULL, FALSE)
expect_equal(decoded$frequency, 12.5)
expect_equal(decoded$shift, -(-5 + 11) / 12.5)
})Change 2: warn about trailing bytes only once, not for every chunkThe new "Ignoring N trailing bytes" warning fires on every call. GGIR reads a file in chunks and passes A simple rule: warn only when the header is being read, which in GGIR happens once per file. File: Before: if (trailingBytes > 0) {After: if (trailingBytes > 0 && is.null(header)) {The PR's own test expects the warning on cached-header reads, so it needs to change too. File: Before: expect_warning(actual <- readAxivity(filename, start = start, end = Inf,
configtz = "UTC", header = expected$header,
interpolationType = interpolation),
message)After: # The warning is given only when the header is read, not for every chunk.
expect_no_warning(actual <- readAxivity(filename, start = start, end = Inf,
configtz = "UTC", header = expected$header,
interpolationType = interpolation)) |
|
Thanks @vincentvanhees I tested the newest pr with 10 institutional Ax3 CWA files (7 day, very similar protocol as UK Biobank, sorry I can't share those files) and all previously mentioned files under both Win 10 (R 4.6.1 & Rtools45) and Linux (docker) Memory usage is a recurring issue on our institutional servers, so I paid particular attention to peak mem/thread use and potential leaks. What I find 1 C++ standards:
2 Full reads are 4.3-5.8x faster on Linux:
3 I ran GGIR 3.3.6 parts 1-2 with default settings on 4 of the institutional ax3 files, once with each GGIRread build: Epoch timestamps, non-wear scores and the part 2 row structure are identical 4 Memory leak under exiting warning handlers (I would suggest fix before merge) Problem.
Suggested fix. Use base R's change GGIRread/src/processAxivityBlocks.cpp Line 133 in 820e113 to Function(Environment::base_env()["warning"])(tfm::format("Skipping corrupt block #%.0f", ii), _["call."] = false);Please also require Reproduction Use corrupt_one.cwa (1.07M rows, one bad block), 30 calls increased R mem usage by 1723 MB on the PR vs 4.7 MB with the suggested fix. I reproduced this with both an exiting Scope I did not see this leak under GGIR's default warning handling. I reproduced it with exiting warning handlers and Ref Rcpp #1474 and #1475 address the related |
Use the firmware's truncated integer frequency to undo the timestamp shim while preserving fractional rates for sample timing. Warn about trailing bytes only when reading the header. Add regression tests and update NEWS.
Emit corrupt-block warnings through Rcpp's protected call to base R so output buffers are released when warning handlers exit or warn = 2. Require Rcpp >= 1.0.10, add memory regression tests and update NEWS.
Preserve sensor columns in single-sample reads, enforce corrupt-block limits during endpoint searches and reject unsupported packing formats. Use a named warning callback to fix macOS compilation while retaining unwind protection. Update documentation, NEWS and regression tests.
|
Thanks everyone for the detailed feedback and testing. I have addressed the findings in three commits:
I am not treating differences from the previous GGIRread version as regressions by themselves. We have corrected known bugs, so some changes in output are inevitable. Comparisons against manufacturer specifications and independently established expected results remain valuable. I plan to do more testing over the coming days. If no further problems emerge, I intend to merge this PR next week. |
@vincentvanhees No, this combination should not be found in the wild. The bit patterns, as produced by the AX6:
The missing combinations (8000/4000/62.5 dps) can be harmlessly treated as extending the pattern, and I think the scaling was specifically chosen to allow the range to be extended either side, if required in the future.
|
Move Axivity packet decoding and the block-processing loop into C++, including timestamp caching, integrity checks, imputation, and resampling. Preserve the existing R interface and timezone conversion, and use bounded input buffering for large recordings.
Fix
readAxivity(start = n - 1, end = n)failing with object 'prevStart' not found. When a chunk contains only the final block, initialise interpolation state and use that block’s recorded sampling frequency. This also handles corrupt starting blocks being skipped, leaving only the final valid block. Fixes #98Add regression tests for decoding, buffer boundaries, corruption, device resets, calendar changes, and final-block reads. Include a repeatable benchmark script, excluded from package builds.
Validation
One existing corrupt-tail case now reports an explicit error instead of failing because last was undefined.
Checklist before merging