Data Analysis - Prototype M.gigas Gaping Pilot Project Round 2

2026
Author

Sam White

Published

October 19, 2026

INTRO

On 20261008 I deployed our second gape logger prototype in the lab (the first one got fried).

The primary goal here was to get an actual extended data set, as well as to assess how long the battery might last.

Here’s an example of how it looked (not an image from this round):

Gape monitor prototype. On the left is a closed, waterproof Pelican Case with yellow insert, a blue battery pack, and an Adafruit Feather Adalogger RP2040 microcontroller, connected to a real-time clock module via wires. On the right, the wire connects to an Adafruit 3-axis magnetometer. It is sealed in a small plastic bag with white marine sealant/adhesive. The bag is attached to the top valve of a Pacific oyster using that same marine sealant/adhesive. The green LED on the magnetometer is illuminated. The oyster and magnetometer are in a small plastic container with seawater.

MATERIALS & METHODS

Important

This was prototype Unit A. This is important because the analysis requires the correct calibration info for the corresponding unit magnetometer!

Gape monitor prototype consists of the following:

  • a 3.7V 6600mAh lithium ion battery pack (Adafruit)

  • Feather Adalogger RP2040 microcontroller with 16GB microSD card (Adafruit)

  • real-time clock (RTC) breakout board (Adafruit; DS3231)

  • Three-axis magnetometer (Adafruit; TLV493D).

Connections between components are I2C STEMMA QT. This type of connection is limited to a max of 400mm between components and severely limits the range for deployment. This is definitely a serious shortcoming of this prototype.

The magnetometer was sealed in a small plastic bag with white 3M marine adhesive sealant. The bag was attached to the top valve of a Pacific oyster using that same 3M marine adhesive sealant.

Three magnets (4mm x 2mm neodymium) were secured near the leading edge of the bottom valve.

Note

The sealant for the magnetometer and the magnets was allowed to cure at room temp in the fume hood for ~2hrs prior to returning the oyster to water.

The oyster was kept in a small plastic container with Instant Ocean at room temp (~10C) in FTR 210 (this is a cold room). It did not get fed, nor have water flow-through.

I did change the water around 12PM on 10/16/2026. I used water from an existing tank of animals, instead of straight Instant Ocean.

RESULTS

Prototype Performance

The logger ran continuously from 2026-10-08 12:30 to 2026-10-19 01:25 (~10.5 days), recording 910,493 samples at 1 Hz with no gaps longer than 2 s. The RTC supplied absolute timestamps for the entire deployment. The battery (3.7V 6600mAh) was depleted after this duration.

The main limitation is the magnetometer’s resolution: the TLV493D reports each axis in 98 µT steps, so raw readings jump by ~100 µT from sample to sample. The analysis smooths the signal (15 s rolling median) before classifying open/closed, which cleanly separates closures from openings.

Gaping Behavior

Percent opening over time (1-minute means) for the ~10.5 day deployment. 0% = closed baseline, 100% = widest opening seen. Through Oct 14, the oyster mostly sat around 30-50% open, interrupted by brief full closures near 0%. From Oct 14-15 onward, it opened much wider (70-100%), still with periodic full closures. A red dashed line marks the open/closed detection threshold at ~11%.
  • Mostly open. The oyster was open ~85% of the recording. Closures were distinct and complete (dropping to near 0% opening) rather than gradual.
  • 82 gape events. Median event duration was ~2.7 h; the longest was ~15.5 h (Oct 16 19:23 - Oct 17 10:52). Closures between events were usually short (median ~2.6 min), with the longest lasting ~6 h.
  • Opening increased mid-deployment. Through Oct 13, the oyster sat at ~30-45% opening when open. Starting Oct 14-15, it opened much wider (~85% when open) and stayed that way through the end of the recording. This shift started before the water change on 10/16 (~12 PM), so the water change doesn’t explain it. It may reflect stress from no feeding/no water flow.
  • Gape distance. Calibrated magnet-to-sensor distance ranged from ~14.2 mm (closed) to ~16.4 mm (widest opening, 1st percentile of field magnitude), a swing of ~2.2 mm at the magnets’ position near the leading edge.

OVERVIEW

This notebook reads the CSV log(s) produced by the CircuitPython datalogger running on the Feather RP2040 Adalogger (gapelog_001.csv, gapelog_002.csv, … on the SD card), reconstructs real-world timestamps, detects gape events, and plots the resulting time series.

Data format on disk (wide form): one row per sample, with the three simultaneous magnetic field components as columns:

column meaning
timestamp YYYY-MM-DDTHH:MM:SS from the onboard DS3231 RTC, or T+123.456s (seconds since that session started) if the RTC wasn’t detected at boot
rp2040_uid the Adalogger’s factory-burned unique ID - used to pick that board’s calibration curve
x_uT magnetic field, X axis, microtesla
y_uT magnetic field, Y axis, microtesla
z_uT magnetic field, Z axis, microtesla

Wide form is kept on the device and on disk because X/Y/Z are simultaneous components of one field vector, not independent long-form observations - computing field magnitude, per-axis deltas, etc. all need all three values from the same row. We pivot to long form only for specific plots below.

One file per power-on. The logger starts a new numbered file every time the board boots, so each gapelog_NNN.csv is exactly one session. Session identity comes from the filename - no need to infer reboots from the data.

Onboard real-time clock. The board has a battery-backed DS3231 RTC, so timestamp is normally an absolute wall-clock value straight off the device - no manual bookkeeping needed. If the RTC isn’t detected at boot (and that failure is recorded in eventlog.csv on the card, see below), the script falls back to a relative T+123.456s marker for that whole session, and this notebook derives elapsed_s from either form so the downstream analysis code is the same regardless.

SETUP

Libraries

Analysis parameters

All user-configurable values are assigned here - edit this chunk for your deployment and nothing else below should need changing.

Code 
# ---- Data location -------------------------------------------------------
# Directory holding the CSV files copied off the SD card, and the filename
# pattern to match. Every matching file is treated as one logging session.
# `post_dir` (set above, in the libraries chunk) is this document's own
# folder as an absolute path - reliable regardless of the working directory
# chunks happen to run with.
data_dir      <- post_dir
file_pattern  <- "^gapelog_[0-9]+\\.csv$"

# Timezone the onboard RTC was set in (see set_rtc.py) - used to parse the
# absolute timestamps it writes.
deployment_tz <- "America/Los_Angeles"

# ---- Gape event detection -----------------------------------------------
# The magnet sits on one valve, the magnetometer on the other: as the shell
# OPENS the magnet moves away and field magnitude FALLS. So a gape event is
# a sustained DROP in magnitude below a threshold.
#
# The threshold can be set three ways:
#   "noise"    - computed per session from the sensor noise: threshold =
#                closed_mean - n_sd_threshold * noise_sd, where noise_sd is
#                the noise left in the smoothed magnitude. This ties
#                the cutoff to how noisy that particular session's signal
#                actually is, so a "gape" is a drop that's statistically
#                distinguishable from background noise, not an arbitrary
#                fractional drop. RECOMMENDED DEFAULT.
#   "relative" - computed per session as a fixed fractional drop below the
#                closed-shell baseline value (kept for comparison /
#                backward compatibility - doesn't account for how noisy the
#                signal is).
#   "absolute" - a fixed microtesla value you supply.
threshold_method <- "noise"

# Used for both "noise" and "relative" methods: the quantile of magnitude
# taken as representative of "fully closed" (magnet closest to sensor).
baseline_quantile <- 0.90

# Used when threshold_method == "noise":
# The closed-baseline mean is estimated from smoothed samples at/above
# `baseline_quantile`; the noise SD from lagged differences of the smoothed
# signal (see the classify-state chunk). n_sd_threshold sets how many noise
# SDs below the closed mean a smoothed sample must fall to count as "open" -
# larger values are more conservative (fewer false-positive events from
# sensor noise), smaller values are more sensitive. 3 SD is a common
# starting point for "clearly not noise".
n_sd_threshold <- 3

# Used when threshold_method == "relative":
gape_drop_frac <- 0.10   # open if magnitude < baseline * (1 - this)

# Used when threshold_method == "absolute":
gape_threshold_uT <- 40

# ---- Smoothing ----------------------------------------------------------
# Rolling-median window (seconds) applied to field magnitude before it is
# compared to the threshold. The TLV493D reports each axis in 98 uT steps, so
# single 1 Hz samples jump by ~100 uT from quantization alone; a short
# rolling median suppresses that without blurring real gapes, which last
# tens of seconds or more. Set to 1 (or less) to disable smoothing.
smooth_window_s <- 15

# ---- Percent opening ----------------------------------------------------
# Each session is rescaled so 0% = closed baseline and 100% = widest opening
# seen. "Widest" is this low quantile of smoothed magnitude rather than the
# single minimum, so one outlier sample can't set the scale.
open_ref_quantile <- 0.01

# ---- Event filtering ----------------------------------------------------
# Reject blips and bridge momentary re-closures, so sensor noise doesn't
# inflate the event count.
min_event_duration_s <- 30   # discard "open" runs shorter than this
min_gap_between_s    <- 30   # merge two events separated by less than this

# ---- Gape distance calibration -------------------------------------------
# Convert field magnitude to an actual gape distance in mm, using the curves
# fit by the per-unit calibration notebooks. Set FALSE to skip it entirely -
# every field-magnitude-based result below still works without calibration.
apply_calibration <- TRUE

# Coefficient CSV(s) written by `calibration_unit-*.Rmd`. Pass as many as you
# like (e.g. all three units); the right row is selected per session by
# matching rp2040_uid AND magnet count. Built from `data_dir` for the same
# reason `data_dir` itself is an absolute path - see the libraries chunk.
coef_path <- file.path(data_dir, "calibration-coefficients_Unit-A.csv")

# Deployment metadata: which magnet count each log file was recorded with.
# The logger records the board UID but NOT the magnet count, so this is the
# only place that association exists.
#
# The general version of this (see `read-deployment-metadata` below) reads a
# `deployment_metadata.csv` with one row per log file - built for a real
# deployment with many oysters/boards/sessions, where magnet count can
# differ file to file. This pilot is a single oyster on a single log file,
# so that whole CSV collapses to one value: just set it here.
#
# If this notebook ever covers more than one oyster/magnet-count at once,
# switch `read-deployment-metadata` back to reading a CSV keyed by
# `source_file` instead of broadcasting these two values to every file.
n_magnets_used  <- 3          # magnets glued to this oyster's valve
pilot_oyster_id <- "oyster_1" # cosmetic label only - not used in any calculation

# Which fitted model to apply: "best" (lowest AIC per magnet count, as marked
# in the coefficient CSV), or force "offset" / "power" to compare them.
calibration_model <- "best"

DATA IMPORT

Locate log files

Code 
log_files <- list.files(
  path = data_dir, pattern = file_pattern, full.names = TRUE
) %>%
  sort()

stopifnot("No log files matched `file_pattern` in `data_dir`." =
            length(log_files) > 0)

tibble(
  file     = basename(log_files),
  size_kb  = round(file.size(log_files) / 1024, 1)
) %>%
  kable(caption = "Log files found")
Log files found
file size_kb
gapelog_083.csv 55821.8

Read and combine

Code 
raw <- log_files %>%
  set_names(basename(.)) %>%
  map(
    read_csv,
    col_types = cols(
      timestamp  = col_character(),
      rp2040_uid = col_character(),
      x_uT       = col_double(),
      y_uT       = col_double(),
      z_uT       = col_double()
    )
  ) %>%
  list_rbind(names_to = "source_file") %>%
  mutate(session_id = dense_rank(source_file), .after = source_file) %>%
  rename(timestamp_raw = timestamp)

glimpse(raw)
Rows: 910,493
Columns: 7
$ source_file   <chr> "gapelog_083.csv", "gapelog_083.csv", "gapelog_083.csv",…
$ session_id    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ timestamp_raw <chr> "2026-10-08T12:30:50", "2026-10-08T12:30:51", "2026-10-0…
$ rp2040_uid    <chr> "df6470a31b8a2f2e", "df6470a31b8a2f2e", "df6470a31b8a2f2…
$ x_uT          <dbl> -98, -196, 0, -196, 0, -196, 0, 0, 0, 0, -98, -98, 0, 0,…
$ y_uT          <dbl> -196, -196, 0, 0, -196, -196, 0, -196, -98, -196, 0, -19…
$ z_uT          <dbl> 0, -196, 0, -196, 0, -98, 0, 0, 98, -98, 0, 98, 0, -196,…

Reconstruct timestamps

timestamp_raw holds either an absolute RTC reading (YYYY-MM-DDTHH:MM:SS) or, when the DS3231 wasn’t detected at boot, a relative marker (T+123.456s = seconds since that session started - see eventlog.csv on the card for the matching RTC_INIT_ERROR row). A session is treated as absolute-time only if every row in it parsed as a real RTC reading; otherwise elapsed_s is derived from the relative marker instead, so every downstream chunk below (all written in terms of elapsed_s) works the same either way.

Code 
n_sessions <- n_distinct(raw$session_id)

is_relative <- str_starts(raw$timestamp_raw, "T\\+")

gaping <- raw %>%
  mutate(
    abs_time = if_else(
      is_relative, as.POSIXct(NA, tz = deployment_tz),
      ymd_hms(timestamp_raw, tz = deployment_tz, quiet = TRUE)
    ),
    relative_s = if_else(
      is_relative,
      as.numeric(str_match(timestamp_raw, "^T\\+([0-9.]+)s$")[, 2]),
      NA_real_
    )
  ) %>%
  group_by(session_id) %>%
  mutate(
    timestamp = if (all(!is.na(abs_time))) abs_time else as.POSIXct(NA, tz = deployment_tz),
    elapsed_s = if (all(!is.na(abs_time))) {
      as.numeric(difftime(abs_time, min(abs_time), units = "secs"))
    } else {
      relative_s - min(relative_s, na.rm = TRUE)
    }
  ) %>%
  ungroup() %>%
  select(-abs_time, -relative_s) %>%
  mutate(
    date = as_date(timestamp),
    time = format(timestamp, "%H:%M:%S")
  )

# Require EVERY session to have real RTC timestamps before plotting on a
# datetime axis - a partial mix would render some panels on real time and
# others on nothing, which is easy to misread. `elapsed_s` (above) has no
# such restriction and is always populated, so it's the safe common axis
# to fall back on for the combined plots when sessions are mixed.
has_real_time <- !any(is.na(gaping$timestamp))

cat(sprintf(
  "%d session(s), %d total samples. Absolute timestamps for ALL sessions: %s\n",
  n_sessions, nrow(gaping),
  if (has_real_time) "YES (from onboard RTC)" else "NO (falling back to elapsed_s for combined plots)"
))
1 session(s), 910493 total samples. Absolute timestamps for ALL sessions: YES (from onboard RTC)

Some experiments mix RTC-timestamped and elapsed-only files (e.g. the RTC wasn’t detected on one boot). timestamp/date/time stay real per row for any session that has them, and elapsed_s is always populated for every session either way - only the combined multi-panel plots below fall back to elapsed_s for everyone when the mix isn’t uniform. Check which sessions are which here:

Code 
gaping %>%
  group_by(session_id, source_file) %>%
  summarise(has_real_time = !any(is.na(timestamp)), .groups = "drop") %>%
  kable(caption = "Timestamp availability per session")
Timestamp availability per session
session_id source_file has_real_time
1 gapelog_083.csv TRUE

DERIVED METRICS

Field magnitude

Field magnitude sqrt(x^2 + y^2 + z^2) is the primary gape proxy: it is rotation-invariant, so it doesn’t depend on exactly how the sensor is oriented relative to the magnet.

What is computed here is the raw magnitude, field_raw_uT - straight off the sensor, ambient field included. The quantity everything downstream actually uses, field_magnitude_uT, is that same magnitude with the ambient vector removed, and it cannot be computed yet: the background vector is a per-unit property that arrives with the calibration coefficients, so it is not known until the GAPE DISTANCE section has matched a curve to each session. See “Background correction” there.

Code 
gaping <- gaping %>%
  mutate(field_raw_uT = sqrt(x_uT^2 + y_uT^2 + z_uT^2))

gaping %>%
  select(x_uT, y_uT, z_uT, field_raw_uT) %>%
  summary()
      x_uT           y_uT              z_uT         field_raw_uT   
 Min.   :-686   Min.   :-1274.0   Min.   :-686.0   Min.   :   0.0  
 1st Qu.:-196   1st Qu.: -784.0   1st Qu.:-294.0   1st Qu.: 635.1  
 Median :-196   Median : -686.0   Median :-196.0   Median : 752.8  
 Mean   :-181   Mean   : -690.9   Mean   :-148.1   Mean   : 753.4  
 3rd Qu.: -98   3rd Qu.: -588.0   3rd Qu.: -98.0   3rd Qu.: 859.9  
 Max.   : 294   Max.   :   98.0   Max.   : 490.0   Max.   :1336.5

Sampling interval

Code 
gaping <- gaping %>%
  group_by(session_id) %>%
  arrange(elapsed_s, .by_group = TRUE) %>%
  mutate(sample_dt_s = elapsed_s - lag(elapsed_s)) %>%
  ungroup()

gaping %>%
  group_by(session_id, source_file) %>%
  summarise(
    n_samples     = n(),
    duration_min  = (max(elapsed_s) - min(elapsed_s)) / 60,
    median_dt_s   = median(sample_dt_s, na.rm = TRUE),
    max_gap_s     = max(sample_dt_s, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  kable(digits = 2, caption = "Per-session sampling summary")
Per-session sampling summary
session_id source_file n_samples duration_min median_dt_s max_gap_s
1 gapelog_083.csv 910493 15174.92 1 2

GAPE DISTANCE

Field magnitude answers “is the shell open?”. To answer “how far open?” the magnitude has to be pushed back through the calibration curve fit for that board, with that number of magnets:

\[ d = a B^{\,b} + c \]

Magnet count is part of the curve, not a detail: more magnets mean a stronger field at every distance, so the same reading maps to a different distance. Applying a 1-magnet curve to 3-magnet data reads a true 10 mm gape as roughly 6 mm, silently. The board UID is recorded on every data row, but the magnet count is not - it only exists in the deployment metadata below.

Everything in this section is skipped when apply_calibration = FALSE, and the field-magnitude results further down are unaffected either way.

Deployment metadata

One row per log file. n_magnets is the number that matters; oyster_id and notes are carried through for convenience but not used in any calculation.

This pilot is one oyster on one board, so there’s nothing to look up per file - n_magnets_used/pilot_oyster_id from the parameters chunk are just broadcast to every log file found above. (A multi-oyster deployment would replace this with read_csv(deployment_path, ...) keyed by source_file, same column names, so nothing downstream would need to change.)

Code 
deployment <- if (!apply_calibration) {
  NULL
} else {
  tibble(
    source_file = basename(log_files),
    n_magnets   = n_magnets_used,
    oyster_id   = pilot_oyster_id,
    notes       = NA_character_
  )
}

if (!is.null(deployment)) {
  kable(deployment, caption = "Deployment metadata as read")
} else {
  cat("`apply_calibration = FALSE` - working in field magnitude only.\n")
}
Deployment metadata as read
source_file n_magnets oyster_id notes
gapelog_083.csv 3 oyster_1 NA

Calibration coefficients

Code 
curves <- if (!apply_calibration) {
  NULL
} else {
  stopifnot("One or more `coef_path` files do not exist." = all(file.exists(coef_path)))

  all_coefs <- read_csv(coef_path, show_col_types = TRUE, col_types = cols(
    rp2040_uid = col_character(), n_magnets = col_integer(),
    model = col_character(), is_best_model = col_logical()
  ))

  # These curves were fit to background-corrected field, so the background
  # vector has to come with them. A coefficient file written before that was
  # exported would otherwise be applied to raw field and silently misread
  # every distance - the failure this stop exists to prevent.
  bg_cols <- c("bg_x_uT", "bg_y_uT", "bg_z_uT")
  missing_bg <- setdiff(bg_cols, names(all_coefs))
  if (length(missing_bg) > 0) {
    stop("Coefficient file(s) are missing the background column(s): ",
         paste(missing_bg, collapse = ", "),
         "\nThese were written by an older version of the calibration ",
         "notebooks. Re-knit `calibration_unit-*.Rmd` to regenerate them.")
  }

  if (calibration_model == "best") {
    all_coefs %>% filter(is_best_model)
  } else {
    all_coefs %>% filter(model == calibration_model)
  }
}

if (!is.null(curves)) {
  curves %>%
    select(unit, rp2040_uid, n_magnets, model, a, b, c, rmse_mm,
           fit_min_field_uT, fit_max_field_uT) %>%
    mutate(across(where(is.numeric), ~ signif(.x, 5))) %>%
    kable(caption = sprintf("Available calibration curves (model = \"%s\")",
                            calibration_model))
}
Available calibration curves (model = “best”)
unit rp2040_uid n_magnets model a b c rmse_mm fit_min_field_uT fit_max_field_uT
Unit-A df6470a31b8a2f2e 1 offset 201.52 -0.34691 -3.5245 0.34705 559.26 49745
Unit-A df6470a31b8a2f2e 2 offset 143.92 -0.25382 -7.7859 0.08991 643.39 65709
Unit-A df6470a31b8a2f2e 3 offset 154.93 -0.25477 -8.6720 0.11877 742.64 55898

Match one curve to each session

Three things have to line up: the session’s UID (from the data), its magnet count (from the metadata), and a curve fit for that combination. Each is checked separately so a failure says which one is missing rather than just producing NA distances.

Code 
session_curves <- NULL

if (apply_calibration) {
  # UID comes from the data itself - one board per log file.
  session_uid <- gaping %>%
    group_by(session_id, source_file) %>%
    summarise(rp2040_uid = paste(unique(rp2040_uid), collapse = "|"),
              .groups = "drop")

  mixed_uid <- session_uid %>% filter(str_detect(rp2040_uid, fixed("|")))
  if (nrow(mixed_uid) > 0) {
    stop("More than one rp2040_uid within a single log file:\n  ",
         paste(mixed_uid$source_file, mixed_uid$rp2040_uid, collapse = "\n  "))
  }

  # Magnet count comes from the metadata - it is nowhere in the data. With
  # `deployment` now built by broadcasting `n_magnets_used`/`pilot_oyster_id`
  # to every file in `log_files` (see `read-deployment-metadata`), this
  # should never actually trigger - it's left in as a guard in case that
  # assumption is ever changed.
  missing_meta <- setdiff(session_uid$source_file, deployment$source_file)
  if (length(missing_meta) > 0) {
    stop("No deployment metadata row (and so no magnet count) for:\n  ",
         paste(missing_meta, collapse = "\n  "),
         "\nCheck `n_magnets_used`/`pilot_oyster_id` in the parameters chunk.")
  }

  session_curves <- session_uid %>%
    left_join(deployment, by = "source_file") %>%
    left_join(
      curves %>% select(rp2040_uid, n_magnets, unit, model,
                        cal_a = a, cal_b = b, cal_c = c, rmse_mm,
                        bg_x_uT, bg_y_uT, bg_z_uT,
                        fit_min_field_uT, fit_max_field_uT),
      by = c("rp2040_uid", "n_magnets")
    )

  unmatched <- session_curves %>% filter(is.na(cal_a))
  if (nrow(unmatched) > 0) {
    stop(
      "No calibration curve for these UID / magnet-count combinations:\n  ",
      paste(sprintf("%s (uid %s, %d magnet(s))", unmatched$source_file,
                    unmatched$rp2040_uid, unmatched$n_magnets),
            collapse = "\n  "),
      "\nEither the wrong `coef_path` is loaded, or that magnet count was ",
      "never calibrated on that board."
    )
  }

  session_curves %>%
    select(session_id, source_file, unit, rp2040_uid, n_magnets, oyster_id,
           model, cal_a, cal_b, cal_c, rmse_mm) %>%
    mutate(across(where(is.numeric), ~ signif(.x, 5))) %>%
    kable(caption = "Calibration curve matched to each session")
}
Calibration curve matched to each session
session_id source_file unit rp2040_uid n_magnets oyster_id model cal_a cal_b cal_c rmse_mm
1 gapelog_083.csv Unit-A df6470a31b8a2f2e 3 oyster_1 offset 154.93 -0.25477 -8.672 0.11877

Background correction

The calibration curves were fit to field with the ambient vector removed, so the same has to be done here before a curve is applied - otherwise the curve is being handed a different quantity than the one it was built from, and it answers confidently in the wrong units.

Ambient field adds as a vector, so it is subtracted per axis and the magnitude recomputed from the corrected axes. Subtracting it from the magnitude would be wrong.

field_magnitude_uT is defined here, and it is what the rest of the notebook uses - event detection included. A fixed vector offset distorts magnitude unevenly as the magnet moves, so the corrected signal is the better one to detect gapes on, not just to convert to millimetres.

With apply_calibration = FALSE there is no background to subtract (it arrives with the coefficients), so field_magnitude_uT falls back to the raw magnitude and every field-magnitude result below still works - it is just uncorrected, and the printout says so.

Code 
if (apply_calibration) {
  gaping <- gaping %>%
    left_join(
      session_curves %>% select(session_id, n_magnets, cal_a, cal_b, cal_c,
                                bg_x_uT, bg_y_uT, bg_z_uT,
                                fit_min_field_uT, fit_max_field_uT),
      by = "session_id"
    ) %>%
    mutate(field_magnitude_uT = sqrt((x_uT - bg_x_uT)^2 +
                                     (y_uT - bg_y_uT)^2 +
                                     (z_uT - bg_z_uT)^2))

  gaping %>%
    group_by(session_id, source_file) %>%
    summarise(
      bg_uT          = sqrt(first(bg_x_uT)^2 + first(bg_y_uT)^2 + first(bg_z_uT)^2),
      mean_raw_uT    = mean(field_raw_uT),
      mean_corr_uT   = mean(field_magnitude_uT),
      shift_uT       = mean_corr_uT - mean_raw_uT,
      .groups = "drop"
    ) %>%
    mutate(across(where(is.numeric), ~ round(.x, 1))) %>%
    kable(caption = paste("Background removed per session (bg_uT is the",
                          "magnitude of the subtracted vector)"))
} else {
  gaping <- gaping %>% mutate(field_magnitude_uT = field_raw_uT)
  cat("`apply_calibration = FALSE`: no background available, so",
      "`field_magnitude_uT` is the RAW magnitude (ambient field included).\n")
}
Background removed per session (bg_uT is the magnitude of the subtracted vector)
session_id source_file bg_uT mean_raw_uT mean_corr_uT shift_uT
1 gapelog_083.csv 1231.6 753.4 1520.4 767

Convert field magnitude to distance

Code 
if (apply_calibration) {
  gaping <- gaping %>%
    mutate(gape_distance_mm = cal_a * field_magnitude_uT^cal_b + coalesce(cal_c, 0))
} else {
  gaping <- gaping %>% mutate(gape_distance_mm = NA_real_)
}

gaping %>%
  group_by(session_id, source_file) %>%
  summarise(
    min_distance_mm    = min(gape_distance_mm),
    median_distance_mm = median(gape_distance_mm),
    max_distance_mm    = max(gape_distance_mm),
    .groups = "drop"
  ) %>%
  kable(digits = 2, caption = "Gape distance range per session")
Gape distance range per session
session_id source_file min_distance_mm median_distance_mm max_distance_mm
1 gapelog_083.csv 13.03 15.33 18.58

Extrapolation check

A power law extrapolates confidently and wrongly. Readings outside the field range the curve was fit over produce a number that looks fine and isn’t, so they are counted here rather than left to be discovered later.

Negative distances are the one failure that announces itself: they mean the field was stronger than anything seen during calibration, which usually means the deployed magnet sits closer than the closest calibration hold, or that the wrong magnet count is recorded in the metadata.

Code 
if (apply_calibration) {
  extrapolation <- gaping %>%
    group_by(session_id, source_file) %>%
    summarise(
      n                 = n(),
      pct_below_range   = 100 * mean(field_magnitude_uT < fit_min_field_uT),
      pct_above_range   = 100 * mean(field_magnitude_uT > fit_max_field_uT),
      pct_negative_dist = 100 * mean(gape_distance_mm < 0),
      .groups = "drop"
    )

  extrapolation %>%
    kable(digits = 2,
          caption = "Samples falling outside the calibrated field range")

  worst <- max(extrapolation$pct_below_range + extrapolation$pct_above_range)
  if (worst > 5) {
    cat(sprintf(
      paste0("WARNING: up to %.1f%% of samples in a session fall outside the ",
             "calibrated field range.\nThose distances are extrapolated - ",
             "treat them as indicative only.\n"), worst))
  } else {
    cat(sprintf("At most %.1f%% of samples per session are extrapolated.\n", worst))
  }
}
At most 0.0% of samples per session are extrapolated.

Cross-check the magnet count

The extrapolation check above catches a magnet count that is uncalibrated or wildly wrong. It does not catch a plausible one: recording 1 magnet in the metadata when 2 were deployed puts every reading comfortably inside the 1-magnet curve’s fitted range, so nothing errors and nothing is flagged - the distances are just quietly compressed. In a test case, a true 7.0 mm gape read as 4.9 mm and a 1.0 mm closed baseline read as 0.14 mm.

There is no way to recover the magnet count from the data alone, so the check has to be a human one. Below, the session’s closed and open field levels are pushed through every curve calibrated for that board. Compare the closed_mm column against how far the magnet actually sat from the sensor when the shell was shut: if the selected row disagrees with the bench measurement and another row matches it, the metadata is wrong.

Code 
if (apply_calibration) {
  session_levels <- gaping %>%
    group_by(session_id, source_file, rp2040_uid) %>%
    summarise(
      closed_field_uT = quantile(field_magnitude_uT, baseline_quantile),
      open_field_uT   = quantile(field_magnitude_uT, 1 - baseline_quantile),
      .groups = "drop"
    )

  session_levels %>%
    left_join(session_curves %>% select(session_id, selected_n = n_magnets),
              by = "session_id") %>%
    inner_join(
      curves %>% select(rp2040_uid, n_magnets, alt_a = a, alt_b = b, alt_c = c),
      by = "rp2040_uid", relationship = "many-to-many"
    ) %>%
    mutate(
      selected  = if_else(n_magnets == selected_n, "<- selected", ""),
      closed_mm = alt_a * closed_field_uT^alt_b + coalesce(alt_c, 0),
      open_mm   = alt_a * open_field_uT^alt_b + coalesce(alt_c, 0),
      swing_mm  = open_mm - closed_mm
    ) %>%
    select(session_id, source_file, n_magnets, selected,
           closed_mm, open_mm, swing_mm) %>%
    arrange(session_id, n_magnets) %>%
    kable(digits = 2,
          caption = "What each calibrated magnet count would imply for this session")
}
What each calibrated magnet count would imply for this session
session_id source_file n_magnets selected closed_mm open_mm swing_mm
1 gapelog_083.csv 1 11.39 13.37 1.99
1 gapelog_083.csv 2 13.63 15.68 2.05
1 gapelog_083.csv 3 <- selected 14.22 16.42 2.20

GAP EVENT DETECTION

Threshold and open/closed state

Field magnitude is first smoothed with a rolling median (smooth_window_s) to suppress the sensor’s quantization noise; all open/closed calls are made on the smoothed signal (field_smooth_uT).

The "noise" method then needs two numbers per session:

  • Closed baseline (closed_mean_uT): the mean of smoothed samples at/above baseline_quantile. This only assumes the animal spends at least 1 - baseline_quantile of the recording fully closed.
  • Noise (noise_sd_uT): the noise left in the smoothed signal, estimated from differences between smoothed values one smoothing window apart (consecutive smoothed values share most of their window, so their differences would understate it). Differencing cancels the slow gape signal and leaves ~sqrt(2) x the noise; using the MAD keeps real opening/closing steps from inflating it. This replaces the SD of the top-quantile samples, which describes only a truncated slice of the distribution and doesn’t measure noise.

The threshold sits n_sd_threshold noise SDs below the closed baseline, so a gape event is a drop that’s statistically distinct from noise - not just a fixed percentage below baseline. The noise is measured on the same smoothed signal the threshold is applied to.

Code 
# Smoothing window in samples at the session's typical sampling interval,
# forced odd as runmed requires; 1 = no smoothing.
smooth_k <- function(elapsed_s, window_s) {
  k <- round(window_s / median(diff(elapsed_s)))
  if (is.na(k) || k < 3) return(1)
  if (k %% 2 == 0) k <- k + 1
  k
}

smooth_magnitude <- function(x, k) {
  if (k < 3) return(x)
  as.numeric(runmed(x, k, endrule = "median"))
}

# Noise SD from differences `lag` samples apart.
noise_sd <- function(x, lag) {
  d <- diff(x, lag = lag)
  s <- mad(d, na.rm = TRUE) / sqrt(2)
  # A heavily quantized signal can give MAD = 0; fall back to the SD then.
  if (s == 0) s <- sd(d, na.rm = TRUE) / sqrt(2)
  s
}

gaping <- gaping %>%
  arrange(session_id, elapsed_s) %>%
  group_by(session_id) %>%
  mutate(
    k               = smooth_k(elapsed_s, smooth_window_s),
    field_smooth_uT = smooth_magnitude(field_magnitude_uT, first(k)),
    baseline_uT    = quantile(field_smooth_uT, baseline_quantile, na.rm = TRUE),
    closed_mean_uT = mean(field_smooth_uT[field_smooth_uT >= baseline_uT], na.rm = TRUE),
    noise_sd_uT    = noise_sd(field_smooth_uT, lag = first(k)),
    threshold_uT = if (threshold_method == "noise") {
      closed_mean_uT - n_sd_threshold * noise_sd_uT
    } else if (threshold_method == "relative") {
      baseline_uT * (1 - gape_drop_frac)
    } else if (threshold_method == "absolute") {
      gape_threshold_uT
    } else {
      NA_real_
    },
    is_open = field_smooth_uT < threshold_uT
  ) %>%
  ungroup() %>%
  select(-k)

stopifnot("`threshold_method` must be one of \"noise\", \"relative\", \"absolute\"." =
            !any(is.na(gaping$threshold_uT)))

gaping %>%
  group_by(session_id, source_file) %>%
  summarise(
    smooth_window_s  = smooth_window_s,
    closed_mean_uT   = first(closed_mean_uT),
    noise_sd_uT      = first(noise_sd_uT),
    threshold_uT     = first(threshold_uT),
    pct_samples_open = 100 * mean(is_open),
    .groups = "drop"
  ) %>%
  kable(digits = 2, caption = "Detection threshold per session")
Detection threshold per session
session_id source_file smooth_window_s closed_mean_uT noise_sd_uT threshold_uT pct_samples_open
1 gapelog_083.csv 15 1837 24.94 1762.19 85.81

Percent opening

A continuous measure of how far open the animal is, independent of any open/closed cutoff: each session is rescaled so 0% = closed baseline (closed_mean_uT) and 100% = widest opening seen (the open_ref_quantile of smoothed magnitude). It’s computed on gape distance when calibration is applied, since distance is linear in valve gape; without calibration it falls back to field magnitude, which is only a rough proxy because field falls off non-linearly with distance. Values are clamped to 0-100%.

Code 
# Puts field magnitude on the scale percent opening is measured on, oriented
# so larger = more open. The calibration coefficients are only evaluated when
# calibration is applied.
opening_scale <- function(field_uT, a, b, c) {
  if (apply_calibration) a * field_uT^b + coalesce(c, 0) else -field_uT
}

gaping <- gaping %>%
  group_by(session_id) %>%
  mutate(
    open_ref_uT = quantile(field_smooth_uT, open_ref_quantile, na.rm = TRUE),
    s_closed    = opening_scale(closed_mean_uT, cal_a, cal_b, cal_c),
    s_open      = opening_scale(open_ref_uT, cal_a, cal_b, cal_c),
    pct_opening = 100 * (opening_scale(field_smooth_uT, cal_a, cal_b, cal_c) - s_closed) /
                    (s_open - s_closed),
    pct_opening = pmin(pmax(pct_opening, 0), 100),
    # Where the open/closed threshold falls on the same scale.
    threshold_pct_opening = 100 * (opening_scale(threshold_uT, cal_a, cal_b, cal_c) - s_closed) /
                              (s_open - s_closed)
  ) %>%
  select(-s_closed, -s_open) %>%
  ungroup()

gaping %>%
  group_by(session_id, source_file) %>%
  summarise(
    open_ref_uT           = first(open_ref_uT),
    threshold_pct_opening = first(threshold_pct_opening),
    mean_pct_opening      = mean(pct_opening),
    median_pct_opening    = median(pct_opening),
    .groups = "drop"
  ) %>%
  kable(digits = 1, caption = "Percent opening per session")
Percent opening per session
session_id source_file open_ref_uT threshold_pct_opening mean_pct_opening median_pct_opening
1 gapelog_083.csv 1264.7 10.7 51.9 45.2

Build the event table

Contiguous runs of is_open become candidate events; runs shorter than min_event_duration_s are dropped, and events separated by less than min_gap_between_s are merged into one.

Code 
# 1. Contiguous runs of the same open/closed state, per session.
runs <- gaping %>%
  group_by(session_id) %>%
  mutate(run_id = cumsum(is_open != lag(is_open, default = first(is_open)))) %>%
  group_by(session_id, source_file, run_id, is_open) %>%
  summarise(
    start_s   = min(elapsed_s),
    end_s     = max(elapsed_s),
    start_time = min(timestamp),
    end_time   = max(timestamp),
    n_samples = n(),
    mean_magnitude_uT = mean(field_magnitude_uT),
    min_magnitude_uT  = min(field_magnitude_uT),
    .groups   = "drop"
  ) %>%
  mutate(duration_s = end_s - start_s)

# 2. Keep open runs long enough to count as real events.
events <- runs %>%
  filter(is_open, duration_s >= min_event_duration_s) %>%
  arrange(session_id, start_s)

# 3. Merge events separated by a shorter closure than min_gap_between_s.
events <- events %>%
  group_by(session_id) %>%
  mutate(
    gap_before_s = start_s - lag(end_s),
    new_event    = is.na(gap_before_s) | gap_before_s >= min_gap_between_s,
    event_group  = cumsum(new_event)
  ) %>%
  group_by(session_id, source_file, event_group) %>%
  summarise(
    start_s           = min(start_s),
    end_s             = max(end_s),
    start_time        = min(start_time),
    end_time          = max(end_time),
    # Weighted mean must be computed BEFORE n_samples is collapsed to a
    # single sum, otherwise the weights vector no longer matches the values.
    mean_magnitude_uT = weighted.mean(mean_magnitude_uT, n_samples),
    min_magnitude_uT  = min(min_magnitude_uT),
    n_samples         = sum(n_samples),
    .groups = "drop"
  ) %>%
  group_by(session_id) %>%
  mutate(
    event_id    = row_number(),
    duration_s  = end_s - start_s,
    duration_min = duration_s / 60
  ) %>%
  ungroup()

cat(sprintf("Detected %d gape event(s) across %d session(s).\n",
            nrow(events), n_sessions))
Detected 82 gape event(s) across 1 session(s).
Code 
events %>%
  select(session_id, event_id, start_s, end_s, duration_s, duration_min,
         mean_magnitude_uT, min_magnitude_uT, n_samples) %>%
  kable(digits = 2, caption = "Detected gape events")
Detected gape events
session_id event_id start_s end_s duration_s duration_min mean_magnitude_uT min_magnitude_uT n_samples
1 1 0 13068 13068 217.80 1484.57 965.40 13069
1 2 13151 26502 13351 222.52 1588.96 1198.25 13312
1 3 26598 38286 11688 194.80 1587.87 1198.25 11682
1 4 38435 49225 10790 179.83 1588.01 1142.29 10784
1 5 62662 74159 11497 191.62 1592.52 1217.48 11498
1 6 74293 80491 6198 103.30 1540.78 1146.76 6199
1 7 80527 84545 4018 66.97 1601.37 1264.73 4019
1 8 84742 84773 31 0.52 1725.73 1541.65 32
1 9 84809 99795 14986 249.77 1586.58 1174.50 14983
1 10 99910 112530 12620 210.33 1585.25 1205.33 12617
1 11 112670 122671 10001 166.68 1567.36 1198.25 10002
1 12 122816 135224 12408 206.80 1582.80 1205.33 12405
1 13 135411 148574 13163 219.38 1587.65 1126.32 13161
1 14 148819 148893 74 1.23 1702.74 1541.65 75
1 15 148926 164615 15689 261.48 1596.21 1202.91 15680
1 16 164828 177128 12300 205.00 1618.16 1261.68 12286
1 17 177313 178410 1097 18.28 1678.35 1419.63 1098
1 18 178441 188256 9815 163.58 1612.04 1201.07 9796
1 19 188429 198144 9715 161.92 1613.73 1201.07 9673
1 20 219604 236074 16470 274.50 1577.36 1146.76 16471
1 21 236242 243208 6966 116.10 1562.25 1205.33 6967
1 22 243347 255912 12565 209.42 1576.14 1198.25 12561
1 23 256121 264929 8808 146.80 1572.45 1193.97 8809
1 24 265165 283511 18346 305.77 1574.46 1193.97 18346
1 25 283668 284383 715 11.92 1670.06 1353.66 694
1 26 284427 301823 17396 289.93 1553.58 1174.50 17397
1 27 301988 315801 13813 230.22 1568.84 1126.32 13806
1 28 316039 327281 11242 187.37 1552.82 1146.76 11243
1 29 327444 338640 11196 186.60 1579.60 1146.76 11175
1 30 338800 343487 4687 78.12 1601.93 1264.73 4688
1 31 343670 351381 7711 128.52 1595.87 1204.80 7711
1 32 351651 361368 9717 161.95 1598.45 1198.25 9706
1 33 376133 388665 12532 208.87 1553.21 1180.31 12533
1 34 388801 406569 17768 296.13 1542.66 1126.32 17769
1 35 406702 423439 16737 278.95 1565.44 1146.76 16732
1 36 423709 431605 7896 131.60 1576.83 1205.33 7896
1 37 431802 441737 9935 165.58 1606.40 1142.29 9928
1 38 454342 469407 15065 251.08 1543.75 1146.76 15065
1 39 469491 478686 9195 153.25 1525.37 1197.73 9185
1 40 491698 510021 18323 305.38 1449.04 1059.63 18323
1 41 524847 546053 21206 353.43 1406.94 987.97 21207
1 42 555133 559454 4321 72.02 1365.30 1055.39 4322
1 43 559601 571036 11435 190.58 1341.57 987.97 11435
1 44 578028 590586 12558 209.30 1341.21 987.97 12558
1 45 594668 605771 11103 185.05 1309.33 921.68 11103
1 46 608857 621593 12736 212.27 1315.87 916.80 12737
1 47 621689 638286 16597 276.62 1345.85 916.80 16598
1 48 642654 682880 40226 670.43 1304.03 926.03 40223
1 49 683014 684473 1459 24.32 1488.48 1126.32 1460
1 50 685126 710552 25426 423.77 1349.81 916.80 25427
1 51 711554 711584 30 0.50 1760.82 1609.75 31
1 52 715934 771724 55790 929.83 1338.76 931.53 55791
1 53 771763 792768 21005 350.08 1341.13 987.97 20991
1 54 792869 792904 35 0.58 1754.59 1566.07 36
1 55 792989 807373 14384 239.73 1346.08 987.97 14314
1 56 807443 808035 592 9.87 1720.99 1446.88 535
1 57 808128 808486 358 5.97 1726.79 1476.92 297
1 58 808533 808566 33 0.55 1711.91 1495.57 35
1 59 808606 808653 47 0.78 1715.54 1541.65 48
1 60 808690 808804 114 1.90 1731.19 1476.92 114
1 61 808861 808949 88 1.47 1738.09 1490.29 89
1 62 809015 809218 203 3.38 1728.84 1498.98 176
1 63 809271 809386 115 1.92 1731.77 1544.97 109
1 64 809437 809514 77 1.28 1730.57 1476.92 78
1 65 809560 809759 199 3.32 1735.90 1554.47 178
1 66 809818 809974 156 2.60 1733.57 1484.58 156
1 67 810092 838480 28388 473.13 1336.68 930.86 28387
1 68 838518 856741 18223 303.72 1324.74 916.80 18178
1 69 856806 871207 14401 240.02 1327.71 922.36 14350
1 70 871267 891551 20284 338.07 1315.78 987.97 20283
1 71 891584 891622 38 0.63 1739.12 1609.75 39
1 72 891659 891696 37 0.62 1733.66 1556.64 38
1 73 891741 891784 43 0.72 1753.94 1618.18 43
1 74 891951 891982 31 0.52 1745.48 1544.97 31
1 75 892271 892316 45 0.75 1727.25 1556.64 46
1 76 892511 892551 40 0.67 1754.45 1545.68 42
1 77 893490 893546 56 0.93 1762.46 1615.02 57
1 78 894736 894781 45 0.75 1767.02 1544.97 46
1 79 895092 895129 37 0.62 1765.06 1618.18 38
1 80 896877 896907 30 0.50 1771.45 1618.18 31
1 81 897187 897233 46 0.77 1758.51 1612.92 47
1 82 897556 910495 12939 215.65 1323.90 987.97 12940

Gape distance per event

Distances are attached to the existing event table rather than folded into its construction, so event detection stays entirely field-magnitude-based and identical whether or not calibration is applied.

max_gape_mm is the furthest the valve opened during the event (the minimum field magnitude - field falls as the magnet moves away). gape_amplitude_mm measures that against the session’s closed baseline, which is usually the more meaningful number: it is how far the animal opened relative to shut, rather than the absolute magnet-to-sensor separation.

Code 
if (apply_calibration && nrow(events) > 0) {
  # Distance at the closed-shell baseline, one value per session.
  closed_distance <- gaping %>%
    group_by(session_id) %>%
    summarise(
      closed_distance_mm = first(cal_a) * first(closed_mean_uT)^first(cal_b) +
                           coalesce(first(cal_c), 0),
      .groups = "drop"
    )

  # Samples inside each event window. A range (non-equi) join only ever
  # materializes the samples that actually fall inside an event - an equi
  # join on session_id followed by a filter would fan out to
  # (samples x events) first, which for ~900k samples and ~1,300 events is
  # over a billion rows and crashes R.
  event_samples <- gaping %>%
    select(session_id, elapsed_s, gape_distance_mm) %>%
    inner_join(events %>% select(session_id, event_id, start_s, end_s),
               by = join_by(session_id, between(elapsed_s, start_s, end_s)))

  event_distance <- event_samples %>%
    group_by(session_id, event_id) %>%
    summarise(
      max_gape_mm  = max(gape_distance_mm),
      mean_gape_mm = mean(gape_distance_mm),
      .groups = "drop"
    ) %>%
    left_join(closed_distance, by = "session_id") %>%
    mutate(gape_amplitude_mm = max_gape_mm - closed_distance_mm)

  events <- events %>%
    left_join(event_distance, by = c("session_id", "event_id"))

  events %>%
    select(session_id, event_id, duration_min, closed_distance_mm,
           max_gape_mm, mean_gape_mm, gape_amplitude_mm) %>%
    kable(digits = 2, caption = "Gape distance per detected event (mm)")
} else if (nrow(events) == 0) {
  cat("No events detected - no per-event distances to compute.\n")
} else {
  cat("`apply_calibration = FALSE` - no per-event distances computed.\n")
}
Gape distance per detected event (mm)
session_id event_id duration_min closed_distance_mm max_gape_mm mean_gape_mm gape_amplitude_mm
1 1 217.80 14.16 18.23 15.46 4.07
1 2 222.52 14.16 16.79 15.03 2.63
1 3 194.80 14.16 16.79 15.04 2.63
1 4 179.83 14.16 17.10 15.04 2.94
1 5 191.62 14.16 16.68 15.02 2.52
1 6 103.30 14.16 17.07 15.22 2.91
1 7 66.97 14.16 16.44 14.99 2.28
1 8 0.52 14.16 15.20 14.54 1.04
1 9 249.77 14.16 16.92 15.05 2.76
1 10 210.33 14.16 16.75 15.05 2.59
1 11 166.68 14.16 16.79 15.12 2.63
1 12 206.80 14.16 16.75 15.06 2.59
1 13 219.38 14.16 17.19 15.04 3.03
1 14 1.23 14.16 15.20 14.62 1.04
1 15 261.48 14.16 16.76 15.01 2.60
1 16 205.00 14.16 16.45 14.92 2.29
1 17 18.28 14.16 15.71 14.70 1.55
1 18 163.58 14.16 16.77 14.95 2.61
1 19 161.92 14.16 16.77 14.94 2.61
1 20 274.50 14.16 17.07 15.08 2.91
1 21 116.10 14.16 16.75 15.14 2.59
1 22 209.42 14.16 16.79 15.09 2.63
1 23 146.80 14.16 16.81 15.10 2.65
1 24 305.77 14.16 16.81 15.09 2.65
1 25 11.92 14.16 16.01 14.72 1.85
1 26 289.93 14.16 16.92 15.17 2.76
1 27 230.22 14.16 17.19 15.11 3.03
1 28 187.37 14.16 17.07 15.18 2.91
1 29 186.60 14.16 17.07 15.07 2.91
1 30 78.12 14.16 16.44 14.99 2.28
1 31 128.52 14.16 16.75 15.01 2.59
1 32 161.95 14.16 16.79 15.00 2.63
1 33 208.87 14.16 16.88 15.17 2.72
1 34 296.13 14.16 17.19 15.22 3.03
1 35 278.95 14.16 17.07 15.13 2.91
1 36 131.60 14.16 16.75 15.08 2.59
1 37 165.58 14.16 17.10 14.97 2.94
1 38 251.08 14.16 17.07 15.21 2.91
1 39 153.25 14.16 16.79 15.29 2.63
1 40 305.38 14.16 17.60 15.61 3.44
1 41 353.43 14.16 18.07 15.79 3.91
1 42 72.02 14.16 17.62 15.98 3.46
1 43 190.58 14.16 18.07 16.09 3.91
1 44 209.30 14.16 18.07 16.09 3.91
1 45 185.05 14.16 18.55 16.24 4.39
1 46 212.27 14.16 18.58 16.21 4.42
1 47 276.62 14.16 18.58 16.08 4.42
1 48 670.43 14.16 18.51 16.27 4.35
1 49 24.32 14.16 17.19 15.45 3.03
1 50 423.77 14.16 18.58 16.05 4.42
1 51 0.50 14.16 14.94 14.42 0.78
1 52 929.83 14.16 18.47 16.10 4.31
1 53 350.08 14.16 18.07 16.09 3.91
1 54 0.58 14.16 15.11 14.44 0.95
1 55 239.73 14.16 18.07 16.06 3.91
1 56 9.87 14.16 15.59 14.55 1.43
1 57 5.97 14.16 15.47 14.51 1.31
1 58 0.55 14.16 15.39 14.59 1.23
1 59 0.78 14.16 15.20 14.57 1.04
1 60 1.90 14.16 15.47 14.52 1.31
1 61 1.47 14.16 15.41 14.50 1.25
1 62 3.38 14.16 15.37 14.51 1.21
1 63 1.92 14.16 15.19 14.50 1.03
1 64 1.28 14.16 15.47 14.52 1.31
1 65 3.32 14.16 15.15 14.48 0.99
1 66 2.60 14.16 15.43 14.51 1.27
1 67 473.13 14.16 18.48 16.11 4.32
1 68 303.72 14.16 18.58 16.17 4.42
1 69 240.02 14.16 18.54 16.16 4.38
1 70 338.07 14.16 18.07 16.21 3.91
1 71 0.63 14.16 14.94 14.49 0.78
1 72 0.62 14.16 15.14 14.51 0.98
1 73 0.72 14.16 14.92 14.45 0.76
1 74 0.52 14.16 15.19 14.47 1.03
1 75 0.75 14.16 15.14 14.53 0.98
1 76 0.67 14.16 15.19 14.44 1.03
1 77 0.93 14.16 14.92 14.41 0.76
1 78 0.75 14.16 15.19 14.40 1.03
1 79 0.62 14.16 14.91 14.40 0.75
1 80 0.50 14.16 14.91 14.38 0.75
1 81 0.77 14.16 14.93 14.42 0.77
1 82 215.65 14.16 18.07 16.17 3.91

Gape distance summary

Code 
if (apply_calibration && nrow(events) > 0) {
  events %>%
    group_by(session_id, source_file) %>%
    summarise(
      n_events             = n(),
      closed_distance_mm   = first(closed_distance_mm),
      mean_max_gape_mm     = mean(max_gape_mm),
      largest_gape_mm      = max(max_gape_mm),
      mean_amplitude_mm    = mean(gape_amplitude_mm),
      largest_amplitude_mm = max(gape_amplitude_mm),
      .groups = "drop"
    ) %>%
    kable(digits = 2, caption = "Gape distance summary by session (mm)")
}
Gape distance summary by session (mm)
session_id source_file n_events closed_distance_mm mean_max_gape_mm largest_gape_mm mean_amplitude_mm largest_amplitude_mm
1 gapelog_083.csv 82 14.16 16.63 18.58 2.47 4.42

Event frequency

Code 
event_frequency <- events %>%
  group_by(session_id, source_file) %>%
  summarise(n_events = n(), .groups = "drop") %>%
  left_join(
    gaping %>%
      group_by(session_id) %>%
      summarise(
        recording_hr = (max(elapsed_s) - min(elapsed_s)) / 3600,
        .groups = "drop"
      ),
    by = "session_id"
  ) %>%
  mutate(
    events_per_hr  = n_events / recording_hr,
    events_per_day = events_per_hr * 24
  )

event_frequency %>%
  kable(digits = 3, caption = "Gape event frequency per session")
Gape event frequency per session
session_id source_file n_events recording_hr events_per_hr events_per_day
1 gapelog_083.csv 82 252.915 0.324 7.781

Event durations

Code 
duration_stats <- events %>%
  group_by(session_id, source_file) %>%
  summarise(
    n_events        = n(),
    mean_duration_s = mean(duration_s),
    sd_duration_s   = sd(duration_s),
    min_duration_s  = min(duration_s),
    median_duration_s = median(duration_s),
    max_duration_s  = max(duration_s),
    total_open_s    = sum(duration_s),
    .groups = "drop"
  ) %>%
  left_join(
    gaping %>%
      group_by(session_id) %>%
      summarise(
        recording_s = max(elapsed_s) - min(elapsed_s), .groups = "drop"
      ),
    by = "session_id"
  ) %>%
  mutate(pct_time_open = 100 * total_open_s / recording_s)

duration_stats %>%
  kable(digits = 2, caption = "Gape event duration statistics (seconds)")
Gape event duration statistics (seconds)
session_id source_file n_events mean_duration_s sd_duration_s min_duration_s median_duration_s max_duration_s total_open_s recording_s pct_time_open
1 gapelog_083.csv 82 9445.95 9633.54 30 9875 55790 774568 910495 85.07

Time between events

Inter-event interval = time from the end of one event to the start of the next (i.e. how long the animal stayed closed between gapes).

Code 
intervals <- events %>%
  group_by(session_id, source_file) %>%
  arrange(start_s, .by_group = TRUE) %>%
  mutate(interval_s = start_s - lag(end_s)) %>%
  filter(!is.na(interval_s)) %>%
  ungroup()

interval_stats <- intervals %>%
  group_by(session_id, source_file) %>%
  summarise(
    n_intervals      = n(),
    mean_interval_s  = mean(interval_s),
    sd_interval_s    = sd(interval_s),
    shortest_interval_s = min(interval_s),
    median_interval_s   = median(interval_s),
    longest_interval_s  = max(interval_s),
    .groups = "drop"
  )

if (nrow(interval_stats) > 0) {
  interval_stats %>%
    kable(digits = 2,
          caption = "Time between consecutive gape events (seconds)")
} else {
  cat("Fewer than 2 events per session - no inter-event intervals to report.\n")
}
Time between consecutive gape events (seconds)
session_id source_file n_intervals mean_interval_s sd_interval_s shortest_interval_s median_interval_s longest_interval_s
1 gapelog_083.csv 81 1678.11 4157.51 31 157 21460

Combined summary

Code 
gape_summary <- event_frequency %>%
  select(session_id, source_file, n_events, recording_hr,
         events_per_hr, events_per_day) %>%
  left_join(
    duration_stats %>%
      select(session_id, mean_duration_s, min_duration_s, max_duration_s,
             total_open_s, pct_time_open),
    by = "session_id"
  ) %>%
  left_join(
    interval_stats %>%
      select(session_id, mean_interval_s, shortest_interval_s,
             longest_interval_s),
    by = "session_id"
  )

gape_summary %>%
  kable(digits = 2, caption = "Gape event summary, all metrics by session")
Gape event summary, all metrics by session
session_id source_file n_events recording_hr events_per_hr events_per_day mean_duration_s min_duration_s max_duration_s total_open_s pct_time_open mean_interval_s shortest_interval_s longest_interval_s
1 gapelog_083.csv 82 252.92 0.32 7.78 9445.95 30 55790 774568 85.07 1678.11 31 21460

PLOTS

Field magnitude with detected events

Shaded bands are detected gape events; the dashed line is the detection threshold.

Code 
x_var  <- if (has_real_time) "timestamp" else "elapsed_s"
x_lab  <- if (has_real_time) "Time" else "Elapsed time (s)"

p_magnitude <- ggplot(gaping, aes(x = .data[[x_var]], y = field_magnitude_uT))

if (nrow(events) > 0) {
  event_bands <- events %>%
    mutate(
      xmin = if (has_real_time) start_time else start_s,
      xmax = if (has_real_time) end_time   else end_s
    )
  p_magnitude <- p_magnitude +
    geom_rect(
      data = event_bands,
      aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
      inherit.aes = FALSE, fill = "steelblue", alpha = 0.20
    )
}

p_magnitude <- p_magnitude +
  geom_line(linewidth = 0.3, colour = "grey60") +
  geom_line(aes(y = field_smooth_uT), linewidth = 0.3) +
  # One threshold row per session - inheriting `gaping` here would draw the
  # same line once per sample.
  geom_hline(
    data = gaping %>% distinct(source_file, threshold_uT),
    aes(yintercept = threshold_uT), linetype = "dashed", color = "firebrick"
  ) +
  facet_wrap(~source_file, scales = "free_x", ncol = 1) +
  labs(
    title = "Field magnitude over time, with detected gape events shaded",
    subtitle = paste0("Background-corrected; magnitude falls as the shell opens ",
                      "(magnet moves away from sensor).\nGrey = raw, black = ",
                      smooth_window_s, " s rolling median, red dashed = open/closed threshold"),
    x = x_lab, y = "Field magnitude, background-corrected (uT)"
  ) +
  theme_minimal()

p_magnitude

Gape distance over time

The same trace as above, in millimetres instead of microtesla - and the other way up, since distance rises as the shell opens. The dashed line is the closed-shell baseline distance.

Code 
if (apply_calibration) {
  p_distance <- ggplot(gaping, aes(x = .data[[x_var]], y = gape_distance_mm))

  if (nrow(events) > 0) {
    p_distance <- p_distance +
      geom_rect(
        data = event_bands,
        aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
        inherit.aes = FALSE, fill = "steelblue", alpha = 0.20
      )
  }

  p_distance <- p_distance +
    geom_line(linewidth = 0.3) +
    geom_hline(
      data = gaping %>%
        group_by(session_id, source_file) %>%
        summarise(closed_mm = first(cal_a) * first(closed_mean_uT)^first(cal_b) +
                              coalesce(first(cal_c), 0), .groups = "drop"),
      aes(yintercept = closed_mm), linetype = "dashed", colour = "firebrick"
    ) +
    facet_wrap(~source_file, scales = "free_x", ncol = 1) +
    labs(
      title = "Gape distance over time, with detected gape events shaded",
      subtitle = "Calibrated magnet-to-sensor distance; dashed line is the closed baseline",
      x = x_lab, y = "Gape distance (mm)"
    ) +
    theme_minimal()

  p_distance
} else {
  cat("`apply_calibration = FALSE` - no distance plot.\n")
}

How far each gape opened

Code 
if (apply_calibration && nrow(events) > 0) {
  events %>%
    mutate(x_pos = if (has_real_time) start_time else start_s) %>%
    ggplot(aes(x = x_pos, y = gape_amplitude_mm)) +
    geom_segment(aes(xend = x_pos, yend = 0), colour = "steelblue", linewidth = 0.4) +
    geom_point(size = 2, colour = "steelblue") +
    facet_wrap(~source_file, scales = "free_x", ncol = 1) +
    labs(
      title = "Gape amplitude per event",
      subtitle = "Maximum opening of each event, relative to the closed-shell baseline",
      x = x_lab, y = "Gape amplitude (mm)"
    ) +
    theme_minimal()
} else {
  cat("Nothing to plot (calibration off, or no events detected).\n")
}

Per-axis time series

Pivoting to long form here is purely a plotting convenience - it is not how the data is stored.

Code 
gaping_long <- gaping %>%
  select(session_id, source_file, all_of(x_var), x_uT, y_uT, z_uT) %>%
  pivot_longer(
    cols = c(x_uT, y_uT, z_uT), names_to = "axis", values_to = "field_uT"
  )

p_axes <- ggplot(gaping_long, aes(x = .data[[x_var]], y = field_uT, color = axis)) +
  geom_line(linewidth = 0.3) +
  facet_grid(axis ~ source_file, scales = "free") +
  labs(
    title = "Magnetic field components over time",
    x = x_lab, y = "Field (uT)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

p_axes

Event duration distribution

Durations span several orders of magnitude (tens of seconds to hours), so the x-axis is log10-scaled. Dashed line = median.

Code 
# Human-readable tick marks for log-scaled time axes (values in seconds)
time_breaks_s <- c(30, 60, 300, 600, 1800, 3600, 3 * 3600, 10 * 3600, 24 * 3600)
time_labels   <- c("30 s", "1 min", "5 min", "10 min", "30 min",
                   "1 h", "3 h", "10 h", "24 h")

# Format seconds as s / min / h for on-plot labels
fmt_duration <- function(s) {
  case_when(
    s < 60   ~ sprintf("%.0f s", s),
    s < 3600 ~ sprintf("%.1f min", s / 60),
    TRUE     ~ sprintf("%.1f h", s / 3600)
  )
}

duration_medians <- events %>%
  group_by(source_file) %>%
  summarise(median_s = median(duration_s), .groups = "drop")

if (nrow(events) > 0) {
  p_durations <- ggplot(events, aes(x = duration_s)) +
    geom_histogram(bins = 40, fill = "steelblue", color = "white") +
    geom_vline(
      data = duration_medians,
      aes(xintercept = median_s), linetype = "dashed"
    ) +
    geom_text(
      data = duration_medians,
      aes(x = median_s, y = Inf, label = paste("median =", fmt_duration(median_s))),
      hjust = -0.1, vjust = 1.5, size = 3.5
    ) +
    scale_x_log10(breaks = time_breaks_s, labels = time_labels) +
    facet_wrap(~source_file, ncol = 1) +
    labs(
      title = "Distribution of gape event durations",
      subtitle = "Dashed line = median event duration",
      x = "Event duration (log scale)", y = "Count of events"
    ) +
    theme_minimal()
  p_durations
} else {
  cat("No events detected - nothing to plot.\n")
}

Event duration over time

Shows whether gape bouts get longer or shorter across the deployment.

Code 
if (nrow(events) > 0) {
  p_dur_time <- events %>%
    mutate(x_pos = if (has_real_time) start_time else start_s) %>%
    ggplot(aes(x = x_pos, y = duration_min)) +
    geom_point(size = 2, color = "steelblue") +
    geom_line(linewidth = 0.3, alpha = 0.5, color = "steelblue") +
    facet_wrap(~source_file, scales = "free_x", ncol = 1) +
    labs(
      title = "Gape event duration over the course of the deployment",
      x = x_lab, y = "Event duration (minutes)"
    ) +
    theme_minimal()
  p_dur_time
} else {
  cat("No events detected - nothing to plot.\n")
}

Inter-event intervals

Log10-scaled x-axis for the same reason as the duration plot. Dashed line = median.

Code 
if (nrow(intervals) > 0) {
  interval_medians <- intervals %>%
    group_by(source_file) %>%
    summarise(median_s = median(interval_s), .groups = "drop")

  p_intervals <- ggplot(intervals, aes(x = interval_s)) +
    geom_histogram(bins = 40, fill = "darkorange", color = "white") +
    geom_vline(
      data = interval_medians,
      aes(xintercept = median_s), linetype = "dashed"
    ) +
    geom_text(
      data = interval_medians,
      aes(x = median_s, y = Inf, label = paste("median =", fmt_duration(median_s))),
      hjust = -0.1, vjust = 1.5, size = 3.5
    ) +
    scale_x_log10(breaks = time_breaks_s, labels = time_labels) +
    facet_wrap(~source_file, ncol = 1) +
    labs(
      title = "Time between consecutive gape events",
      subtitle = "Dashed line = median interval between events",
      x = "Interval between events (log scale)", y = "Count"
    ) +
    theme_minimal()
  p_intervals
} else {
  cat("Fewer than 2 events per session - no intervals to plot.\n")
}

Events per hour over time

Number of gape events starting in each hour of the deployment. Shows whether gaping activity changes over time (e.g. day/night rhythm, settling in after deployment). Partial hours at the start and end of a session will read low.

Code 
if (nrow(events) > 0) {
  p_rate_time <- events %>%
    mutate(x_pos = if (has_real_time) start_time else start_s / 3600) %>%
    ggplot(aes(x = x_pos)) +
    geom_histogram(
      binwidth = if (has_real_time) 3600 else 1,  # 1 h bins (POSIXct is in s)
      boundary = 0, fill = "steelblue", color = "white"
    ) +
    facet_wrap(~source_file, scales = "free_x", ncol = 1) +
    labs(
      title = "Gape events per hour across the deployment",
      x = if (has_real_time) "Time" else "Elapsed time (h)",
      y = "Events per hour"
    ) +
    theme_minimal()
  p_rate_time
} else {
  cat("No events detected - nothing to plot.\n")
}

Events per hour by session

Overall event rate for each log file. Only drawn when there is more than one session to compare.

Code 
if (n_sessions > 1) {
  p_rate <- ggplot(event_frequency, aes(x = source_file, y = events_per_hr)) +
    geom_col(fill = "steelblue") +
    labs(
      title = "Gape event frequency by session",
      x = "Log file (session)", y = "Events per hour"
    ) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  p_rate
} else {
  cat("Single session - see the events-per-hour table above for its overall rate.\n")
}
Single session - see the events-per-hour table above for its overall rate.

Percent opening over time

One-minute means of percent opening (0% = closed baseline, 100% = widest opening seen).

Code 
pct_by_min <- gaping %>%
  mutate(minute = floor(elapsed_s / 60)) %>%
  group_by(source_file, minute) %>%
  summarise(
    timestamp   = min(timestamp),
    elapsed_s   = min(elapsed_s),
    pct_opening = mean(pct_opening),
    .groups = "drop"
  )

pct_thresholds <- gaping %>% distinct(source_file, threshold_pct_opening)

p_pct_time <- ggplot(pct_by_min, aes(x = .data[[x_var]], y = pct_opening)) +
  geom_line(linewidth = 0.3, colour = "steelblue") +
  geom_hline(
    data = pct_thresholds,
    aes(yintercept = threshold_pct_opening), linetype = "dashed", colour = "firebrick"
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  facet_wrap(~source_file, scales = "free_x", ncol = 1) +
  labs(
    title = "Percent opening over time (1-minute means)",
    subtitle = "0% = closed baseline, 100% = widest opening seen; red dashed = open/closed threshold",
    x = x_lab, y = "Opening (%)"
  ) +
  theme_minimal()

p_pct_time

Percent opening distribution

How much of the recording was spent at each degree of opening (1 s samples). A clean two-state animal would show two peaks, at 0% and near 100%.

Code 
p_pct_dist <- ggplot(gaping, aes(x = pct_opening)) +
  geom_histogram(binwidth = 2.5, boundary = 0, fill = "steelblue", color = "white") +
  geom_vline(
    data = pct_thresholds,
    aes(xintercept = threshold_pct_opening), linetype = "dashed", colour = "firebrick"
  ) +
  geom_text(
    data = pct_thresholds,
    aes(x = threshold_pct_opening, y = Inf,
        label = sprintf("threshold = %.0f%%", threshold_pct_opening)),
    hjust = -0.1, vjust = 1.5, size = 3.5, colour = "firebrick"
  ) +
  facet_wrap(~source_file, ncol = 1) +
  labs(
    title = "Distribution of percent opening",
    subtitle = "Red dashed = open/closed detection threshold on the same scale",
    x = "Opening (%)", y = "Count of samples"
  ) +
  theme_minimal()

p_pct_dist

NOTES / NEXT STEPS

  • Gape distance depends on the magnet count, not just the board. The UID check ties a log to the right unit; nothing in the data ties it to the right magnet count. A wrong-but-calibrated count passes silently and compresses every distance, so treat the magnet-count cross-check table as a step to actually read, not a formality - and record the closed-shell magnet-to-sensor distance at deployment so there is something to check it against.
  • Distances are only as good as the calibration range. rmse_mm in the matched-curve table is the fit error over the bench range; the extrapolation check reports how much of the deployment fell outside it.
  • Threshold calibration. The noise and relative methods take the top 1 - baseline_quantile fraction of smoothed samples (default 10%) as “closed”. That only assumes the animal spends at least that fraction of the recording fully closed - not most of it. An animal can be open most of the time and still give a good baseline, as long as it closes fully now and then; check that the highest magnitudes are steady across the deployment. If the animal rarely or never closes fully, the baseline shifts toward open and the threshold with it. In that case, characterize known open/closed magnitudes on the bench and switch to threshold_method = "absolute".
  • Noise-based threshold (n_sd_threshold). noise_sd_uT comes from lagged differences of the smoothed signal, so it doesn’t depend on how often the animal is closed, and it shrinks as smooth_window_s grows. The default of 3 SD below the closed baseline is a starting point, not a calibrated value. Vary n_sd_threshold and smooth_window_s and confirm the detected event count and pct_samples_open are stable rather than an artifact of one setting.
  • Open/closed is a definition, not a measurement. If the percent-opening distribution is one broad hump rather than two peaks, there’s no natural break between open and closed, and the time-open fraction depends heavily on where the cutoff sits (check how pct_samples_open moves as n_sd_threshold is varied). In that case, report percent opening as the primary continuous measure and state the open/closed cutoff explicitly.
  • Event filtering. min_event_duration_s and min_gap_between_s set how aggressively blips are rejected and near-adjacent gapes are merged. Vary them and re-render to check the event count is stable rather than an artifact of one particular setting.
  • Timestamps. The onboard DS3231 RTC supplies absolute time directly, so no manual start-time entry is needed. If a session shows “NO (relative T+ elapsed time only)” in the DATA IMPORT summary, check that session’s RTC_INIT_ERROR row in eventlog.csv to see why the RTC wasn’t detected at boot.