OSCR

The Crunchometer, a low-cost, open-source acoustic analysis of feeding microstructure.

Code ↔ Paper

18 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 18 matches
  1. [1] § Results › The Crunchometer: an open-source sound-based method for studying the microstructure of feeding behavior ↔ crunchometer_python/crunchometer/core/analyzer.py, lines 38–87 · score 0.79 · dynamic thresholds, fixed threshold, frequency band, feeding behavior, crunch, 950 Hz
  2. [2] § Materials and methods › SVM model of the Crunchometer: automatic detection of biting from audio ↔ Xsembles/Linkage_JP.m, lines 1–60 · score 0.76 · Euclidean distances, hierarchical clustering, linkage, Ward, algorithm, vector
  3. [3] § Materials and methods › Calcium imaging experiments and analysis ↔ Xsembles_2P.m, lines 14–43 · score 0.70 · peak signal, noise ratio, motion correction, algorithm, Raw, neuronal
  4. [4] § Materials and methods › Tracking trajectory and running speed ↔ crunchometer_python/motion_energy.py, lines 14–89 · score 0.67 · binary mask, mass center, grayscale, speed, subtracted, mp4
  5. [5] § Materials and methods › Statistical analysis ↔ crunchometer_python/group_analysis.py, lines 456–507 · score 0.67 · Mann Whitney, post hoc, way ANOVA, metrics, Crunchometer
  6. [6] § Materials and methods › Calcium imaging experiments and analysis ↔ Deconvolution Caiman/foopsi_oasisAR2.m, lines 1–102 · score 0.66 · fluorescence traces, Neural activity, Calcium imaging, signal, events
  7. [7] § Materials and methods › Tracking trajectory and running speed ↔ crunchometer_python/crunchometer/core/tracker.py, lines 572–702 · score 0.65 · darker, Foreground, grayscale, morphological, static, speed
  8. [8] § Results › The Crunchometer: an open-source sound-based method for studying the microstructure of feeding behavior ↔ crunchometer_python/crunchometer/core/classifier.py, lines 23–53 · score 0.64 · high fat diet, gnawing behavior, motion energy, Feeding bouts, optional, noise
  9. [9] § Results › The Crunchometer system distinguishes meal patterns of fed and fasted mice ↔ crunchometer_python/group_analysis.py, lines 1121–1184 · score 0.62 · Inter Bout Interval, cumulative feeding, feeding rate, hour, SEM, IBI
  10. [10] § Results › Semaglutide suppresses feeding and reduces preference for a high-fat diet ↔ crunchometer_python/group_analysis.py, lines 456–507 · score 0.61 · Mann Whitney, post hoc, way ANOVA, treatment, Crunchometer
  11. [11] § Materials and methods › Identification of Chow vs. HFD consumption ↔ crunchometer_python/motion_energy.py, lines 167–226 · score 0.61 · Motion energy, HFD ROI, Chow ROI, ROIs, classified, threshold
  12. [12] § Results › The Crunchometer: an open-source sound-based method for studying the microstructure of feeding behavior ↔ crunchometer_python/crunchometer/core/resnet_detector.py, lines 1–39 · score 0.58 · sound bite events, chewing, spectrum, feeding bout, pipeline, noise
  13. [13] § Materials and methods › SVM model of the Crunchometer: automatic detection of biting from audio ↔ crunchometer_python/crunchometer/core/svm_detector.py, lines 252–335 · score 0.57 · frequency resolution, spectral, nfft, SVM, segments, overlapping
  14. [14] § Results › Event-locked alignment of calcium recordings to acoustically detected feeding bouts via the Crunchometer ↔ Xsembles_2P_Viewer.m, lines 833–920 · score 0.57 · OFFsemble, ONsemble, population, ensemble, activity, events
  15. [15] § Materials and methods › Signal processing pipeline of the Crunchometer ↔ crunchometer_python/crunchometer/core/svm_detector.py, lines 252–335 · score 0.55 · power spectrum, audio signals, overlap, spectrogram, window, bin
  16. [16] § Results › The Crunchometer system distinguishes meal patterns of fed and fasted mice ↔ crunchometer_python/group_analysis.py, lines 1–60 · score 0.55 · Mann Whitney, feeding rate, ANOVA, IBIs, feeding bouts, interaction
  17. [17] § Materials and methods › Signal processing pipeline of the Crunchometer ↔ crunchometer_python/crunchometer/core/audio.py, lines 67–132 · score 0.54 · 500–950 Hz, frequency band, sounds, spectrogram, power, bin
  18. [18] § Materials and methods › Identification of Chow vs. HFD consumption ↔ MoussionEnergy.m, lines 1–59 · score 0.54 · MoussionEnergy, Motion energy, ROI, frames, video

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,783 lines · 76 KB · MIT · 4 matches

  1. #!/usr/bin/env python3
  2. """
  3. Crunchometer Group Analysis v5.0
  4. =================================
  5. Complete pipeline for group-level analysis of feeding behavior data
  6. generated by CrunchometerV2.
  8. Ported from Crunchometer_Analysis_v5.m (MATLAB, 1710 lines)
  10. Features:
  11. - Pre-scan experiment folder structure (groups + subjects)
  12. - Load ground truth intake data (table or matrix format)
  13. - Compute feeding rates, bout durations, IBIs, and meal clusters
  14. - Calculate 12 per-subject metrics
  15. - Automatic statistical testing (t-test / Mann-Whitney / ANOVA / Kruskal-Wallis)
  16. - Outlier detection (2-SD) and CV reporting
  17. - Generate 5 publication-ready matplotlib figures
  18. - Export 2 Excel files (Summary + DeepAnalysis, 9 sheets total)
  19. - Save analysis log for traceability
  21. Matrix_Data columns (0-indexed in Python):
  22. Col 0: Start time of feeding bout (seconds)
  23. Col 1: End time of feeding bout (seconds)
  24. Col 2: Duration of feeding bout (seconds)
  25. Col 3: Number of bins crossing detection threshold
  26. Col 4: Power (dB) of each bite/event
  27. Col 5: Classification code (0=Artifact, 1=Chow, 2=HFD, 3=Gnawing)
  29. Original Author: Benjamin Arroyo (2025)
  30. Enhanced by: Ranier Gutierrez and AI Assistant (2026)
  31. """
  33. import os
  34. import re
  35. import warnings
  36. from dataclasses import dataclass, field
  37. from datetime import datetime
  38. from typing import Callable, Dict, List, Optional, Tuple
  40. import matplotlib
  41. # Use QtAgg backend when available (for popup figures in PyQt GUI),
  42. # fall back to Agg (non-interactive) for headless / scripted usage.
  43. try:
  44. matplotlib.use('QtAgg')
  45. except Exception:
  46. matplotlib.use('Agg')
  47. import matplotlib.pyplot as plt
  48. import matplotlib.patches as mpatches
  49. import numpy as np
  50. import pandas as pd
  51. from scipy.io import loadmat
  52. from scipy.optimize import curve_fit
  53. from scipy import stats
  56. # =========================================================================
  57. # CONFIGURATION
  58. # =========================================================================
  60. @dataclass
  61. class GroupAnalysisConfig:
  62. """Configuration parameters for group analysis pipeline."""
  63. bin_size: int = 60 # Time bin resolution (seconds)
  64. ibi_window: int = 600 # IBI averaging window (seconds)
  65. exp_duration: int = 7200 # Total experiment duration (seconds)
  66. meal_gap: int = 300 # Max gap between bouts in same meal (seconds)
  67. alpha_level: float = 0.05 # Significance threshold for statistics
  68. col_code: int = 5 # Classification column index (0-based)
  69. outlier_sd: float = 2.0 # SD threshold for outlier detection
  72. @dataclass
  73. class GroupData:
  74. """Container for a single treatment group's data."""
  75. name: str
  76. treat_idx: int # Index for legacy matrix GT format
  77. n_subs: int
  78. color: np.ndarray
  79. subjects: List[str]
  80. DM: List[Optional[np.ndarray]] = field(default_factory=list)
  81. intake: Optional[np.ndarray] = None # Feeding rate (all food)
  82. intake_chow: Optional[np.ndarray] = None # Chow only
  83. intake_hfd: Optional[np.ndarray] = None # HFD only
  84. bout_durations: List[Optional[np.ndarray]] = field(default_factory=list)
  85. ibi_values: List[Optional[np.ndarray]] = field(default_factory=list)
  86. meals: List[Optional[np.ndarray]] = field(default_factory=list)
  89. @dataclass
  90. class AnalysisResults:
  91. """Container for all analysis results."""
  92. figures: List[plt.Figure] = field(default_factory=list)
  93. output_path: str = ""
  94. summary_excel: str = ""
  95. deep_excel: str = ""
  96. log_file: str = ""
  99. # =========================================================================
  100. # DEFAULT COLORS
  101. # =========================================================================
  103. # Default group colors matching MATLAB palette
  104. DEFAULT_COLORS = np.array([
  105. [0.00, 0.50, 0.00], # Green
  106. [0.49, 0.18, 0.56], # Purple
  107. [0.85, 0.33, 0.10], # Orange
  108. [0.00, 0.45, 0.74], # Blue
  109. [0.30, 0.75, 0.93], # Light blue
  110. [0.64, 0.08, 0.18], # Dark red
  111. [0.47, 0.67, 0.19], # Olive
  112. [0.87, 0.49, 0.00], # Amber
  113. ])
  115. # Food type colors (consistent across all plots)
  116. FOOD_COLORS = {
  117. 'chow': np.array([0.93, 0.69, 0.13]), # Gold
  118. 'hfd': np.array([1.00, 0.00, 0.00]), # Red
  119. 'gnawing': np.array([0.00, 0.00, 0.00]), # Black
  120. }
  122. # Metrics definitions
  123. METRICS_LABELS = [
  124. 'Intake_g', 'Total_Time_s', 'Num_Bouts', 'Bout_Size_s',
  125. 'Latency_s', 'Mean_IBI_s', 'Num_Meals', 'Meal_Size_s',
  126. 'Pref_Ratio', 'Satiety_Ratio', 'Chow_Time_s', 'HFD_Time_s'
  127. ]
  129. METRICS_DISPLAY = [
  130. 'Intake (g)', 'Total Time (s)', '# Bouts', 'Bout Size (s)',
  131. 'Latency (s)', 'Mean IBI (s)', '# Meals', 'Meal Size (s)',
  132. 'Pref. Ratio', 'Satiety Ratio', 'Chow Time (s)', 'HFD Time (s)'
  133. ]
  135. N_METRICS = len(METRICS_LABELS)
  138. # =========================================================================
  139. # MAIN ANALYZER CLASS
  140. # =========================================================================
  142. class GroupAnalyzer:
  143. """
  144. Main group analysis engine for Crunchometer data.
  146. Orchestrates the full pipeline: scan → load → process → stats → figures → export.
  148. Parameters
  149. ----------
  150. config : GroupAnalysisConfig
  151. Analysis configuration parameters.
  152. experiment_path : str
  153. Path to the experiment root folder.
  154. intake_path : str, optional
  155. Path to manual intake Excel/CSV file.
  156. """
  158. def __init__(self, config: GroupAnalysisConfig, experiment_path: str,
  159. intake_path: Optional[str] = None):
  160. self.cfg = config
  161. self.experiment_path = experiment_path
  162. self.intake_path = intake_path or ""
  164. # Derived values
  165. self.experiment_name = re.sub(r'\s+', '_',
  166. os.path.basename(experiment_path.rstrip(os.sep)))
  167. self.results_path = os.path.join(experiment_path, 'Results')
  168. self.n_bins = int(np.ceil(self.cfg.exp_duration / self.cfg.bin_size))
  170. # Data containers (populated during processing)
  171. self.groups: List[GroupData] = []
  172. self.treatment_names: List[str] = []
  173. self.num_treatments: int = 0
  174. self.col: np.ndarray = np.array([]) # Group color matrix
  175. self.c_ibi_all: Optional[np.ndarray] = None # IBI temporal matrix
  177. # Results containers
  178. self.all_metrics: List[np.ndarray] = []
  179. self.stats_results: List[dict] = []
  180. self.outlier_flags: List[np.ndarray] = []
  181. self.cv_table: Optional[np.ndarray] = None
  182. self.GT_table = None
  183. self.has_intake: bool = False
  184. self.gt_format: str = 'none'
  186. # Log entries for traceability
  187. self.log_entries: List[str] = []
  189. # -----------------------------------------------------------------
  190. # STEP 1: SCAN EXPERIMENT STRUCTURE
  191. # -----------------------------------------------------------------
  192. def scan_experiment(self) -> None:
  193. """
  194. Discover treatment groups and subjects from folder structure.
  196. Expects: experiment_path / GroupName / Subject.mat
  197. Excludes: '.', '..', 'Results', '__pycache__'
  198. """
  199. self._log(f"Crunchometer Analysis v5.0 - {datetime.now():%Y-%m-%d %H:%M:%S}")
  200. self._log(f"Experiment: {self.experiment_name}")
  201. self._log(f"Path: {self.experiment_path}")
  202. self._log(f"Config: BinSize={self.cfg.bin_size}, IBIWindow={self.cfg.ibi_window}, "
  203. f"ExpDuration={self.cfg.exp_duration}, MealGap={self.cfg.meal_gap}")
  205. # Find treatment subfolders
  206. items = sorted(os.listdir(self.experiment_path))
  207. sub_folders = [
  208. d for d in items
  209. if os.path.isdir(os.path.join(self.experiment_path, d))
  210. and d not in ('.', '..', 'Results', '__pycache__')
  211. ]
  213. self.treatment_names = sub_folders
  214. self.num_treatments = len(sub_folders)
  216. if self.num_treatments == 0:
  217. raise ValueError(f"No treatment group folders found in: {self.experiment_path}")
  219. # Pre-scan all groups to find subjects
  220. subject_names = []
  221. max_subs = 0
  223. for i, treatment in enumerate(self.treatment_names):
  224. treat_path = os.path.join(self.experiment_path, treatment)
  225. mat_files = sorted([f for f in os.listdir(treat_path) if f.endswith('.mat')])
  226. subject_names.append(mat_files)
  227. n_subs = len(mat_files)
  228. if n_subs > max_subs:
  229. max_subs = n_subs
  230. self._log(f" {treatment}: {n_subs} subjects")
  232. # Assign colors
  233. if self.num_treatments <= len(DEFAULT_COLORS):
  234. self.col = DEFAULT_COLORS[:self.num_treatments]
  235. else:
  236. # Generate evenly spaced colors from a colormap
  237. cmap = plt.cm.get_cmap('tab10', self.num_treatments)
  238. self.col = np.array([cmap(i)[:3] for i in range(self.num_treatments)])
  240. # Create Results output folder
  241. os.makedirs(self.results_path, exist_ok=True)
  243. # Initialize group data structures
  244. self.groups = []
  245. for i, treatment in enumerate(self.treatment_names):
  246. n_s = len(subject_names[i])
  247. grp = GroupData(
  248. name=treatment,
  249. treat_idx=i,
  250. n_subs=n_s,
  251. color=self.col[i],
  252. subjects=subject_names[i],
  253. DM=[None] * n_s,
  254. intake=np.zeros((n_s, self.n_bins)),
  255. intake_chow=np.zeros((n_s, self.n_bins)),
  256. intake_hfd=np.zeros((n_s, self.n_bins)),
  257. bout_durations=[None] * n_s,
  258. ibi_values=[None] * n_s,
  259. meals=[None] * n_s,
  260. )
  261. self.groups.append(grp)
  263. # Pre-allocate IBI temporal matrix
  264. self.c_ibi_all = np.full((max_subs, self.cfg.exp_duration, self.num_treatments), np.nan)
  266. self._log(f"Groups found: {self.num_treatments}")
  268. # -----------------------------------------------------------------
  269. # STEP 2: LOAD GROUND TRUTH
  270. # -----------------------------------------------------------------
  271. def load_ground_truth(self) -> None:
  272. """
  273. Load manual intake ground truth data from Excel/CSV.
  275. Supports two formats:
  276. - 'table': Excel with Subject / Treatment / Intake_g columns (preferred)
  277. - 'matrix': Legacy numeric matrix (positional matching)
  278. """
  279. self.GT_table = None
  280. self.has_intake = False
  281. self.gt_format = 'none'
  283. manual_path = self.intake_path
  285. # Auto-search if no path provided
  286. if not manual_path:
  287. keywords = ['intake', 'manual', 'ground', 'gt', 'consumo', 'weight']
  288. for f in os.listdir(self.experiment_path):
  289. fpath = os.path.join(self.experiment_path, f)
  290. if os.path.isdir(fpath):
  291. continue
  292. _, ext = os.path.splitext(f)
  293. if ext.lower() == '.mat':
  294. continue # Skip .mat files
  295. fname_lower = f.lower()
  296. if any(kw in fname_lower for kw in keywords):
  297. manual_path = fpath
  298. break
  300. if not manual_path or not os.path.isfile(manual_path):
  301. self._log("Ground truth: NOT LOADED")
  302. return
  304. # Load the file
  305. try:
  306. _, ext = os.path.splitext(manual_path)
  308. if ext.lower() == '.mat':
  309. # Verify it's NOT a Crunchometer data file
  310. mat_data = loadmat(manual_path)
  311. if 'Matrix_Data' in mat_data:
  312. warnings.warn(f"File {manual_path} contains Matrix_Data. Skipping as intake file.")
  313. return
  314. # Use first non-system variable
  315. for key in mat_data:
  316. if not key.startswith('__'):
  317. self.GT_table = mat_data[key]
  318. self.gt_format = 'matrix'
  319. break
  320. else:
  321. # Try reading as table (preferred format)
  322. try:
  323. df = pd.read_excel(manual_path) if ext.lower() in ('.xlsx', '.xls') \
  324. else pd.read_csv(manual_path)
  325. required_cols = {'Subject', 'Treatment', 'Intake_g'}
  326. if required_cols.issubset(set(df.columns)):
  327. self.GT_table = df
  328. self.gt_format = 'table'
  329. else:
  330. # Fall back to matrix format
  331. self.GT_table = df.values
  332. self.gt_format = 'matrix'
  333. except Exception:
  334. # Last resort: try reading as numeric matrix
  335. df = pd.read_excel(manual_path) if ext.lower() in ('.xlsx', '.xls') \
  336. else pd.read_csv(manual_path)
  337. self.GT_table = df.values
  338. self.gt_format = 'matrix'
  340. self.has_intake = True
  341. self._log(f"Ground truth loaded: {manual_path} (format: {self.gt_format})")
  343. except Exception as e:
  344. warnings.warn(f"Error reading intake file: {e}")
  346. # -----------------------------------------------------------------
  347. # STEP 3: PROCESS ALL DATA
  348. # -----------------------------------------------------------------
  349. def process_all(self) -> None:
  350. """
  351. Load .mat files for each group/subject and compute:
  352. - Feeding rate time series (all food, chow-only, HFD-only)
  353. - Bout durations
  354. - Inter-bout intervals (IBIs)
  355. - Meal clusters
  356. """
  357. for i, grp in enumerate(self.groups):
  358. treat_path = os.path.join(self.experiment_path, grp.name)
  360. for j in range(grp.n_subs):
  361. f_path = os.path.join(treat_path, grp.subjects[j])
  362. try:
  363. mat = loadmat(f_path)
  364. if 'Matrix_Data' not in mat:
  365. self._log(f"WARNING: {grp.name}/{grp.subjects[j]} - no Matrix_Data")
  366. continue
  368. mdata = mat['Matrix_Data'].astype(float)
  370. # Validate Matrix_Data
  371. mdata = self._validate_matrix_data(mdata, grp.subjects[j])
  372. if mdata is None or len(mdata) == 0:
  373. continue
  375. grp.DM[j] = mdata
  376. code_col = mdata[:, self.cfg.col_code]
  378. # --- A: Feeding rate (combined, chow-only, hfd-only) ---
  379. for food_type in range(3): # 0=all, 1=chow, 2=hfd
  380. if food_type == 0:
  381. valid_codes = [1, 2]
  382. elif food_type == 1:
  383. valid_codes = [1]
  384. else:
  385. valid_codes = [2]
  387. valid_rows = np.isin(code_col, valid_codes)
  388. valid_evs = mdata[valid_rows][:, [0, 2]] # [start, duration]
  390. if len(valid_evs) > 0:
  391. # Force minimum 1s duration
  392. valid_evs[valid_evs[:, 1] <= 0, 1] = 1
  393. xt = np.zeros(self.cfg.exp_duration, dtype=float)
  394. for dt in range(len(valid_evs)):
  395. s_idx = max(0, int(round(valid_evs[dt, 0])))
  396. e_idx = min(self.cfg.exp_duration,
  397. s_idx + int(round(valid_evs[dt, 1])))
  398. xt[s_idx:e_idx] = 1
  400. for b in range(self.n_bins):
  401. idx_start = self.cfg.bin_size * b
  402. idx_end = min((b + 1) * self.cfg.bin_size,
  403. self.cfg.exp_duration)
  404. binned_val = np.sum(xt[idx_start:idx_end])
  405. if food_type == 0:
  406. grp.intake[j, b] = binned_val
  407. elif food_type == 1:
  408. grp.intake_chow[j, b] = binned_val
  409. else:
  410. grp.intake_hfd[j, b] = binned_val
  412. # --- B: Bout durations (from column 2 = duration) ---
  413. feed_rows = np.isin(code_col, [1, 2])
  414. if np.any(feed_rows):
  415. grp.bout_durations[j] = mdata[feed_rows, 2]
  417. # --- C: Inter-Bout Intervals ---
  418. det_f = mdata[feed_rows][:, 0:2] # [start, end]
  419. if len(det_f) > 1:
  420. ibis = det_f[1:, 0] - det_f[:-1, 1]
  421. grp.ibi_values[j] = ibis
  423. # Store in temporal matrix for IBI time course
  424. mid_times = det_f[:-1, 1] + ibis / 2
  425. idx_t = np.clip(np.round(mid_times).astype(int),
  426. 0, self.cfg.exp_duration - 1)
  427. for k in range(len(idx_t)):
  428. self.c_ibi_all[j, idx_t[k], i] = ibis[k]
  430. # --- D: Meal clustering ---
  431. if np.any(feed_rows):
  432. grp.meals[j] = self._cluster_bouts_into_meals(
  433. det_f, self.cfg.meal_gap)
  435. except Exception as e:
  436. self._log(f"ERROR: {grp.name}/{grp.subjects[j]} - {e}")
  438. # -----------------------------------------------------------------
  439. # STEP 4: COMPUTE METRICS
  440. # -----------------------------------------------------------------
  441. def compute_metrics(self) -> None:
  442. """
  443. Calculate 12 per-subject metrics for each treatment group.
  445. Metrics: Intake_g, Total_Time_s, Num_Bouts, Bout_Size_s,
  446. Latency_s, Mean_IBI_s, Num_Meals, Meal_Size_s,
  447. Pref_Ratio, Satiety_Ratio, Chow_Time_s, HFD_Time_s
  448. """
  449. self.all_metrics = []
  450. for grp in self.groups:
  451. mvals = np.full((grp.n_subs, N_METRICS), np.nan)
  452. for s in range(grp.n_subs):
  453. mvals[s, :] = self._compute_subject_metrics(s, grp)
  454. self.all_metrics.append(mvals)
  456. # -----------------------------------------------------------------
  457. # STEP 5: STATISTICAL ANALYSIS
  458. # -----------------------------------------------------------------
  459. def run_statistics(self) -> None:
  460. """
  461. Run automatic statistical tests on all 12 metrics.
  463. Selects appropriate test based on number of groups,
  464. sample sizes, and normality assumptions:
  465. - 2 groups, normal: Two-sample t-test
  466. - 2 groups, non-normal: Mann-Whitney U
  467. - 3+ groups, normal: One-way ANOVA (+ Tukey post-hoc)
  468. - 3+ groups, non-normal: Kruskal-Wallis (+ Dunn post-hoc)
  469. """
  470. self.stats_results = []
  471. for m in range(N_METRICS):
  472. data_per_group = []
  473. for i in range(self.num_treatments):
  474. vals = self.all_metrics[i][:, m]
  475. data_per_group.append(vals[~np.isnan(vals)])
  477. result = self._run_auto_stats(data_per_group, METRICS_DISPLAY[m])
  478. self.stats_results.append(result)
  479. self._log(f" {METRICS_DISPLAY[m]}: {result['test_name']}, p={result['p_value']:.4f}")
  481. # -----------------------------------------------------------------
  482. # STEP 6: OUTLIER DETECTION
  483. # -----------------------------------------------------------------
  484. def detect_outliers(self) -> None:
  485. """
  486. Flag outliers (>2 SD from mean) and compute coefficient of variation.
  487. """
  488. self.outlier_flags = []
  489. self.cv_table = np.zeros((self.num_treatments, N_METRICS))
  491. for i in range(self.num_treatments):
  492. mvals = self.all_metrics[i]
  493. flags = np.zeros(mvals.shape, dtype=bool)
  494. for m in range(N_METRICS):
  495. vals = mvals[:, m]
  496. valid = vals[~np.isnan(vals)]
  497. if len(valid) >= 3:
  498. mu = np.mean(valid)
  499. sd = np.std(valid, ddof=1)
  500. flags[:, m] = np.abs(vals - mu) > self.cfg.outlier_sd * sd
  501. if mu != 0:
  502. self.cv_table[i, m] = (sd / abs(mu)) * 100
  503. n_out = int(np.sum(flags[:, m]))
  504. if n_out > 0:
  505. self._log(f"OUTLIER: {self.treatment_names[i]}/"
  506. f"{METRICS_DISPLAY[m]} - {n_out} subject(s)")
  507. self.outlier_flags.append(flags)
  509. # -----------------------------------------------------------------
  510. # STEP 7: FIGURE GENERATION
  511. # -----------------------------------------------------------------
  512. def generate_figures(self) -> List[plt.Figure]:
  513. """
  514. Generate 5 publication-ready matplotlib figures.
  516. Returns
  517. -------
  518. list of matplotlib.figure.Figure
  519. The 5 generated figures.
  520. """
  521. figs = []
  522. figs.append(self._fig1_ethograms_cumulative())
  523. figs.append(self._fig2_core_metrics())
  524. figs.append(self._fig3_food_type())
  525. figs.append(self._fig4_microstructure())
  526. figs.append(self._fig5_advanced())
  527. return figs
  529. # -----------------------------------------------------------------
  530. # STEP 8: EXCEL EXPORT
  531. # -----------------------------------------------------------------
  532. def export_excel(self) -> Tuple[str, str]:
  533. """
  534. Export results to 2 Excel files:
  535. - Results_Summary.xlsx (4 sheets: Individual, Group Stats, Tests, CV)
  536. - DeepAnalysis_Data.xlsx (5 sheets: Events, Cumulative, Rate, IBI, Meals)
  538. Returns
  539. -------
  540. tuple of (str, str)
  541. Paths to the two generated Excel files.
  542. """
  543. xls_summary = os.path.join(
  544. self.results_path, f"{self.experiment_name}_Results_Summary.xlsx")
  545. xls_deep = os.path.join(
  546. self.results_path, f"{self.experiment_name}_DeepAnalysis_Data.xlsx")
  548. self._export_summary_excel(xls_summary)
  549. self._export_deep_excel(xls_deep)
  551. return xls_summary, xls_deep
  553. # -----------------------------------------------------------------
  554. # STEP 9: SAVE LOG
  555. # -----------------------------------------------------------------
  556. def save_log(self) -> str:
  557. """Save analysis log to text file. Returns log file path."""
  558. log_file = os.path.join(
  559. self.results_path, f"{self.experiment_name}_AnalysisLog.txt")
  560. self._log(f"Analysis completed: {datetime.now():%Y-%m-%d %H:%M:%S}")
  561. with open(log_file, 'w') as f:
  562. for entry in self.log_entries:
  563. f.write(entry + '\n')
  564. return log_file
  566. # -----------------------------------------------------------------
  567. # FULL PIPELINE ORCHESTRATOR
  568. # -----------------------------------------------------------------
  569. def run_full_pipeline(self, progress_callback: Optional[Callable] = None) -> AnalysisResults:
  570. """
  571. Run the complete analysis pipeline.
  573. Parameters
  574. ----------
  575. progress_callback : callable, optional
  576. Function(current, total, message) for progress updates.
  578. Returns
  579. -------
  580. AnalysisResults
  581. Container with figures, file paths, etc.
  582. """
  583. results = AnalysisResults()
  585. def _progress(step, total, msg):
  586. if progress_callback:
  587. progress_callback(step, total, msg)
  589. total_steps = 9
  591. # Step 1: Scan experiment structure
  592. _progress(1, total_steps, "Scanning experiment structure...")
  593. self.scan_experiment()
  595. # Step 2: Load ground truth
  596. _progress(2, total_steps, "Loading ground truth data...")
  597. self.load_ground_truth()
  599. # Step 3: Process all data
  600. _progress(3, total_steps, "Processing feeding data...")
  601. self.process_all()
  603. # Step 4: Compute metrics
  604. _progress(4, total_steps, "Computing per-subject metrics...")
  605. self.compute_metrics()
  607. # Step 5: Statistics
  608. _progress(5, total_steps, "Running statistical tests...")
  609. self.run_statistics()
  611. # Step 6: Outlier detection
  612. _progress(6, total_steps, "Detecting outliers...")
  613. self.detect_outliers()
  615. # Step 7: Generate figures
  616. _progress(7, total_steps, "Generating figures...")
  617. results.figures = self.generate_figures()
  619. # Step 8: Export Excel
  620. _progress(8, total_steps, "Exporting Excel reports...")
  621. results.summary_excel, results.deep_excel = self.export_excel()
  623. # Step 9: Save log
  624. _progress(9, total_steps, "Saving analysis log...")
  625. results.log_file = self.save_log()
  626. results.output_path = self.results_path
  628. return results
  630. # =====================================================================
  631. # PRIVATE HELPER METHODS
  632. # =====================================================================
  634. def _log(self, message: str) -> None:
  635. """Append a message to the analysis log."""
  636. self.log_entries.append(message)
  638. def _validate_matrix_data(self, mdata: np.ndarray, filename: str) -> Optional[np.ndarray]:
  639. """
  640. Validate Matrix_Data integrity.
  641. Returns cleaned array or None if invalid.
  642. """
  643. if mdata.ndim != 2 or mdata.shape[1] <= self.cfg.col_code:
  644. self._log(f"WARNING: {filename} has fewer than {self.cfg.col_code + 1} columns")
  645. return None
  647. # Remove rows with negative times
  648. bad_rows = (mdata[:, 0] < 0) | (mdata[:, 1] < 0)
  649. if np.any(bad_rows):
  650. self._log(f"WARNING: {filename}: Removed {int(np.sum(bad_rows))} rows "
  651. "with negative timestamps")
  652. mdata = mdata[~bad_rows]
  654. # Remove rows where start > end
  655. if len(mdata) > 0:
  656. bad_order = mdata[:, 0] > mdata[:, 1]
  657. if np.any(bad_order):
  658. self._log(f"WARNING: {filename}: Removed {int(np.sum(bad_order))} rows "
  659. "where start > end")
  660. mdata = mdata[~bad_order]
  662. if len(mdata) == 0:
  663. self._log(f"WARNING: {filename}: No valid data after validation")
  664. return None
  666. # Recalculate duration from col1 - col0
  667. mdata[:, 2] = mdata[:, 1] - mdata[:, 0]
  669. # Remove rows with zero or negative duration
  670. bad_dur = mdata[:, 2] <= 0
  671. if np.any(bad_dur):
  672. self._log(f"WARNING: {filename}: Removed {int(np.sum(bad_dur))} rows "
  673. "with duration <= 0")
  674. mdata = mdata[~bad_dur]
  676. # Remove events that START after the experiment duration
  677. # (these are entirely outside the analysis window)
  678. if len(mdata) > 0:
  679. after_exp = mdata[:, 0] >= self.cfg.exp_duration
  680. if np.any(after_exp):
  681. self._log(f"WARNING: {filename}: Removed {int(np.sum(after_exp))} "
  682. f"events starting after experiment duration "
  683. f"({self.cfg.exp_duration}s)")
  684. mdata = mdata[~after_exp]
  686. # Clip events that straddle the experiment boundary
  687. # (start < exp_duration < end)
  688. if len(mdata) > 0:
  689. over_dur = mdata[:, 1] > self.cfg.exp_duration
  690. if np.any(over_dur):
  691. mdata[over_dur, 1] = self.cfg.exp_duration
  692. mdata[over_dur, 2] = mdata[over_dur, 1] - mdata[over_dur, 0]
  693. # Safety: remove any that became <= 0 after clipping
  694. bad_after_clip = mdata[:, 2] <= 0
  695. if np.any(bad_after_clip):
  696. mdata = mdata[~bad_after_clip]
  698. if len(mdata) == 0:
  699. self._log(f"WARNING: {filename}: No valid data after validation")
  700. return None
  702. return mdata
  704. def _compute_subject_metrics(self, s: int, group: GroupData) -> np.ndarray:
  705. """
  706. Calculate all 12 metrics for a single subject.
  707. Returns a row vector of length N_METRICS.
  708. """
  709. mvals = np.full(N_METRICS, np.nan)
  710. data = group.DM[s]
  711. if data is None or data.ndim != 2 or data.shape[1] <= self.cfg.col_code:
  712. return mvals
  714. code_col = data[:, self.cfg.col_code]
  715. feed_rows = np.isin(code_col, [1, 2])
  716. evs = data[feed_rows]
  718. # 1: Intake (g) from ground truth
  719. if self.has_intake:
  720. weights = self._get_gt_for_group(group, group.treat_idx)
  721. if s < len(weights):
  722. mvals[0] = weights[s]
  724. # 2: Total feeding time (s)
  725. mvals[1] = np.sum(group.intake[s, :])
  727. # 3: Number of bouts
  728. mvals[2] = len(evs)
  730. # 4: Mean bout size (s)
  731. if len(evs) > 0:
  732. mvals[3] = np.mean(evs[:, 2])
  734. # 5: Latency to first bout (s)
  735. if len(evs) > 0:
  736. mvals[4] = evs[0, 0]
  738. # 6: Mean IBI (s)
  739. ibi_vals = group.ibi_values[s]
  740. if ibi_vals is not None and len(ibi_vals) > 0:
  741. mvals[5] = np.mean(ibi_vals)
  743. # 7: Number of meals
  744. meals = group.meals[s]
  745. if meals is not None and len(meals) > 0:
  746. mvals[6] = len(meals)
  748. # 8: Mean meal size (s)
  749. if meals is not None and len(meals) > 0:
  750. mvals[7] = np.mean(meals[:, 2])
  752. # 9: Preference ratio (HFD / total)
  753. chow_time = np.sum(group.intake_chow[s, :])
  754. hfd_time = np.sum(group.intake_hfd[s, :])
  755. total_time = chow_time + hfd_time
  756. if total_time > 0:
  757. mvals[8] = hfd_time / total_time
  759. # 10: Satiety ratio (mean post-bout IBI / mean bout duration)
  760. if (ibi_vals is not None and len(ibi_vals) > 0 and
  761. len(evs) > 1):
  762. mean_ibi = np.mean(ibi_vals)
  763. mean_bout = np.mean(evs[:, 2])
  764. if mean_bout > 0:
  765. mvals[9] = mean_ibi / mean_bout
  767. # 11: Chow time (s)
  768. mvals[10] = np.sum(group.intake_chow[s, :])
  770. # 12: HFD time (s)
  771. mvals[11] = np.sum(group.intake_hfd[s, :])
  773. return mvals
  775. def _get_gt_for_group(self, group: GroupData, treat_idx: int) -> np.ndarray:
  776. """
  777. Extract ground truth intake values for a specific group.
  778. Supports both 'table' (name-matched) and 'matrix' (positional) formats.
  779. """
  780. n_s = group.n_subs
  781. weights = np.full(n_s, np.nan)
  783. if self.GT_table is None:
  784. return weights
  786. if self.gt_format == 'table':
  787. df = self.GT_table
  788. for s in range(n_s):
  789. sname = os.path.splitext(group.subjects[s])[0]
  790. mask = (df['Subject'].astype(str) == sname) & \
  791. (df['Treatment'].astype(str) == group.name)
  792. matches = df.loc[mask, 'Intake_g']
  793. if len(matches) > 0:
  794. weights[s] = matches.iloc[0]
  796. elif self.gt_format == 'matrix':
  797. gt = self.GT_table
  798. if treat_idx < gt.shape[1]:
  799. n_valid = min(n_s, gt.shape[0])
  800. weights[:n_valid] = gt[:n_valid, treat_idx]
  802. return weights
  804. @staticmethod
  805. def _cluster_bouts_into_meals(bout_intervals: np.ndarray,
  806. meal_gap: float) -> Optional[np.ndarray]:
  807. """
  808. Group feeding bouts into meals based on inter-bout gaps.
  810. Parameters
  811. ----------
  812. bout_intervals : ndarray, shape (N, 2)
  813. [start_time, end_time] for each bout.
  814. meal_gap : float
  815. Maximum gap (seconds) between bouts in the same meal.
  817. Returns
  818. -------
  819. ndarray or None
  820. Shape (M, 4): [meal_start, meal_end, meal_duration, n_bouts]
  821. """
  822. if bout_intervals is None or len(bout_intervals) == 0:
  823. return None
  825. # Sort by start time
  826. sorted_bouts = bout_intervals[bout_intervals[:, 0].argsort()]
  827. n = len(sorted_bouts)
  829. meals = []
  830. meal_start = sorted_bouts[0, 0]
  831. meal_end = sorted_bouts[0, 1]
  832. n_bouts_in_meal = 1
  834. for k in range(1, n):
  835. gap = sorted_bouts[k, 0] - meal_end
  836. if gap <= meal_gap:
  837. # Continue current meal
  838. meal_end = sorted_bouts[k, 1]
  839. n_bouts_in_meal += 1
  840. else:
  841. # Close current meal, start new one
  842. meals.append([meal_start, meal_end,
  843. meal_end - meal_start, n_bouts_in_meal])
  844. meal_start = sorted_bouts[k, 0]
  845. meal_end = sorted_bouts[k, 1]
  846. n_bouts_in_meal = 1
  848. # Close last meal
  849. meals.append([meal_start, meal_end,
  850. meal_end - meal_start, n_bouts_in_meal])
  852. return np.array(meals)
  854. def _run_auto_stats(self, data_per_group: List[np.ndarray],
  855. metric_name: str) -> dict:
  856. """
  857. Automatically select and run appropriate statistical test.
  859. Returns dict with: test_name, statistic, p_value, decision, note, posthoc
  860. """
  861. result = {
  862. 'test_name': 'None',
  863. 'statistic': np.nan,
  864. 'p_value': np.nan,
  865. 'decision': 'N/A',
  866. 'note': '',
  867. 'posthoc': []
  868. }
  870. # Filter out empty groups
  871. valid_groups = [g for g in data_per_group if len(g) >= 1]
  872. if len(valid_groups) < 2:
  873. result['test_name'] = 'Descriptive only'
  874. result['note'] = 'Fewer than 2 valid groups for comparison.'
  875. return result
  877. # Check minimum sample sizes
  878. min_n = min(len(g) for g in valid_groups)
  879. if min_n < 3:
  880. result['test_name'] = 'Descriptive only'
  881. result['note'] = f'Insufficient sample size (min n={min_n}). Need n>=3 per group.'
  882. return result
  884. # Check normality (Shapiro-Wilk as Python alternative to lillietest)
  885. is_normal = True
  886. for g in valid_groups:
  887. if len(g) >= 4:
  888. try:
  889. _, p_shapiro = stats.shapiro(g)
  890. if p_shapiro < 0.05:
  891. is_normal = False
  892. break
  893. except Exception:
  894. is_normal = False
  895. break
  897. # Check homogeneity of variance (Bartlett or Levene)
  898. homo_var = True
  899. if len(valid_groups) >= 2 and is_normal:
  900. try:
  901. _, p_bart = stats.bartlett(*valid_groups)
  902. if p_bart < 0.05:
  903. homo_var = False
  904. except Exception:
  905. homo_var = True
  907. # Select and run test
  908. n_active = len(valid_groups)
  910. if n_active == 2:
  911. g1, g2 = valid_groups[0], valid_groups[1]
  913. if is_normal and homo_var:
  914. stat_val, p_val = stats.ttest_ind(g1, g2)
  915. result['test_name'] = 'Two-sample t-test'
  916. result['statistic'] = float(stat_val)
  917. result['p_value'] = float(p_val)
  918. result['note'] = 'Normality: passed. Equal variance: passed.'
  919. else:
  920. stat_val, p_val = stats.mannwhitneyu(g1, g2, alternative='two-sided')
  921. result['test_name'] = 'Mann-Whitney U'
  922. result['statistic'] = float(stat_val)
  923. result['p_value'] = float(p_val)
  924. reasons = []
  925. if not is_normal:
  926. reasons.append('Normality failed')
  927. if not homo_var:
  928. reasons.append('Unequal variance')
  929. result['note'] = f"Non-parametric ({'; '.join(reasons)})."
  931. elif n_active > 2:
  932. if is_normal and homo_var:
  933. stat_val, p_val = stats.f_oneway(*valid_groups)
  934. result['test_name'] = 'One-way ANOVA'
  935. result['statistic'] = float(stat_val)
  936. result['p_value'] = float(p_val)
  937. result['note'] = 'Normality: passed. Equal variance: passed.'
  939. # Post-hoc Tukey if significant
  940. if p_val < 0.05:
  941. try:
  942. tukey = stats.tukey_hsd(*valid_groups)
  943. for ii in range(len(valid_groups)):
  944. for jj in range(ii + 1, len(valid_groups)):
  945. ph_p = tukey.pvalue[ii, jj]
  946. sig = '*' if ph_p < 0.05 else 'ns'
  947. result['posthoc'].append(
  948. f"{self.treatment_names[ii]} vs "
  949. f"{self.treatment_names[jj]}: p={ph_p:.4f} {sig}")
  950. except Exception:
  951. pass
  952. else:
  953. stat_val, p_val = stats.kruskal(*valid_groups)
  954. result['test_name'] = 'Kruskal-Wallis'
  955. result['statistic'] = float(stat_val)
  956. result['p_value'] = float(p_val)
  957. reasons = []
  958. if not is_normal:
  959. reasons.append('Normality failed')
  960. if not homo_var:
  961. reasons.append('Unequal variance')
  962. result['note'] = f"Non-parametric ({'; '.join(reasons)})."
  964. # Post-hoc Dunn (simplified pairwise Mann-Whitney with Bonferroni)
  965. if p_val < 0.05:
  966. n_comp = n_active * (n_active - 1) // 2
  967. for ii in range(len(valid_groups)):
  968. for jj in range(ii + 1, len(valid_groups)):
  969. try:
  970. _, ph_p = stats.mannwhitneyu(
  971. valid_groups[ii], valid_groups[jj],
  972. alternative='two-sided')
  973. # Bonferroni correction
  974. ph_p_adj = min(ph_p * n_comp, 1.0)
  975. sig = '*' if ph_p_adj < 0.05 else 'ns'
  976. result['posthoc'].append(
  977. f"{self.treatment_names[ii]} vs "
  978. f"{self.treatment_names[jj]}: p={ph_p_adj:.4f} {sig}")
  979. except Exception:
  980. pass
  982. # Decision
  983. if not np.isnan(result['p_value']):
  984. p = result['p_value']
  985. if p < 0.001:
  986. result['decision'] = 'Significant (p<0.001)'
  987. elif p < 0.01:
  988. result['decision'] = 'Significant (p<0.01)'
  989. elif p < 0.05:
  990. result['decision'] = 'Significant (p<0.05)'
  991. else:
  992. result['decision'] = 'Not significant'
  994. # Low N warning
  995. if min_n < 5:
  996. result['note'] += ' WARNING: Low sample size (n<5), interpret with caution.'
  998. return result
  1000. # =====================================================================
  1001. # FIGURE GENERATION METHODS
  1002. # =====================================================================
  1004. def _shaded_error_bar(self, ax, x, y, err, color, linewidth=1.5,
  1005. alpha=0.15, marker=None, markersize=4, label=None):
  1006. """
  1007. Plot mean line with shaded SEM region.
  1008. Returns the main line handle.
  1009. """
  1010. y = np.asarray(y).ravel()
  1011. err = np.asarray(err).ravel()
  1012. x = np.asarray(x).ravel()
  1014. kwargs = dict(color=color, linewidth=linewidth, label=label)
  1015. if marker:
  1016. kwargs['marker'] = marker
  1017. kwargs['markersize'] = markersize
  1018. kwargs['markerfacecolor'] = color
  1020. line, = ax.plot(x, y, **kwargs)
  1021. ax.fill_between(x, y - err, y + err, color=color, alpha=alpha)
  1022. return line
  1024. def _add_significance_bracket(self, ax, x1, x2, p_val):
  1025. """Draw a significance bracket between two bar positions."""
  1026. yl = ax.get_ylim()
  1027. y_top = yl[1] * 0.90
  1028. y_bracket = y_top * 0.95
  1030. ax.plot([x1, x1, x2, x2],
  1031. [y_top * 0.88, y_bracket, y_bracket, y_top * 0.88],
  1032. '-k', linewidth=1)
  1034. if p_val < 0.001:
  1035. stars = '***'
  1036. elif p_val < 0.01:
  1037. stars = '**'
  1038. elif p_val < 0.05:
  1039. stars = '*'
  1040. else:
  1041. stars = 'ns'
  1043. ax.text((x1 + x2) / 2, y_bracket * 1.02, stars,
  1044. ha='center', fontsize=10, fontweight='bold')
  1046. def _ethogram_plot(self, ax, matrix_data, sub_idx, duration, code_col_idx):
  1047. """
  1048. Draw color-coded patches for each behavioral event.
  1049. Code 1 (Chow) = Gold, Code 2 (HFD) = Red, Code 3 (Gnawing) = Black
  1050. """
  1051. y_center = -(sub_idx + 1) # Offset each subject downward
  1052. beh_codes = [1, 2, 3]
  1053. beh_colors = [
  1054. (0.93, 0.69, 0.13), # Gold (Chow)
  1055. (1.0, 0.0, 0.0), # Red (HFD)
  1056. (0.0, 0.0, 0.0), # Black (Gnawing)
  1057. ]
  1059. if matrix_data.shape[1] > code_col_idx:
  1060. codes = matrix_data[:, code_col_idx]
  1061. for k, code in enumerate(beh_codes):
  1062. rows = codes == code
  1063. if not np.any(rows):
  1064. continue
  1065. intervals = matrix_data[rows][:, 0:2]
  1066. for iv in intervals:
  1067. rect = mpatches.Rectangle(
  1068. (iv[0], y_center - 0.4), iv[1] - iv[0], 0.8,
  1069. facecolor=beh_colors[k], edgecolor='none')
  1070. ax.add_patch(rect)
  1072. def _save_figure(self, fig, fig_label: str) -> None:
  1073. """Save figure in PNG format to results folder."""
  1074. base_name = f"{self.experiment_name}_{fig_label}"
  1075. full_base = os.path.join(self.results_path, base_name)
  1076. fig.savefig(f"{full_base}.png", dpi=300, bbox_inches='tight',
  1077. facecolor='white')
  1079. # ----- FIGURE 1: Ethograms + Cumulative Feeding -----
  1080. def _fig1_ethograms_cumulative(self) -> plt.Figure:
  1081. """Generate Figure 1: Ethograms and cumulative feeding per group."""
  1082. # Compute global Y limit for cumulative plots
  1083. global_max_y = 0
  1084. for grp in self.groups:
  1085. if grp.intake is not None and grp.intake.size > 0:
  1086. cum_d = np.cumsum(grp.intake, axis=1)
  1087. n_s = cum_d.shape[0]
  1088. peak = np.max(np.mean(cum_d, axis=0) +
  1089. np.std(cum_d, axis=0, ddof=0) / np.sqrt(max(n_s, 1)))
  1090. if peak > global_max_y:
  1091. global_max_y = peak
  1092. if global_max_y == 0:
  1093. global_max_y = 10
  1094. global_ylim = (0, global_max_y * 1.15)
  1096. fig, axes = plt.subplots(2, self.num_treatments,
  1097. figsize=(5 * self.num_treatments, 6.5),
  1098. squeeze=False)
  1100. for i, grp in enumerate(self.groups):
  1101. # Top row: ethograms
  1102. ax_top = axes[0, i]
  1103. ax_top.set_title(grp.name.replace('_', ' '), color=grp.color,
  1104. fontweight='bold', fontsize=11)
  1105. valid_subs = 0
  1106. for sub in range(len(grp.DM)):
  1107. if grp.DM[sub] is not None:
  1108. self._ethogram_plot(ax_top, grp.DM[sub], sub,
  1109. self.cfg.exp_duration, self.cfg.col_code)
  1110. valid_subs = sub + 1
  1111. ax_top.set_xlim(0, self.cfg.exp_duration)
  1112. if valid_subs > 0:
  1113. ax_top.set_ylim(-valid_subs - 1, 0)
  1114. ax_top.set_yticks([])
  1115. ax_top.set_xlabel('Time (s)')
  1116. if i == 0:
  1117. ax_top.set_ylabel('Subjects')
  1118. ax_top.spines['top'].set_visible(False)
  1119. ax_top.spines['right'].set_visible(False)
  1121. # Bottom row: cumulative feeding
  1122. ax_bot = axes[1, i]
  1123. data_mat = grp.intake
  1124. if data_mat is not None and data_mat.size > 0:
  1125. cum_data = np.cumsum(data_mat, axis=1)
  1126. mu = np.mean(cum_data, axis=0)
  1127. sigma = np.std(cum_data, axis=0, ddof=0) / np.sqrt(max(cum_data.shape[0], 1))
  1128. x_axis = np.arange(1, cum_data.shape[1] + 1) * self.cfg.bin_size
  1129. x_hours = x_axis / 3600.0
  1130. self._shaded_error_bar(ax_bot, x_hours, mu, sigma, grp.color)
  1131. ax_bot.set_xlim(0, self.cfg.exp_duration / 3600.0)
  1132. ax_bot.set_ylim(global_ylim)
  1133. ax_bot.set_xlabel('Time (h)')
  1134. if i == 0:
  1135. ax_bot.set_ylabel('Cum. Feeding Time (s)')
  1136. ax_bot.spines['top'].set_visible(False)
  1137. ax_bot.spines['right'].set_visible(False)
  1139. fig.suptitle(f"{self.experiment_name.replace('_', ' ')} — Ethograms & Cumulative Feeding",
  1140. fontsize=13, fontweight='bold')
  1141. fig.tight_layout(rect=[0, 0, 1, 0.95])
  1142. self._save_figure(fig, 'Fig1_Ethogram')
  1143. return fig
  1145. # ----- FIGURE 2: Core Metrics -----
  1146. def _fig2_core_metrics(self) -> plt.Figure:
  1147. """Generate Figure 2: IBI time course, validation, feeding rate, bar charts."""
  1148. fig = plt.figure(figsize=(14, 9))
  1149. gs = fig.add_gridspec(3, 6, hspace=0.45, wspace=0.5)
  1151. # --- 2A: IBI Time Course ---
  1152. ax_ibi = fig.add_subplot(gs[0, 0:4])
  1153. n_bins_ibi = int(self.cfg.exp_duration // self.cfg.ibi_window)
  1154. handles_ibi = []
  1155. for i, grp in enumerate(self.groups):
  1156. raw_ibi = self.c_ibi_all[:, :, i]
  1157. ibi_binned_mean = np.zeros(n_bins_ibi)
  1158. ibi_binned_sem = np.zeros(n_bins_ibi)
  1159. x_ibi = np.zeros(n_bins_ibi)
  1161. for b in range(n_bins_ibi):
  1162. idx_s = b * self.cfg.ibi_window
  1163. idx_e = (b + 1) * self.cfg.ibi_window
  1164. block = raw_ibi[:, idx_s:idx_e]
  1165. x_ibi[b] = (idx_s + idx_e) / 2.0
  1166. vals = block[~np.isnan(block)]
  1167. if len(vals) > 0:
  1168. ibi_binned_mean[b] = np.mean(vals)
  1169. ibi_binned_sem[b] = np.std(vals, ddof=0) / np.sqrt(len(vals))
  1171. x_hours = x_ibi / 3600.0
  1172. h = ax_ibi.errorbar(x_hours, ibi_binned_mean, yerr=ibi_binned_sem,
  1173. fmt='-o', color=grp.color, markerfacecolor='w',
  1174. linewidth=1.2, markersize=5,
  1175. label=grp.name.replace('_', ' '))
  1176. handles_ibi.append(h)
  1178. ax_ibi.set_ylabel('IBI (s)')
  1179. ax_ibi.set_xlabel('Time (h)')
  1180. ax_ibi.set_title('Inter-Bout Interval', fontsize=10)
  1181. ax_ibi.set_xlim(0, self.cfg.exp_duration / 3600.0)
  1182. ax_ibi.legend(fontsize=7, frameon=False, loc='best')
  1183. ax_ibi.spines['top'].set_visible(False)
  1184. ax_ibi.spines['right'].set_visible(False)
  1186. # --- 2B: Validation (Regression) ---
  1187. ax_reg = fig.add_subplot(gs[0, 4:6])
  1188. if self.has_intake:
  1189. all_weights = []
  1190. all_times = []
  1191. all_colors_reg = []
  1192. for i, grp in enumerate(self.groups):
  1193. total_time = np.sum(grp.intake, axis=1)
  1194. weights = self._get_gt_for_group(grp, i)
  1195. n_valid = min(len(total_time), len(weights))
  1196. w = weights[:n_valid]
  1197. t = total_time[:n_valid]
  1198. keep = ~np.isnan(w) & ~np.isnan(t)
  1199. w = w[keep]
  1200. t = t[keep]
  1201. if len(w) > 0:
  1202. ax_reg.plot(w, t, 'o', color=grp.color,
  1203. markerfacecolor=grp.color, markersize=6,
  1204. label=grp.name.replace('_', ' '))
  1205. all_weights.extend(w)
  1206. all_times.extend(t)
  1208. all_weights = np.array(all_weights)
  1209. all_times = np.array(all_times)
  1210. valid_idx = ~np.isnan(all_weights) & ~np.isnan(all_times)
  1211. if np.any(valid_idx):
  1212. w_v = all_weights[valid_idx]
  1213. t_v = all_times[valid_idx]
  1214. # Linear regression
  1215. coeffs = np.polyfit(w_v, t_v, 1)
  1216. y_fit = np.polyval(coeffs, w_v)
  1217. ss_res = np.sum((t_v - y_fit) ** 2)
  1218. ss_tot = np.sum((t_v - np.mean(t_v)) ** 2)
  1219. r_sq = 1 - ss_res / ss_tot if ss_tot > 0 else 0
  1221. x_line = np.linspace(np.min(w_v), np.max(w_v), 100)
  1222. ax_reg.plot(x_line, np.polyval(coeffs, x_line), '-r', linewidth=1.5)
  1223. ax_reg.text(0.05, 0.95, f'R² = {r_sq:.3f}',
  1224. transform=ax_reg.transAxes, color='r',
  1225. fontsize=10, fontweight='bold', va='top')
  1227. ax_reg.set_xlabel('Manual Intake (g)')
  1228. ax_reg.set_ylabel('Feeding Time (s)')
  1229. ax_reg.set_title('Regression Validation', fontsize=10)
  1230. ax_reg.legend(fontsize=7, frameon=False, loc='best')
  1231. ax_reg.grid(True, alpha=0.3)
  1232. else:
  1233. ax_reg.text(0.5, 0.5, 'No Intake Data', ha='center', va='center',
  1234. fontsize=11, color='gray', transform=ax_reg.transAxes)
  1235. ax_reg.set_axis_off()
  1237. # --- 2B2: Bland-Altman ---
  1238. ax_ba = fig.add_subplot(gs[1, 4:6])
  1239. if self.has_intake and len(all_weights) > 0:
  1240. w_v = np.array(all_weights)
  1241. t_v = np.array(all_times)
  1242. valid = ~np.isnan(w_v) & ~np.isnan(t_v)
  1243. if np.sum(valid) > 1:
  1244. w_z = (w_v[valid] - np.mean(w_v[valid])) / (np.std(w_v[valid], ddof=1) + 1e-10)
  1245. t_z = (t_v[valid] - np.mean(t_v[valid])) / (np.std(t_v[valid], ddof=1) + 1e-10)
  1246. avg_ba = (w_z + t_z) / 2
  1247. diff_ba = t_z - w_z
  1248. mean_diff = np.mean(diff_ba)
  1249. sd_diff = np.std(diff_ba, ddof=1)
  1251. # Build color array per point
  1252. color_list = []
  1253. for i, grp in enumerate(self.groups):
  1254. total_time = np.sum(grp.intake, axis=1)
  1255. weights = self._get_gt_for_group(grp, i)
  1256. n_valid = min(len(total_time), len(weights))
  1257. w = weights[:n_valid]
  1258. t = total_time[:n_valid]
  1259. keep = ~np.isnan(w) & ~np.isnan(t)
  1260. color_list.extend([grp.color] * int(np.sum(keep)))
  1261. color_arr = np.array(color_list) if color_list else np.array([[0, 0, 0]])
  1263. ax_ba.scatter(avg_ba, diff_ba, c=color_arr[:len(avg_ba)],
  1264. s=50, edgecolors='k', linewidths=0.5)
  1265. ax_ba.axhline(mean_diff, color='b', linewidth=1.2)
  1266. ax_ba.axhline(mean_diff + 1.96 * sd_diff, color='r',
  1267. linestyle='--', linewidth=1)
  1268. ax_ba.axhline(mean_diff - 1.96 * sd_diff, color='r',
  1269. linestyle='--', linewidth=1)
  1270. ax_ba.set_xlabel('Mean (z-score)')
  1271. ax_ba.set_ylabel('Difference (z-score)')
  1272. ax_ba.set_title('Bland-Altman (normalized)', fontsize=10)
  1273. ax_ba.text(0.05, 0.95, f'Bias = {mean_diff:.3f}',
  1274. transform=ax_ba.transAxes, color='b', fontsize=8,
  1275. fontweight='bold', va='top')
  1276. ax_ba.grid(True, alpha=0.3)
  1277. else:
  1278. ax_ba.text(0.5, 0.5, 'No Intake Data', ha='center', va='center',
  1279. fontsize=11, color='gray', transform=ax_ba.transAxes)
  1280. ax_ba.set_axis_off()
  1282. # --- 2C: Feeding Rate Time Course ---
  1283. ax_tc = fig.add_subplot(gs[1, 0:4])
  1284. for i, grp in enumerate(self.groups):
  1285. raw = grp.intake
  1286. if raw is None or raw.size == 0:
  1287. continue
  1288. group_factor = 10
  1289. new_cols = raw.shape[1] // group_factor
  1290. if new_cols == 0:
  1291. continue
  1292. binned_vals = np.zeros((raw.shape[0], new_cols))
  1293. for k in range(new_cols):
  1294. binned_vals[:, k] = np.sum(
  1295. raw[:, k * group_factor:(k + 1) * group_factor], axis=1)
  1296. mu = np.nanmean(binned_vals, axis=0)
  1297. err = np.nanstd(binned_vals, axis=0, ddof=0) / np.sqrt(max(binned_vals.shape[0], 1))
  1298. x_mins = np.linspace(0, self.cfg.exp_duration / 60.0, len(mu))
  1299. self._shaded_error_bar(ax_tc, x_mins, mu, err, grp.color,
  1300. marker='o', markersize=4,
  1301. label=grp.name.replace('_', ' '))
  1303. ax_tc.set_xlabel('Time (min)')
  1304. ax_tc.set_ylabel('Feeding Time (s / 10 min)')
  1305. ax_tc.set_title('Feeding Rate Time Course', fontsize=10)
  1306. ax_tc.legend(fontsize=7, frameon=False, loc='best')
  1307. ax_tc.spines['top'].set_visible(False)
  1308. ax_tc.spines['right'].set_visible(False)
  1310. # --- 2D: Bar charts for core metrics (6 panels) ---
  1311. core_metrics_idx = [0, 1, 2, 3, 4, 5] # 0-based
  1312. for m_i in range(6):
  1313. ax = fig.add_subplot(gs[2, m_i])
  1314. m = core_metrics_idx[m_i]
  1315. for ii in range(self.num_treatments):
  1316. vals = self.all_metrics[ii][:, m]
  1317. vals_clean = vals[~np.isnan(vals)]
  1318. if len(vals_clean) == 0:
  1319. continue
  1320. ax.bar(ii + 1, np.mean(vals_clean), color=self.col[ii],
  1321. edgecolor='none', alpha=0.6)
  1322. ax.errorbar(ii + 1, np.mean(vals_clean),
  1323. yerr=np.std(vals_clean, ddof=0) / np.sqrt(len(vals_clean)),
  1324. fmt='none', ecolor='k', linewidth=1)
  1325. # Scatter individual points
  1326. jitter = (np.random.rand(len(vals_clean)) - 0.5) * 0.3
  1327. ax.plot(ii + 1 + jitter, vals_clean, '.', color='0.3', markersize=6)
  1329. ax.set_title(METRICS_DISPLAY[m], fontsize=9)
  1330. ax.set_xticks(range(1, self.num_treatments + 1))
  1331. ax.set_xticklabels(
  1332. [n.replace('_', ' ') for n in self.treatment_names],
  1333. rotation=45, ha='right', fontsize=7)
  1334. ax.set_aspect('auto')
  1336. # Significance bracket for 2-group comparisons
  1337. if (len(self.stats_results) > m and
  1338. self.stats_results[m]['p_value'] < 0.05 and
  1339. self.num_treatments == 2):
  1340. self._add_significance_bracket(ax, 1, 2, self.stats_results[m]['p_value'])
  1342. fig.suptitle(f"{self.experiment_name.replace('_', ' ')} — Core Metrics",
  1343. fontsize=13, fontweight='bold')
  1344. fig.tight_layout(rect=[0, 0, 1, 0.95])
  1345. self._save_figure(fig, 'Fig2_Metrics')
  1346. return fig
  1348. # ----- FIGURE 3: Food Type Analysis -----
  1349. def _fig3_food_type(self) -> plt.Figure:
  1350. """Generate Figure 3: Chow vs HFD cumulative, preference ratio, bars."""
  1351. fig, axes = plt.subplots(1, 3, figsize=(12, 5))
  1353. # --- 3A: Chow vs HFD time course ---
  1354. ax = axes[0]
  1355. for i, grp in enumerate(self.groups):
  1356. # Chow cumulative
  1357. cum_c = np.cumsum(grp.intake_chow, axis=1)
  1358. mu_c = np.mean(cum_c, axis=0)
  1359. x_h = np.arange(1, cum_c.shape[1] + 1) * self.cfg.bin_size / 3600.0
  1360. ax.plot(x_h, mu_c, '-', color=grp.color, linewidth=1.5,
  1361. label=f"{grp.name.replace('_', ' ')} Chow")
  1362. # HFD cumulative
  1363. cum_h = np.cumsum(grp.intake_hfd, axis=1)
  1364. mu_h = np.mean(cum_h, axis=0)
  1365. ax.plot(x_h, mu_h, '--', color=grp.color, linewidth=1.5,
  1366. label=f"{grp.name.replace('_', ' ')} HFD")
  1368. ax.set_xlabel('Time (h)')
  1369. ax.set_ylabel('Cum. Feeding Time (s)')
  1370. ax.set_title('Cumulative: Chow vs HFD', fontsize=10)
  1371. ax.legend(fontsize=7, frameon=False, loc='best')
  1372. ax.spines['top'].set_visible(False)
  1373. ax.spines['right'].set_visible(False)
  1375. # --- 3B: Preference Ratio over time ---
  1376. ax = axes[1]
  1377. for i, grp in enumerate(self.groups):
  1378. cum_c = np.cumsum(grp.intake_chow, axis=1)
  1379. cum_h = np.cumsum(grp.intake_hfd, axis=1)
  1380. total = cum_c + cum_h
  1381. pref = cum_h / np.maximum(total, 1) # Avoid division by zero
  1382. mu_p = np.nanmean(pref, axis=0)
  1383. err_p = np.nanstd(pref, axis=0, ddof=0) / np.sqrt(max(pref.shape[0], 1))
  1384. x_h = np.arange(1, pref.shape[1] + 1) * self.cfg.bin_size / 3600.0
  1385. self._shaded_error_bar(ax, x_h, mu_p, err_p, grp.color,
  1386. label=grp.name.replace('_', ' '))
  1388. ax.axhline(0.5, color='k', linestyle='--', linewidth=0.8)
  1389. ax.set_xlabel('Time (h)')
  1390. ax.set_ylabel('HFD Preference Ratio')
  1391. ax.set_title('Preference Ratio (HFD / Total)', fontsize=10)
  1392. ax.set_ylim(0, 1)
  1393. ax.legend(fontsize=7, frameon=False, loc='best')
  1394. ax.spines['top'].set_visible(False)
  1395. ax.spines['right'].set_visible(False)
  1397. # --- 3C: Chow vs HFD bar chart ---
  1398. ax = axes[2]
  1399. bar_width = 0.35
  1400. for i, grp in enumerate(self.groups):
  1401. chow_total = np.sum(grp.intake_chow, axis=1)
  1402. hfd_total = np.sum(grp.intake_hfd, axis=1)
  1403. x_pos = i + 1
  1404. ax.bar(x_pos - bar_width / 2, np.mean(chow_total), bar_width,
  1405. color=FOOD_COLORS['chow'], edgecolor='none', alpha=0.7)
  1406. ax.bar(x_pos + bar_width / 2, np.mean(hfd_total), bar_width,
  1407. color=FOOD_COLORS['hfd'], edgecolor='none', alpha=0.7)
  1408. ax.errorbar(x_pos - bar_width / 2, np.mean(chow_total),
  1409. yerr=np.std(chow_total, ddof=0) / np.sqrt(max(len(chow_total), 1)),
  1410. fmt='none', ecolor='k', linewidth=1)
  1411. ax.errorbar(x_pos + bar_width / 2, np.mean(hfd_total),
  1412. yerr=np.std(hfd_total, ddof=0) / np.sqrt(max(len(hfd_total), 1)),
  1413. fmt='none', ecolor='k', linewidth=1)
  1415. ax.set_xticks(range(1, self.num_treatments + 1))
  1416. ax.set_xticklabels([n.replace('_', ' ') for n in self.treatment_names],
  1417. rotation=45, ha='right', fontsize=9)
  1418. ax.set_ylabel('Total Feeding Time (s)')
  1419. ax.set_title('Chow vs HFD by Group', fontsize=10)
  1420. ax.legend(['Chow', 'HFD'], fontsize=8, frameon=False, loc='best')
  1422. fig.suptitle(f"{self.experiment_name.replace('_', ' ')} — Food Type Analysis",
  1423. fontsize=13, fontweight='bold')
  1424. fig.tight_layout(rect=[0, 0, 1, 0.95])
  1425. self._save_figure(fig, 'Fig3_FoodType')
  1426. return fig
  1428. # ----- FIGURE 4: Microstructure -----
  1429. def _fig4_microstructure(self) -> plt.Figure:
  1430. """Generate Figure 4: Bout/IBI distributions, meal metrics, activity heatmap."""
  1431. fig = plt.figure(figsize=(14, 7))
  1432. gs = fig.add_gridspec(2, 3, hspace=0.4, wspace=0.4)
  1434. # --- 4A: Bout duration distribution (boxplot) ---
  1435. ax = fig.add_subplot(gs[0, 0])
  1436. all_bout_dur = []
  1437. group_label_bout = []
  1438. for i, grp in enumerate(self.groups):
  1439. for s in range(grp.n_subs):
  1440. bd = grp.bout_durations[s]
  1441. if bd is not None and len(bd) > 0:
  1442. all_bout_dur.extend(bd)
  1443. group_label_bout.extend([i + 1] * len(bd))
  1445. if all_bout_dur:
  1446. bp_data = [[] for _ in range(self.num_treatments)]
  1447. for val, lbl in zip(all_bout_dur, group_label_bout):
  1448. bp_data[lbl - 1].append(val)
  1449. bp = ax.boxplot(bp_data, patch_artist=True, sym='.',
  1450. labels=[n.replace('_', ' ') for n in self.treatment_names])
  1451. for patch, color in zip(bp['boxes'], self.col):
  1452. patch.set_facecolor((*color, 0.5))
  1453. ax.set_ylabel('Bout Duration (s)')
  1454. ax.set_title('Bout Duration Distribution', fontsize=10)
  1455. ax.tick_params(axis='x', rotation=45)
  1457. # --- 4B: IBI distribution ---
  1458. ax = fig.add_subplot(gs[0, 1])
  1459. all_ibi = []
  1460. group_label_ibi = []
  1461. for i, grp in enumerate(self.groups):
  1462. for s in range(grp.n_subs):
  1463. ib = grp.ibi_values[s]
  1464. if ib is not None and len(ib) > 0:
  1465. all_ibi.extend(ib)
  1466. group_label_ibi.extend([i + 1] * len(ib))
  1468. if all_ibi:
  1469. bp_data = [[] for _ in range(self.num_treatments)]
  1470. for val, lbl in zip(all_ibi, group_label_ibi):
  1471. bp_data[lbl - 1].append(val)
  1472. bp = ax.boxplot(bp_data, patch_artist=True, sym='.',
  1473. labels=[n.replace('_', ' ') for n in self.treatment_names])
  1474. for patch, color in zip(bp['boxes'], self.col):
  1475. patch.set_facecolor((*color, 0.5))
  1476. ax.set_ylabel('IBI (s)')
  1477. ax.set_title('IBI Distribution', fontsize=10)
  1478. ax.tick_params(axis='x', rotation=45)
  1480. # --- 4C: Meal metrics bars (Num_Meals and Meal_Size) ---
  1481. meal_metrics_idx = [6, 7] # 0-based
  1482. for m_i in range(2):
  1483. ax = fig.add_subplot(gs[0, 2]) if m_i == 0 else fig.add_subplot(gs[1, 0])
  1484. # Fix: use separate subplot for each
  1485. if m_i == 1:
  1486. ax = fig.add_subplot(gs[1, 0])
  1487. m = meal_metrics_idx[m_i]
  1488. for ii in range(self.num_treatments):
  1489. vals = self.all_metrics[ii][:, m]
  1490. vals_clean = vals[~np.isnan(vals)]
  1491. if len(vals_clean) == 0:
  1492. continue
  1493. ax.bar(ii + 1, np.mean(vals_clean), color=self.col[ii],
  1494. edgecolor='none', alpha=0.6)
  1495. ax.errorbar(ii + 1, np.mean(vals_clean),
  1496. yerr=np.std(vals_clean, ddof=0) / np.sqrt(len(vals_clean)),
  1497. fmt='none', ecolor='k', linewidth=1)
  1498. jitter = (np.random.rand(len(vals_clean)) - 0.5) * 0.3
  1499. ax.plot(ii + 1 + jitter, vals_clean, '.', color='0.3', markersize=8)
  1500. ax.set_title(METRICS_DISPLAY[m], fontsize=9)
  1501. ax.set_xticks(range(1, self.num_treatments + 1))
  1502. ax.set_xticklabels(
  1503. [n.replace('_', ' ') for n in self.treatment_names],
  1504. rotation=45, ha='right', fontsize=8)
  1506. # --- 4D: Activity heatmap ---
  1507. ax = fig.add_subplot(gs[1, 1:3])
  1508. all_heatmap_data = []
  1509. heatmap_labels = []
  1510. for i, grp in enumerate(self.groups):
  1511. raw = grp.intake
  1512. for s in range(raw.shape[0]):
  1513. all_heatmap_data.append(raw[s, :])
  1514. sname = os.path.splitext(grp.subjects[s])[0]
  1515. heatmap_labels.append(
  1516. f"{grp.name.replace('_', '')}_{sname}")
  1518. if all_heatmap_data:
  1519. heatmap_arr = np.array(all_heatmap_data)
  1520. im = ax.imshow(heatmap_arr, aspect='auto', cmap='hot',
  1521. interpolation='nearest')
  1522. fig.colorbar(im, ax=ax, fraction=0.02, pad=0.04)
  1523. ax.set_xlabel(f'Time bins ({self.cfg.bin_size} s)')
  1524. ax.set_ylabel('Subject')
  1525. ax.set_title('Feeding Activity Heatmap', fontsize=10)
  1526. ax.set_yticks(range(len(heatmap_labels)))
  1527. ax.set_yticklabels(heatmap_labels, fontsize=6)
  1529. # Add group separators
  1530. cum_n = 0
  1531. for i in range(self.num_treatments - 1):
  1532. cum_n += self.groups[i].n_subs
  1533. ax.axhline(cum_n - 0.5, color='white', linewidth=2)
  1535. fig.suptitle(f"{self.experiment_name.replace('_', ' ')} — Microstructure Analysis",
  1536. fontsize=13, fontweight='bold')
  1537. fig.tight_layout(rect=[0, 0, 1, 0.95])
  1538. self._save_figure(fig, 'Fig4_Microstructure')
  1539. return fig
  1541. # ----- FIGURE 5: Advanced Analytics -----
  1542. def _fig5_advanced(self) -> plt.Figure:
  1543. """Generate Figure 5: Satiety ratio, bout duration vs time, curve fitting."""
  1544. fig, axes = plt.subplots(1, 3, figsize=(12, 5))
  1546. # --- 5A: Satiety Ratio bars ---
  1547. ax = axes[0]
  1548. m = 9 # Satiety_Ratio (0-based)
  1549. for i in range(self.num_treatments):
  1550. vals = self.all_metrics[i][:, m]
  1551. vals_clean = vals[~np.isnan(vals)]
  1552. if len(vals_clean) == 0:
  1553. continue
  1554. ax.bar(i + 1, np.mean(vals_clean), color=self.col[i],
  1555. edgecolor='none', alpha=0.6)
  1556. ax.errorbar(i + 1, np.mean(vals_clean),
  1557. yerr=np.std(vals_clean, ddof=0) / np.sqrt(len(vals_clean)),
  1558. fmt='none', ecolor='k', linewidth=1)
  1559. jitter = (np.random.rand(len(vals_clean)) - 0.5) * 0.3
  1560. ax.plot(i + 1 + jitter, vals_clean, '.', color='0.3', markersize=8)
  1561. ax.set_title('Satiety Ratio', fontsize=10)
  1562. ax.set_ylabel('IBI_post / Bout Duration')
  1563. ax.set_xticks(range(1, self.num_treatments + 1))
  1564. ax.set_xticklabels([n.replace('_', ' ') for n in self.treatment_names],
  1565. rotation=45, ha='right', fontsize=9)
  1567. # --- 5B: Bout duration vs time (scatter) ---
  1568. ax = axes[1]
  1569. for i, grp in enumerate(self.groups):
  1570. all_starts = []
  1571. all_durs = []
  1572. for s in range(grp.n_subs):
  1573. dm = grp.DM[s]
  1574. if dm is None:
  1575. continue
  1576. feed = dm[np.isin(dm[:, self.cfg.col_code], [1, 2])]
  1577. if len(feed) > 0:
  1578. all_starts.extend(feed[:, 0])
  1579. all_durs.extend(feed[:, 2])
  1580. if all_starts:
  1581. ax.scatter(np.array(all_starts) / 3600.0, all_durs,
  1582. s=15, c=[grp.color], alpha=0.4,
  1583. label=grp.name.replace('_', ' '))
  1585. ax.set_xlabel('Time (h)')
  1586. ax.set_ylabel('Bout Duration (s)')
  1587. ax.set_title('Bout Duration Over Time', fontsize=10)
  1588. ax.set_xlim(0, self.cfg.exp_duration / 3600.0)
  1589. ax.legend(fontsize=7, frameon=False, loc='best')
  1590. ax.spines['top'].set_visible(False)
  1591. ax.spines['right'].set_visible(False)
  1593. # --- 5C: Cumulative curve fitting ---
  1594. ax = axes[2]
  1595. for i, grp in enumerate(self.groups):
  1596. raw = grp.intake
  1597. if raw is None or raw.size == 0:
  1598. continue
  1599. cum_mean = np.mean(np.cumsum(raw, axis=1), axis=0)
  1600. x_sec = np.arange(1, len(cum_mean) + 1) * self.cfg.bin_size
  1601. x_h = x_sec / 3600.0
  1602. ax.plot(x_h, cum_mean, 'o', color=grp.color, markersize=3,
  1603. label=grp.name.replace('_', ' '))
  1605. # Fit exponential saturation: y = a * (1 - exp(-b * x))
  1606. try:
  1607. def sat_func(x, a, b):
  1608. return a * (1 - np.exp(-b * x))
  1610. popt, _ = curve_fit(sat_func, x_sec, cum_mean,
  1611. p0=[max(cum_mean) * 1.2, 1 / 1800.0],
  1612. bounds=([0, 0], [np.inf, np.inf]),
  1613. maxfev=5000)
  1614. x_fine = np.linspace(0, self.cfg.exp_duration, 200)
  1615. ax.plot(x_fine / 3600.0, sat_func(x_fine, *popt), '-',
  1616. color=grp.color, linewidth=1.5)
  1617. except Exception:
  1618. pass
  1620. ax.set_xlabel('Time (h)')
  1621. ax.set_ylabel('Cum. Feeding Time (s)')
  1622. ax.set_title('Cumulative Curve Fit', fontsize=10)
  1623. ax.legend(fontsize=7, frameon=False, loc='best')
  1624. ax.spines['top'].set_visible(False)
  1625. ax.spines['right'].set_visible(False)
  1627. fig.suptitle(f"{self.experiment_name.replace('_', ' ')} — Advanced Analytics",
  1628. fontsize=13, fontweight='bold')
  1629. fig.tight_layout(rect=[0, 0, 1, 0.95])
  1630. self._save_figure(fig, 'Fig5_Advanced')
  1631. return fig
  1633. # =====================================================================
  1634. # EXCEL EXPORT METHODS
  1635. # =====================================================================
  1637. def _export_summary_excel(self, xls_path: str) -> None:
  1638. """Export Results Summary Excel (4 sheets)."""
  1639. with pd.ExcelWriter(xls_path, engine='openpyxl') as writer:
  1641. # --- Sheet 1: Individual Data ---
  1642. rows = []
  1643. for i, grp in enumerate(self.groups):
  1644. for s in range(grp.n_subs):
  1645. sname = os.path.splitext(grp.subjects[s])[0]
  1646. row = {'Subject': sname, 'Treatment': grp.name}
  1647. for m in range(N_METRICS):
  1648. row[METRICS_LABELS[m]] = self.all_metrics[i][s, m]
  1649. rows.append(row)
  1650. pd.DataFrame(rows).to_excel(writer, sheet_name='Individual_Data', index=False)
  1652. # --- Sheet 2: Group Stats ---
  1653. rows = []
  1654. for i in range(self.num_treatments):
  1655. row = {'Treatment': self.treatment_names[i]}
  1656. for m in range(N_METRICS):
  1657. vals = self.all_metrics[i][:, m]
  1658. vals_clean = vals[~np.isnan(vals)]
  1659. if len(vals_clean) > 0:
  1660. row[f'{METRICS_LABELS[m]}_Mean'] = np.mean(vals_clean)
  1661. row[f'{METRICS_LABELS[m]}_SEM'] = (
  1662. np.std(vals_clean, ddof=0) / np.sqrt(len(vals_clean)))
  1663. else:
  1664. row[f'{METRICS_LABELS[m]}_Mean'] = np.nan
  1665. row[f'{METRICS_LABELS[m]}_SEM'] = np.nan
  1666. rows.append(row)
  1667. pd.DataFrame(rows).to_excel(writer, sheet_name='Group_Stats', index=False)
  1669. # --- Sheet 3: Statistical Tests ---
  1670. rows = []
  1671. for m in range(N_METRICS):
  1672. sr = self.stats_results[m]
  1673. rows.append({
  1674. 'Metric': METRICS_DISPLAY[m],
  1675. 'Test': sr['test_name'],
  1676. 'Statistic': sr['statistic'],
  1677. 'P_Value': sr['p_value'],
  1678. 'Decision': sr['decision'],
  1679. 'Note': sr['note'],
  1680. })
  1681. pd.DataFrame(rows).to_excel(writer, sheet_name='Statistical_Tests', index=False)
  1683. # --- Sheet 4: CV Table ---
  1684. rows = []
  1685. for i in range(self.num_treatments):
  1686. row = {'Treatment': self.treatment_names[i]}
  1687. for m in range(N_METRICS):
  1688. row[f'{METRICS_LABELS[m]}_CV_pct'] = self.cv_table[i, m]
  1689. rows.append(row)
  1690. pd.DataFrame(rows).to_excel(writer, sheet_name='Variability_CV', index=False)
  1692. def _export_deep_excel(self, xls_path: str) -> None:
  1693. """Export Deep Analysis Data Excel (5 sheets)."""
  1694. with pd.ExcelWriter(xls_path, engine='openpyxl') as writer:
  1696. # --- Sheet 1: Raw Events ---
  1697. rows = []
  1698. for i, grp in enumerate(self.groups):
  1699. for s in range(grp.n_subs):
  1700. data = grp.DM[s]
  1701. if data is None or data.shape[1] <= self.cfg.col_code:
  1702. continue
  1703. codes = data[:, self.cfg.col_code]
  1704. valid = np.isin(codes, [1, 2, 3])
  1705. if not np.any(valid):
  1706. continue
  1707. evs = data[valid]
  1708. sname = os.path.splitext(grp.subjects[s])[0]
  1709. for e in range(len(evs)):
  1710. rows.append({
  1711. 'Subject': sname,
  1712. 'Treatment': grp.name,
  1713. 'Start_s': evs[e, 0],
  1714. 'End_s': evs[e, 1],
  1715. 'Duration_s': evs[e, 2],
  1716. 'Power_dB': evs[e, 4],
  1717. 'Type_Code': int(evs[e, self.cfg.col_code]),
  1718. })
  1719. pd.DataFrame(rows).to_excel(writer, sheet_name='Events_Raw', index=False)
  1721. # --- Sheet 2: Cumulative Time Course ---
  1722. timeline_min = np.arange(1, self.n_bins + 1) * self.cfg.bin_size / 60.0
  1723. cum_dict = {'Time_Min': timeline_min}
  1724. for i, grp in enumerate(self.groups):
  1725. cum_d = np.cumsum(grp.intake, axis=1)
  1726. for s in range(cum_d.shape[0]):
  1727. sname = os.path.splitext(grp.subjects[s])[0]
  1728. col_name = f"{re.sub('[^a-zA-Z0-9]', '', grp.name)}_{sname}"
  1729. cum_dict[col_name] = cum_d[s, :]
  1730. pd.DataFrame(cum_dict).to_excel(
  1731. writer, sheet_name='TimeCourse_Cumulative', index=False)
  1733. # --- Sheet 3: Feeding Rate ---
  1734. rate_dict = {'Time_Min': timeline_min}
  1735. for i, grp in enumerate(self.groups):
  1736. for s in range(grp.intake.shape[0]):
  1737. sname = os.path.splitext(grp.subjects[s])[0]
  1738. col_name = f"{re.sub('[^a-zA-Z0-9]', '', grp.name)}_{sname}"
  1739. rate_dict[col_name] = grp.intake[s, :]
  1740. pd.DataFrame(rate_dict).to_excel(
  1741. writer, sheet_name='TimeCourse_FeedingRate', index=False)
  1743. # --- Sheet 4: IBI Time Course ---
  1744. n_bins_ibi = int(self.cfg.exp_duration // self.cfg.ibi_window)
  1745. time_ibi = (np.arange(n_bins_ibi) * self.cfg.ibi_window +
  1746. self.cfg.ibi_window / 2.0) / 60.0
  1747. ibi_dict = {'Time_Min': time_ibi}
  1748. for i, grp in enumerate(self.groups):
  1749. raw_ibi = self.c_ibi_all[:, :, i]
  1750. for s in range(grp.n_subs):
  1751. subj_ibi_trace = np.full(n_bins_ibi, np.nan)
  1752. for b in range(n_bins_ibi):
  1753. idx_s = b * self.cfg.ibi_window
  1754. idx_e = (b + 1) * self.cfg.ibi_window
  1755. vals = raw_ibi[s, idx_s:idx_e]
  1756. vals = vals[~np.isnan(vals)]
  1757. if len(vals) > 0:
  1758. subj_ibi_trace[b] = np.mean(vals)
  1759. sname = os.path.splitext(grp.subjects[s])[0]
  1760. col_name = f"{re.sub('[^a-zA-Z0-9]', '', grp.name)}_{sname}"
  1761. ibi_dict[col_name] = subj_ibi_trace
  1762. pd.DataFrame(ibi_dict).to_excel(
  1763. writer, sheet_name='TimeCourse_IBI', index=False)
  1765. # --- Sheet 5: Meal Data ---
  1766. rows = []
  1767. for i, grp in enumerate(self.groups):
  1768. for s in range(grp.n_subs):
  1769. meals = grp.meals[s]
  1770. if meals is None or len(meals) == 0:
  1771. continue
  1772. sname = os.path.splitext(grp.subjects[s])[0]
  1773. for meal_num in range(len(meals)):
  1774. rows.append({
  1775. 'Subject': sname,
  1776. 'Treatment': grp.name,
  1777. 'Meal_Num': meal_num + 1,
  1778. 'Meal_Start_s': meals[meal_num, 0],
  1779. 'Meal_End_s': meals[meal_num, 1],
  1780. 'Meal_Duration_s': meals[meal_num, 2],
  1781. 'Num_Bouts_In_Meal': int(meals[meal_num, 3]),
  1782. })
  1783. pd.DataFrame(rows).to_excel(writer, sheet_name='Meal_Clusters', index=False)

group_analysis.py at commit 63b05d6, under MIT · at the source

Overview

Authors: Elvi Gil Lievana1,2, Benjamin Arroyo1,2, Jesús Pérez-Ortega1,2,3, Axel Lopez1,2, Luis Rodriguez-Blanco1,2, Xarenny Diaz1,2, Gustavo Hernandez1,2, Alam Coss1,3, Emily Alway4,5,6, Naama Reicher4,5, Enrique Hernández-Lemus7,8, Maya Kaelberer4,9, Diego V Bohórquez4,5,6,10,11,12,13, Ranier Gutierrez1,2
ORCID iDs: Jesús Pérez-Ortega, Emily Alway, Ranier Gutierrez
13 affiliations
  1. Laboratory Neurobiology of Appetite; Department of Pharmacology, CINVESTAV, Mexico City, Mexico
  2. Laboratory Neurobiology of Appetite; Center for Research on Aging (CIE), CINVESTAV, Mexico City, Mexico
  3. Facultad de Ingeniería, Universidad Nacional Autónoma de México, Mexico City, Mexico
  4. Department of Medicine, Duke University, Durham, United States
  5. Laboratory of Gut-Brain Neurobiology, Duke University, Durham, United States
  6. Department of Neurobiology, Duke University, Durham, United States
  7. Computational Genomics Division, National Institute of Genomic Medicine (INMEGEN), Mexico City, Mexico
  8. Center for Complexity Sciences, Universidad Nacional Autónoma de México, Mexico City, Mexico
  9. Department of Physiology, University of Arizona, Tucson, United States
  10. Department of Molecular Genetics and Microbiology, Duke University, Durham, United States
  11. Department of Pathology, Duke University, Durham, United States
  12. Department of Cell Biology, Duke University, Durham, United States
  13. Duke Institute for Brain Sciences, Duke University, Durham, United States
Institutions: Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional (Mexico); Universidad Nacional Autónoma de México (Mexico); Duke University (United States); National Institute of Genomic Medicine (Mexico); University of Arizona (United States)
Journal: eLife, volume 14, article RP108663
Dates: published online 14 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.108663 · PMID 42444443 · PMCID PMC13368178 · OpenAlex W4415477351
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), methods / tools (subfield)
Methods: Spectral & time-frequency, Statistics, Machine learning, Preprocessing, fMRI & imaging, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions
Keywords: Mouse
MeSH: Acoustics*, Feeding Behavior*, Animals, Eating, Male, Mice, Mice, Inbred C57BL, Neurons, Semaglutide (* major topic)
Topic: Infant Health and Development (Pharmacy, Health Professions), according to OpenAlex
Funding: Secretaría de Ciencia, Humanidades, Tecnología e Innovación (CBF-2026-1938, CF-2023-G-518); NIDDK NIH HHS (F30 DK136229, F32 DK139628, R01 DK131112, R03 DK114500, R01 DK132070, K01 DK131403); NIMH NIH HHS (DP2 MH122402); NCCIH NIH HHS (R21 AT010818)
Citations: cited by 2 papers (Europe PMC); 68 references in the paper

Abstract

Elucidating the neuronal circuits that govern appetite requires precise, high-resolution monitoring of the microstructure of solid food consumption, a need unmet by existing tools, which are either costly or lack the temporal resolution to align feeding events with neuronal activity. To overcome this, we developed the Crunchometer, a low-cost, open-source acoustic system that uses computational algorithms to generate high-resolution feeding ethograms from the sounds produced during solid food consumption. Validation across energy states (hunger/satiety) confirmed its sensitivity to changes in feeding microstructure, and the system reliably detected semaglutide-induced suppression of intake and reduced preference for a high-fat diet. Leveraging its seamless integration with in vivo recordings in freely behaving mice, we paired the Crunchometer with lateral hypothalamus (LH) electrophysiology to identify ‘meal-related’ neurons that track entire meals rather than individual bouts. Calcium imaging further revealed that distinct subsets of LH GABAergic and glutamatergic neurons were tuned to feeding only, to licking only, or to both behaviors. Thus, LH neuronal ensembles differentially encode the consumption of solid food versus liquid sucrose. These findings demonstrate that the Crunchometer is a robust, accessible platform for dissecting the neural correlates of feeding behavior at the resolution of a single bite.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 18 matches between paragraphs and lines of code.

OSF bmkdc

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 6 files
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
At the source:

RanierLabNeurobiologyAppetite2026/TheCrunchometerPy

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 63b05d6d0066b7277f2ef819fa61a238f406c206, 2 March 2026
Languages: Python (44), MATLAB (2), Shell (1)
Size: 87 files, 47 scripts
Software Heritage: not archived
Found in: the references
Holds: README, license file, environment (crunchometer_python/pyproject.toml, crunchometer_python/requirements.txt), continuous integration
Not found: CITATION.cff, tests, documentation
Tools: NumPy (32 files), SciPy (13 files), OpenCV (10 files), Matplotlib (7 files), pandas (7 files), h5py (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
49 files
At the source: github.com/RanierLabNeurobiologyAppetite2026/TheCrunchometerPy

Zenodo 8422691

License: other-open
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
11 files
At the source:

Zenodo 8423311

License: other-open
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (15 files), Image Processing Toolbox (4 files), Signal Processing Toolbox (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
113 files
At the source:

perezortegaj/moussionenergy

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 51a2154f8e1b06dec210dad00f3eb200acaaaad4, 9 October 2023
Languages: MATLAB (10)
Size: 13 files, 10 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: license file
Not found: README, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
11 files
At the source: github.com/PerezOrtegaJ/MoussionEnergy/tree/v1.0.0

perezortegaj/xsembles2p

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 274cdb9bce3c6eec0db909c9a40dd6a2f30cdc01, 24 April 2024
Languages: MATLAB (112)
Size: 122 files, 112 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, license file, CITATION.cff, tests
Not found: environment file, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (15 files), Image Processing Toolbox (4 files), Signal Processing Toolbox (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
114 files
At the source: github.com/PerezOrtegaJ/Xsembles2P/tree/v1.1.0

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 6 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 290 scripts, each with its path and the digest of its content;
  • 18 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Benchmark and Crunchometer software is available on OSF (Arroyo et al., 2025) at https://osf.io/bmkdc/. The OSF repository contains two folders under Files. The TheCrunchometerV2_Matlab folder includes both the classic scripts and an App GUI version, compatible with Windows, macOS (Intel and Apple Silicon), and Linux. This version requires a MATLAB license and FFmpeg, and is limited to videos up to 4 hours long. The TheCrunchometerV2_Python folder contains a standalone macOS Silicon executable (TheCrunchometer-2.0.0-macOS-arm64.dmg); this is the recommended version. It requires no MATLAB license or additional dependencies, supports 24-hour video files, and runs optimally on Apple Silicon hardware (e.g., Mac mini M4). To install, double-click the.dmg file and authorize the application via Apple menu → System Settings → Privacy & Security → Open Anyway. A standalone Windows executable is also available; it handles 24 hour videos but runs more slowly than the macOS build. Alternatively, a Python source version of the Crunchometer is available on GitHub (https://github.com/RanierLabNeurobiologyAppetite2026/TheCrunchometerPy, Arroyo et al., 2026). We recommend directly downloading the standalone software from GitHub using the following link: for MacOS Silicon: TheCrunchometer-2.0.0-macOS-arm64.dmg (https://github.com/RanierLabNeurobiologyAppetite2026/TheCrunchometerPy/releases/download/v2.0.0/TheCrunchometer-2.0.0-macOS-arm64.dmg) (to install, double-click the.dmg file and authorize the application via Apple menu → System Settings → Privacy & Security → Open Anyway); for Windows: Setup_TheCrunchometer-2.0.0-Windows-x64.exe (https://github.com/RanierLabNeurobiologyAppetite2026/TheCrunchometerPy/releases/download/v2.1.0/Setup_TheCrunchometer-2.0.0-Windows-x64.exe). An example folder, CrunchometerRunExample, is also provided, containing Audio_Data (.ogg) and Video_Data (.mkv) subfolders from a 10 minute recording for testing the Crunchometer software. A tutorial video on installing and running the Crunchometer is available at https://youtu.be/cubBnfOoGy8.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 14 authors, 1 keyword, 9 MeSH terms, 4 funders, 66 references.

Cite

This paper

Gil Lievana, E., Arroyo, B., Pérez-Ortega, J., Lopez, A., Rodriguez-Blanco, L., Diaz, X., Hernandez, G., Coss, A., Alway, E., Reicher, N., Hernández-Lemus, E., Kaelberer, M., Bohórquez, D. V., & Gutierrez, R. (2026). The Crunchometer, a low-cost, open-source acoustic analysis of feeding microstructure. eLife, 14, RP108663. https://doi.org/10.7554/elife.108663

BibTeX

@article{gillievana2026crunchometer,
author = {Gil Lievana, Elvi and Arroyo, Benjamin and Pérez-Ortega, Jesús and Lopez, Axel and Rodriguez-Blanco, Luis and Diaz, Xarenny and Hernandez, Gustavo and Coss, Alam and Alway, Emily and Reicher, Naama and Hernández-Lemus, Enrique and Kaelberer, Maya and Bohórquez, Diego V and Gutierrez, Ranier},
title = {{The Crunchometer, a low-cost, open-source acoustic analysis of feeding microstructure}},
journal = {eLife},
year = {2026},
month = jul,
volume = {14},
pages = {RP108663},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.108663},
url = {https://doi.org/10.7554/elife.108663},
pmid = {42444443},
pmcid = {PMC13368178}
}

RIS

TY - JOUR
AU - Gil Lievana, Elvi
AU - Arroyo, Benjamin
AU - Pérez-Ortega, Jesús
AU - Lopez, Axel
AU - Rodriguez-Blanco, Luis
AU - Diaz, Xarenny
AU - Hernandez, Gustavo
AU - Coss, Alam
AU - Alway, Emily
AU - Reicher, Naama
AU - Hernández-Lemus, Enrique
AU - Kaelberer, Maya
AU - Bohórquez, Diego V
AU - Gutierrez, Ranier
TI - The Crunchometer, a low-cost, open-source acoustic analysis of feeding microstructure
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/07/14
VL - 14
SP - RP108663
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.108663
UR - https://doi.org/10.7554/elife.108663
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.108663",
"type": "article-journal",
"title": "The Crunchometer, a low-cost, open-source acoustic analysis of feeding microstructure",
"container-title": "eLife",
"author": [
{
"family": "Gil Lievana",
"given": "Elvi"
},
{
"family": "Arroyo",
"given": "Benjamin"
},
{
"family": "Pérez-Ortega",
"given": "Jesús"
},
{
"family": "Lopez",
"given": "Axel"
},
{
"family": "Rodriguez-Blanco",
"given": "Luis"
},
{
"family": "Diaz",
"given": "Xarenny"
},
{
"family": "Hernandez",
"given": "Gustavo"
},
{
"family": "Coss",
"given": "Alam"
},
{
"family": "Alway",
"given": "Emily"
},
{
"family": "Reicher",
"given": "Naama"
},
{
"family": "Hernández-Lemus",
"given": "Enrique"
},
{
"family": "Kaelberer",
"given": "Maya"
},
{
"family": "Bohórquez",
"given": "Diego V"
},
{
"family": "Gutierrez",
"given": "Ranier"
}
],
"container-title-short": "Elife",
"volume": "14",
"page": "RP108663",
"DOI": "10.7554/elife.108663",
"PMID": "42444443",
"PMCID": "PMC13368178",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.108663",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
14
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76581-6 [code]
Thalamocortical bursts encode reward contingencies and drive associative learning.
Journal: Nature communications
In common: OpenCV, h5py, Signal Processing Toolbox, 7 other tools, mouse
[2] doi:10.1038/s41592-026-03154-2 [code]
Simultaneous single-cell calcium imaging of neuronal population activity and brain-wide BOLD fMRI.
Journal: Nature methods
In common: OpenCV, h5py, Signal Processing Toolbox, 7 other tools, mouse
[3] doi:10.1038/s41593-026-02376-z [code]
A framework for comparative analysis of human and mouse cortical neuron dendrites in corresponding brain regions.
Journal: Nature neuroscience
In common: OpenCV, h5py, Signal Processing Toolbox, 7 other tools, mouse
[4] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: OpenCV, h5py, Signal Processing Toolbox, 7 other tools, mouse
[5] doi:10.1016/j.isci.2026.117375 [code]
Motor priming is associated with widespread recruitment into neural ensembles and more rapid ensemble transitions.
Journal: iScience
In common: OpenCV, h5py, Signal Processing Toolbox, 7 other tools
[6] doi:10.1038/s41467-026-71458-0 [code]
Early differential impact of MeCP2 mutations on functional networks in Rett syndrome patient-derived human cortical organoids.
Journal: Nature communications
In common: OpenCV, h5py, Signal Processing Toolbox, 6 other tools
[7] doi:10.1162/imag.a.1299 [code]
A modular semantic-structural pipeline for visual decoding from primate spiking data via selective temporal integration.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: h5py, Signal Processing Toolbox, Image Processing Toolbox, 6 other tools, methods / tools
[8] doi:10.3389/fendo.2026.1828487 [code]
Castration-induced nigrostriatal deficits are linked to reduced TrkB and loss of mature spines in the dorsal striatum.
Journal: Frontiers in endocrinology
In common: OpenCV, h5py, Signal Processing Toolbox, 6 other tools, mouse
[9] doi:10.1038/s41593-026-02255-7 [code]
Neural circuits encode prior knowledge of temporal statistics.
Journal: Nature neuroscience
In common: Signal Processing Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 5 other tools, mouse, 1 reference
[10] doi:10.1016/j.neuron.2026.03.034 [code]
Dentate gyrus interneurons modulate winner-take-all network dynamics in freely behaving mice.
Journal: Neuron
In common: h5py, Signal Processing Toolbox, Image Processing Toolbox, 6 other tools, mouse

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.