From 021a0446dcbcce5598d4c7a4ba4185c1ee17d672 Mon Sep 17 00:00:00 2001
From: Stephen Royle <s.j.royle@warwick.ac.uk>
Date: Tue, 16 Sep 2025 08:48:33 +0100
Subject: [PATCH 1/4] v.0.3.12 draft

---
 DESCRIPTION               |   7 +-
 NAMESPACE                 |  17 +--
 NEWS.md                   |   5 +
 R/TrackMateR-package.R    |  17 +--
 R/calculateTrackDensity.R |  67 ++++++------
 R/readTrackMateXML.R      | 223 ++++++++++++++------------------------
 6 files changed, 128 insertions(+), 208 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 0d01569..3f90177 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: TrackMateR
 Type: Package
 Title: Working with TrackMate outputs in R
-Version: 0.3.11
+Version: 0.3.12
 Authors@R:
     person(given = "Stephen J",
            family = "Royle",
@@ -21,16 +21,13 @@ VignetteBuilder:
   knitr,
   rmarkdown
 Imports:
-  doParallel,
   dplyr,
-  foreach,
   ggplot2,
   ggforce,
   patchwork,
-  parallelly,
   reshape2,
   utils,
-  XML,
+  xml2,
   zoo
 Suggests:
   knitr,
diff --git a/NAMESPACE b/NAMESPACE
index 4736767..3519d8b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -30,25 +30,14 @@ export(plot_tm_width)
 export(readGTFile)
 export(readTrackMateXML)
 export(reportDataset)
-import(doParallel)
 import(dplyr)
 import(ggplot2)
-import(parallelly)
 import(patchwork)
-importFrom(XML,getNodeSet)
-importFrom(XML,xmlDoc)
-importFrom(XML,xmlGetAttr)
-importFrom(XML,xmlParse)
-importFrom(XML,xpathApply)
-importFrom(XML,xpathSApply)
-importFrom(foreach,"%do%")
-importFrom(foreach,"%dopar%")
-importFrom(foreach,foreach)
 importFrom(ggforce,geom_sina)
 importFrom(graphics,frame)
 importFrom(graphics,hist)
-importFrom(parallelly,availableCores)
 importFrom(reshape2,melt)
+importFrom(stats,aggregate)
 importFrom(stats,approx)
 importFrom(stats,coef)
 importFrom(stats,dist)
@@ -62,4 +51,8 @@ importFrom(utils,data)
 importFrom(utils,install.packages)
 importFrom(utils,read.csv)
 importFrom(utils,write.csv)
+importFrom(xml2,read_xml)
+importFrom(xml2,xml_attr)
+importFrom(xml2,xml_find_all)
+importFrom(xml2,xml_find_first)
 importFrom(zoo,rollmean)
diff --git a/NEWS.md b/NEWS.md
index 16a6323..8d75350 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,10 @@
 # TrackMateR
 
+# TrackMateR 0.3.12
+
+- Moved from {XML} to {xml2} for reading TrackMate XML files. Future-proof but currently slower.
+- Speed improvements in analysis of large datasets
+
 # TrackMateR 0.3.11
 
 - Better handling of more than one channel in TrackMate XML files
diff --git a/R/TrackMateR-package.R b/R/TrackMateR-package.R
index 114adfc..177c01c 100644
--- a/R/TrackMateR-package.R
+++ b/R/TrackMateR-package.R
@@ -1,11 +1,8 @@
-#' @importFrom foreach %do%
-#' @importFrom foreach %dopar%
-#' @importFrom foreach foreach
 #' @importFrom ggforce geom_sina
 #' @importFrom graphics frame
 #' @importFrom graphics hist
-#' @importFrom parallelly availableCores
 #' @importFrom reshape2 melt
+#' @importFrom stats aggregate
 #' @importFrom stats approx
 #' @importFrom stats quantile
 #' @importFrom stats sd
@@ -19,17 +16,13 @@
 #' @importFrom utils install.packages
 #' @importFrom utils read.csv
 #' @importFrom utils write.csv
-#' @importFrom XML xmlParse
-#' @importFrom XML getNodeSet
-#' @importFrom XML xpathSApply
-#' @importFrom XML xmlGetAttr
-#' @importFrom XML xmlDoc
-#' @importFrom XML xpathApply
+#' @importFrom xml2 read_xml
+#' @importFrom xml2 xml_find_all
+#' @importFrom xml2 xml_find_first
+#' @importFrom xml2 xml_attr
 #' @importFrom zoo rollmean
 #' @import ggplot2
 #' @import dplyr
 #' @import patchwork
-#' @import doParallel
-#' @import parallelly
 NULL
 #> NULL
diff --git a/R/calculateTrackDensity.R b/R/calculateTrackDensity.R
index f6aecf0..d040357 100644
--- a/R/calculateTrackDensity.R
+++ b/R/calculateTrackDensity.R
@@ -16,49 +16,46 @@
 #' tdDF <- calculateTrackDensity(dataList = tmObj, radius = 2)
 #' @export
 
-calculateTrackDensity <- function(dataList, radius = 1) {
-  x <- y <- trace <- NULL
 
-  if(inherits(dataList, "list")) {
-    df <- dataList[[1]]
-    calibration <- dataList[[2]]
-  } else {
+calculateTrackDensity <- function(dataList, radius = 1) {
+  # base R optimized version
+  if (!inherits(dataList, "list")) {
     cat("Function requires a list of TrackMate data and calibration data\n")
     return(NULL)
   }
-  # make a list of all traces (this is the filtered list of traces from TrackMate XML)
+  df <- dataList[[1]]
+  calibration <- dataList[[2]]
   traceList <- unique(df$trace)
 
-  # for each trace, find the first frame
-  for (i in traceList) {
-    a <- df %>%
-      filter(trace == i) %>%
-      select(frame, x, y)
-    frame0 <- a$frame[1]
-    x0 <- a$x[1]
-    y0 <- a$y[1]
-    # select first frame for this track
-    a <- df %>%
-      filter(frame == frame0) %>%
-      select(x, y)
-    # calculate the distance from x0, y0 to all other coords
-    distances <- find_distances(x0,y0,a)
-    # count how many are less than search radius (this will include the track itself, so subtract 1)
-    neighbours <- sum(distances <= radius, na.rm = TRUE) - 1
-    # calculate how much of the search circle was inside the frame
-    search_fraction <- find_td_area(r = radius, xy = c(x0,y0), a= c(0,calibration[3,1]), b = c(0,calibration[4,1])) / (pi * radius^2)
-    subdf <- data.frame(trace = i,
-                        neighbours = neighbours,
-                        fraction = search_fraction)
-    if(i == traceList[1]) {
-      dfall <- subdf
-    } else {
-      dfall <- rbind(dfall,subdf)
-    }
+  # Precompute starting positions for all traces
+  first_frame_idx <- match(traceList, df$trace)
+  frame0s <- df$frame[first_frame_idx]
+  x0s <- df$x[first_frame_idx]
+  y0s <- df$y[first_frame_idx]
+
+  # Preallocate result vectors
+  neighbours <- numeric(length(traceList))
+  fractions <- numeric(length(traceList))
+
+  for (j in seq_along(traceList)) {
+    frame0 <- frame0s[j]
+    x0 <- x0s[j]
+    y0 <- y0s[j]
+    # select all tracks in the same starting frame
+    idx <- which(df$frame == frame0)
+    x_all <- df$x[idx]
+    y_all <- df$y[idx]
+    # vectorized distance calculation
+    dists <- sqrt((x_all - x0)^2 + (y_all - y0)^2)
+    neighbours[j] <- sum(dists <= radius, na.rm = TRUE) - 1
+    # search area fraction
+    fractions[j] <- find_td_area(r = radius, xy = c(x0, y0), a = c(0, calibration[3,1]), b = c(0, calibration[4,1])) / (pi * radius^2)
   }
-  # divide the count by the fraction of circle that was inside the frame to give "density" - do this at the end
-  dfall$density <- dfall$neighbours / dfall$fraction
 
+  dfall <- data.frame(trace = traceList,
+                      neighbours = neighbours,
+                      fraction = fractions)
+  dfall$density <- dfall$neighbours / dfall$fraction
   return(dfall)
 }
 
diff --git a/R/readTrackMateXML.R b/R/readTrackMateXML.R
index 48239fc..f5281af 100644
--- a/R/readTrackMateXML.R
+++ b/R/readTrackMateXML.R
@@ -15,147 +15,99 @@
 #' calibrationDF <- tmObj[[2]]
 #' @export
 
-readTrackMateXML<- function(XMLpath, slim = FALSE){
-  # check if XML file exists
-  if(!file.exists(XMLpath)) {
+
+readTrackMateXML <- function(XMLpath, slim = FALSE) {
+  # Requires xml2 package
+  if (!file.exists(XMLpath)) {
     cat("XML file does not exist: ", XMLpath, "\n")
     return(NULL)
   }
-  # get necessary XMLNodeSet
-  e <- xmlParse(XMLpath)
-  track <- getNodeSet(e, "//Track")
-  if(length(track) == 0) {
+
+  e <- read_xml(XMLpath)
+  track_nodes <- xml_find_all(e, ".//Track")
+  if (length(track_nodes) == 0) {
     cat("No tracks found in XML file\n")
     return(NULL)
   }
-  filtered <- getNodeSet(e, "//TrackID")
-  subdoc <- getNodeSet(e,"//AllSpots//SpotsInFrame//Spot")
-  sublist <- getNodeSet(e,"//FeatureDeclarations//SpotFeatures//Feature/@feature")
-  attrName <- c("name",unlist(sublist))
-  # what are the units?
-  attrList <- c("spatialunits","timeunits")
-  unitVec <- sapply(attrList, function(x) xpathSApply(e, "//Model", xmlGetAttr, x))
-  unitVec <- c(unitVec,"widthpixels","heightpixels","ntraces","maxframes")
-  attrList <- c("pixelwidth","timeinterval","width","height")
-  valueVec <- sapply(attrList, function(x) xpathSApply(e, "//ImageData", xmlGetAttr, x))
-  calibrationDF <- data.frame(value = c(as.numeric(valueVec),0,0),
-                              unit = unitVec)
-  # convert the width and height of the image from pixels (0-based so minus 1 from wdth/height) to whatever units the TrackMate file uses
-  calibrationDF[3:4,1] <- (calibrationDF[3:4,1] - 1) * calibrationDF[1,1]
-  # use s for seconds
-  calibrationDF[2,2] <- ifelse(calibrationDF[2,2] == "sec", "s", calibrationDF[2,2])
-  # readout what the units were and warn if spatial units are pixels
-  cat("Units are: ",calibrationDF[1,1],calibrationDF[1,2],"and",calibrationDF[2,1],calibrationDF[2,2],"\n")
-  if(unitVec[1] == "pixel") {
+  filtered_nodes <- xml_find_all(e, ".//TrackID")
+  spot_nodes <- xml_find_all(e, ".//AllSpots//SpotsInFrame//Spot")
+  feature_nodes <- xml_find_all(e, ".//FeatureDeclarations//SpotFeatures//Feature")
+  attrName <- c("name", xml_attr(feature_nodes, "feature"))
+
+  # Units and calibration
+  model_node <- xml_find_first(e, ".//Model")
+  unitVec <- c(xml_attr(model_node, "spatialunits"), xml_attr(model_node, "timeunits"), "widthpixels", "heightpixels", "ntraces", "maxframes")
+  image_node <- xml_find_first(e, ".//ImageData")
+  valueVec <- c(xml_attr(image_node, "pixelwidth"), xml_attr(image_node, "timeinterval"), xml_attr(image_node, "width"), xml_attr(image_node, "height"))
+  calibrationDF <- data.frame(value = as.numeric(c(valueVec, 0, 0)), unit = unitVec)
+  calibrationDF[3:4, 1] <- (calibrationDF[3:4, 1] - 1) * calibrationDF[1, 1]
+  calibrationDF[2, 2] <- ifelse(calibrationDF[2, 2] == "sec", "s", calibrationDF[2, 2])
+  cat("Units are: ", calibrationDF[1, 1], calibrationDF[1, 2], "and", calibrationDF[2, 1], calibrationDF[2, 2], "\n")
+  if (unitVec[1] == "pixel") {
     cat("Spatial units are in pixels - consider transforming to real units\n")
   }
-  # which channel is being tracked?
-  targetVec <-  xpathSApply(e, "//DetectorSettings", xmlGetAttr, "TARGET_CHANNEL")
-  calibrationDF[7,1] <- as.numeric(targetVec)
-  calibrationDF[7,2] <- "channel"
+  detector_node <- xml_find_first(e, ".//DetectorSettings")
+  target_channel <- xml_attr(detector_node, "TARGET_CHANNEL")
+  calibrationDF[7, 1] <- as.numeric(target_channel)
+  calibrationDF[7, 2] <- "channel"
 
-  # if we are doing slim processing, we only get the minimum data required to process the tracks
-  if(slim) {
-    # we only need a subset of possible attributes
-    slimAttr <- c("name", "POSITION_X", "POSITION_Y", "POSITION_Z",
-                  "POSITION_T", "FRAME", "MEAN_INTENSITY",
-                  paste0("MEAN_INTENSITY_CH",targetVec))
+  if (slim) {
+    slimAttr <- c("name", "POSITION_X", "POSITION_Y", "POSITION_Z", "POSITION_T", "FRAME", "MEAN_INTENSITY", paste0("MEAN_INTENSITY_CH", target_channel))
     attrName <- attrName[attrName %in% slimAttr]
   }
 
-  # multicore processing
-  numCores <- parallelly::availableCores()
-
-  if (.Platform[["OS.type"]] == "windows") {
-    ## PSOCK-based parallel processing
-    cl <- parallel::makeCluster(numCores)
-    on.exit(parallel::stopCluster(cl))
-    registerDoParallel(cl = cl)
-  } else {
-    ## Forked parallel processing
-    registerDoParallel(cores = numCores)
-  }
-
-  # perform parallel read
-  # test if we are running on windows
-  if (.Platform[["OS.type"]] == "windows") {
-    cat("Collecting spot data...\n")
-    dtf <- as.data.frame(foreach(i = 1:length(attrName), .packages = c("foreach","XML"), .combine = cbind) %do% {
-      sapply(subdoc, xmlGetAttr, attrName[i])
-    })
-  } else {
-    cat(paste0("Collecting spot data. Using ",numCores," cores\n"))
-    dtf <- as.data.frame(foreach(i = 1:length(attrName), .combine = cbind) %dopar% {
-      sapply(subdoc, xmlGetAttr, attrName[i])
-    })
+  # Extract spot attributes efficiently
+  spot_data <- lapply(attrName, function(attr) xml_attr(spot_nodes, attr))
+  dtf <- as.data.frame(spot_data, stringsAsFactors = FALSE)
+  for (i in 2:length(attrName)) {
+    suppressWarnings(dtf[, i] <- as.numeric(dtf[, i]))
   }
-
-  for (i in 2:length(attrName)){
-    suppressWarnings(dtf[,i] <- as.numeric(as.character(dtf[,i])))
-  }
-  # more R-like headers
   headerNames <- tolower(attrName)
-  # in the case of slim = TRUE, we rename the MEAN_INTENSITY_CHX column to
-  # remove _CH1 or whatever from any header
-  if(slim) {
-    headerNames <- gsub("_ch\\d$","",headerNames,ignore.case = TRUE)
+  if (slim) {
+    headerNames <- gsub("_ch\\d$", "", headerNames, ignore.case = TRUE)
   }
-  # change x y z t
-  headerNames <- gsub("^position\\w","",headerNames,ignore.case = TRUE)
+  headerNames <- gsub("^position\\w", "", headerNames, ignore.case = TRUE)
   names(dtf) <- headerNames
 
   cat("Matching track data...\n")
-
-  # trace is an alternative name for track
-  IDtrace <- data.frame(name = NA, trace = NA, displacement = NA, speed = NA)
-
-  for (i in seq(along = track)){
-    subDoc = xmlDoc(track[[i]])
-    IDvec <- unique(c(unlist(xpathApply(subDoc, "//Edge", xmlGetAttr, "SPOT_SOURCE_ID")),
-                      unlist(xpathApply(subDoc, "//Edge", xmlGetAttr, "SPOT_TARGET_ID"))))
-    traceVec <- rep(sapply(track[i], function(el){xmlGetAttr(el, "TRACK_ID")}), length(IDvec))
-    traceDF <- data.frame(name = paste0("ID",IDvec), trace = traceVec)
-    # now retrieve displacement and speed for target spots in track
-    targetVec <- unlist(xpathApply(subDoc, "//Edge", xmlGetAttr, "SPOT_TARGET_ID"))
-    dispVec <- unlist(xpathApply(subDoc, "//Edge", xmlGetAttr, "DISPLACEMENT"))
-    speedVec <- unlist(xpathApply(subDoc, "//Edge", xmlGetAttr, "SPEED"))
-    # if the user has not added Analyzers, displacement and speed will be NULL
-    if(is.null(dispVec)) {
-      dispVec <- rep(0, length(targetVec))
-    }
-    if(is.null(speedVec)) {
-      speedVec <- rep(0, length(targetVec))
-    }
-    dataDF <- data.frame(name = paste0("ID",targetVec), displacement = as.numeric(dispVec), speed = as.numeric(speedVec))
-    # left join (will give NA for the first spot)
-    allDF <- merge(x = traceDF, y = dataDF, by = "name", all.x = TRUE)
-    allDF$displacement <- ifelse(is.na(allDF$displacement), 0, allDF$displacement)
-    allDF$speed <- ifelse(is.na(allDF$speed), 0, allDF$speed)
-    # grow the dataframe
-    IDtrace <- rbind(IDtrace, allDF)
+  # Preallocate list for track info
+  IDtrace_list <- vector("list", length(track_nodes))
+  for (i in seq_along(track_nodes)) {
+    tr_node <- track_nodes[[i]]
+    edge_nodes <- xml_find_all(tr_node, ".//Edge")
+    source_ids <- xml_attr(edge_nodes, "SPOT_SOURCE_ID")
+    target_ids <- xml_attr(edge_nodes, "SPOT_TARGET_ID")
+    IDvec <- unique(c(source_ids, target_ids))
+    trace_id <- xml_attr(tr_node, "TRACK_ID")
+    traceDF <- data.frame(name = paste0("ID", IDvec), trace = trace_id, stringsAsFactors = FALSE)
+    dispVec <- xml_attr(edge_nodes, "DISPLACEMENT")
+    speedVec <- xml_attr(edge_nodes, "SPEED")
+    # handle missing analyzers
+    if (is.null(dispVec)) dispVec <- rep(0, length(target_ids))
+    if (is.null(speedVec)) speedVec <- rep(0, length(target_ids))
+    dataDF <- data.frame(name = paste0("ID", target_ids), displacement = as.numeric(dispVec), speed = as.numeric(speedVec), stringsAsFactors = FALSE)
+    allDF <- merge(traceDF, dataDF, by = "name", all.x = TRUE)
+    allDF$displacement[is.na(allDF$displacement)] <- 0
+    allDF$speed[is.na(allDF$speed)] <- 0
+    IDtrace_list[[i]] <- allDF
   }
+  IDtrace <- do.call(rbind, IDtrace_list)
 
-  # merge track information with spot data
-  daten <- merge(IDtrace, dtf, by="name")
-  # now we subset for filtered tracks
-  FTvec <- sapply(filtered, xmlGetAttr, "TRACK_ID")
-  daten <- subset(daten, trace %in% FTvec)
-  # sort final data frame
-  daten <- daten[order(daten$trace, daten$t),]
+  daten <- merge(IDtrace, dtf, by = "name")
+  FTvec <- xml_attr(filtered_nodes, "TRACK_ID")
+  daten <- daten[daten$trace %in% FTvec, ]
+  daten <- daten[order(daten$trace, daten$t), ]
 
   cat("Calculating distances...\n")
-
-  # cumulative distance and duration
-  cumdist <- numeric()
+  cumdist <- numeric(nrow(daten))
+  dur <- numeric(nrow(daten))
   cumdist[1] <- 0
-  dur <- numeric()
   dur[1] <- 0
   startdur <- daten$t[1]
-
-  for (i in 2:nrow(daten)){
-    if(daten$trace[i] == daten$trace[i-1]) {
-      cumdist[i] <- cumdist[i-1] + daten$displacement[i]
-    }else{
+  for (i in 2:nrow(daten)) {
+    if (daten$trace[i] == daten$trace[i - 1]) {
+      cumdist[i] <- cumdist[i - 1] + daten$displacement[i]
+    } else {
       cumdist[i] <- 0
       startdur <- daten$t[i]
     }
@@ -164,41 +116,24 @@ readTrackMateXML<- function(XMLpath, slim = FALSE){
   daten$cumulative_distance <- cumdist
   daten$track_duration <- dur
 
-  # users have reported tracks/traces with multiple spots per frame which cannot be processed
-  # possibly caused by track splitting or merging. Produce warning about this.
-  b <- daten %>% group_by(trace,frame) %>%
-    reframe(n = n())
-  # because it is not possible to have n = 0 in this column (as it is reframe) we can do
-  if(sum(b$n) > nrow(b)) {
+  # Check for multiple spots per frame
+  b <- aggregate(name ~ trace + frame, daten, length)
+  if (sum(b$name) > nrow(b)) {
     cat("Warning: Detected multiple spots per frame for one or more tracks.\nTrackMateR will only process single tracks. Subsequent analysis will likely fail!\n")
   }
 
-  # it is possible that xy coords lie outside the image(!)
-  # we can detect xy coords that are less than 0,0 and then use this information to offset *all* coords by this
-  # this is necessary because later code relies on the origin
+  # Offset coordinates if needed
   minx <- min(daten$x)
   miny <- min(daten$y)
-  if(minx < 0) {
-    daten$x <- daten$x - minx
-  }
-  if(miny < 0) {
-    daten$y <- daten$y - miny
-  }
-  # now we need to redefine the size of the "image" because xy coords may lie outside, or they may now lie outside after offsetting
+  if (minx < 0) daten$x <- daten$x - minx
+  if (miny < 0) daten$y <- daten$y - miny
   maxx <- max(daten$x)
   maxy <- max(daten$y)
-  if(maxx > calibrationDF[3,1]) {
-    calibrationDF[3,1] <- ceiling(maxx)
-  }
-  if(maxy > calibrationDF[4,1]) {
-    calibrationDF[4,1] <- ceiling(maxy)
-  }
-  # we need to know how many traces we have and how many frames is the longest one for later
-  calibrationDF[5,1] <- length(unique(daten$trace))
-  calibrationDF[6,1] <- max(daten$track_duration) / calibrationDF[2,1]
-
-  # daten is our dataframe of all data, calibrationDF is the calibration data
-  dfList <- list(daten,calibrationDF)
+  if (maxx > calibrationDF[3, 1]) calibrationDF[3, 1] <- ceiling(maxx)
+  if (maxy > calibrationDF[4, 1]) calibrationDF[4, 1] <- ceiling(maxy)
+  calibrationDF[5, 1] <- length(unique(daten$trace))
+  calibrationDF[6, 1] <- max(daten$track_duration) / calibrationDF[2, 1]
 
+  dfList <- list(daten, calibrationDF)
   return(dfList)
 }

From 934dbca8083bd93ac3164b8f029a28bc4b70fa98 Mon Sep 17 00:00:00 2001
From: Stephen Royle <s.j.royle@warwick.ac.uk>
Date: Tue, 16 Sep 2025 13:54:46 +0100
Subject: [PATCH 2/4] v.0.3.12 overhaul of compareDatasets()

---
 NAMESPACE              |   2 +
 NEWS.md                |   3 +-
 R/TrackMateR-package.R |   2 +
 R/compareDatasets.R    | 186 ++++++++++++++++++++++-------------------
 R/fittingJD.R          |   4 +
 R/makeComparison.R     |   2 +-
 6 files changed, 110 insertions(+), 89 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 3519d8b..39b5776 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -36,6 +36,8 @@ import(patchwork)
 importFrom(ggforce,geom_sina)
 importFrom(graphics,frame)
 importFrom(graphics,hist)
+importFrom(parallel,detectCores)
+importFrom(parallel,mclapply)
 importFrom(reshape2,melt)
 importFrom(stats,aggregate)
 importFrom(stats,approx)
diff --git a/NEWS.md b/NEWS.md
index 8d75350..af3a90f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,7 +3,8 @@
 # TrackMateR 0.3.12
 
 - Moved from {XML} to {xml2} for reading TrackMate XML files. Future-proof but currently slower.
-- Speed improvements in analysis of large datasets
+- Large speed improvement when analysing multiple datasets.
+- Better memory management when analysing large datasets.
 
 # TrackMateR 0.3.11
 
diff --git a/R/TrackMateR-package.R b/R/TrackMateR-package.R
index 177c01c..1d7190b 100644
--- a/R/TrackMateR-package.R
+++ b/R/TrackMateR-package.R
@@ -1,6 +1,8 @@
 #' @importFrom ggforce geom_sina
 #' @importFrom graphics frame
 #' @importFrom graphics hist
+#' @importFrom parallel detectCores
+#' @importFrom parallel mclapply
 #' @importFrom reshape2 melt
 #' @importFrom stats aggregate
 #' @importFrom stats approx
diff --git a/R/compareDatasets.R b/R/compareDatasets.R
index 5a7a205..a4d8527 100644
--- a/R/compareDatasets.R
+++ b/R/compareDatasets.R
@@ -19,7 +19,10 @@
 #' @export
 compareDatasets <- function(...) {
 
-  condition <- value <- dataid <- cumulative_distance <- track_duration <- mean_intensity <- NULL
+  condition <- value <- dataid <- cumulative_distance <- track_duration <- mean_intensity <- calibrationDF <- timeRes <-NULL
+  megamsd <- megaalpha <- megadee <- megatd <- megaspeed <- megafd <- NULL
+  units_vec <- NULL
+  summary_params <- NULL
 
   if(!dir.exists("Data")) {
     # there is no cross-platform way to safely choose directory
@@ -68,74 +71,67 @@ compareDatasets <- function(...) {
       calibrate <- TRUE
     }
     cat(paste0("\n","Processing ",condFolderName,"\n"))
-    bigtm <- bigmsd <- bigalpha <- bigdee <- bigjd <- bigtd <- bigfd <- data.frame()
-    for(j in 1:length(allTrackMateFiles)) {
+    # Use lists for accumulation
+    bigtm_list <- list()
+    bigcalibration_list <- list()
+    bigmsd_list <- list()
+    bigalpha_list <- list()
+    bigdee_list <- list()
+    bigjd_list <- list()
+    bigjdParams_list <- list()
+    bigtd_list <- list()
+    bigfd_list <- list()
+    bigreport_list <- list()
+
+    # Parallelize inner loop over XML files
+    process_file <- function(j) {
       fileName <- allTrackMateFiles[j]
       thisFilePath <- paste0(condFolderPath, "/", fileName)
-      # read dataset
       tmObj <- readTrackMateXML(XMLpath = thisFilePath, slim = TRUE)
       if(is.null(tmObj)) {
         cat(paste0("Skipping ",fileName," - no data found!\n"))
-        next
+        return(NULL)
       }
-      # scale dataset if required
       if(calibrate) {
         calibrationDF <- tmObj[[2]]
-        # scalar for conversion is new / old (units not relevant)
         calibrationXY <- calibDF[1,1] / calibrationDF[1,1]
         calibrationT <- calibDF[2,1] / calibrationDF[2,1]
-        # 0 in calibDF indicates no scaling is to be done
         calibrationXY <- ifelse(calibrationXY == 0, 1, calibrationXY)
         calibrationT <- ifelse(calibrationT == 0, 1, calibrationT)
-        # ignore an error of 2.5%
         calibrationXY <- ifelse(calibrationXY < 1.025 & calibrationXY > 0.975, 1, calibrationXY)
         calibrationT <- ifelse(calibrationT < 1.025 & calibrationT > 0.975, 1, calibrationT)
-        if(calibrationXY != 1 & calibrationT != 1) {
+  if(calibrationXY != 1 && calibrationT != 1) {
           tmObj <- correctTrackMateData(dataList = tmObj, xyscalar = calibrationXY, tscalar = calibrationT, xyunit = calibDF[1,2], tunit = calibDF[2,2])
-        } else if(calibrationXY != 1 & calibrationT == 1) {
+  } else if(calibrationXY != 1 && calibrationT == 1) {
           tmObj <- correctTrackMateData(dataList = tmObj, xyscalar = calibrationXY, xyunit = calibDF[1,2])
-        } else if(calibrationXY == 1 & calibrationT != 1) {
+  } else if(calibrationXY == 1 && calibrationT != 1) {
           tmObj <- correctTrackMateData(dataList = tmObj, tscalar = calibrationT, tunit = calibDF[2,2])
         } else {
-          # the final possibility is nothing needs scaling but units need to change.
-          # do not test if they are the same just flush the units with the values in the csv file
           tmObj <- correctTrackMateData(dataList = tmObj, xyunit = calibDF[1,2], tunit = calibDF[2,2])
         }
       }
-      # we can filter here if required - for example only analyse tracks of certain length
       tmDF <- tmObj[[1]]
       calibrationDF <- tmObj[[2]]
-      # sanity check - probably not needed
       if(is.null(tmDF)) {
         cat(paste0("Skipping ",fileName," - no data found!\n"))
-        next
+        return(NULL)
       }
-      # if the data is not rich enough for a summary we will skip it
       if(calibrationDF[5,1] < 3 & calibrationDF[6,1] < 10) {
         cat(paste0("Skipping ",fileName," as it has less than 3 tracks and longest track has less than 10 frames\n"))
-        next
+        return(NULL)
       }
-      # take the units
-      units <- calibrationDF$unit[1:2]
-
-      ## we need to combine data frames
-      # first add a column to id the data
+  units_vec <- calibrationDF$unit[1:2]
       thisdataid <- paste0(condFolderName,"_",as.character(j))
       tmDF$dataid <- thisdataid
-
-      # calculate MSD
       msdObj <- calculateMSD(tmDF, N = 3, short = 8)
       msdDF <- msdObj[[1]]
       alphaDF <- msdObj[[2]]
       deeDF <- msdObj[[3]]
-      if(!is.null(msdDF) | !is.null(alphaDF) | !is.null(deeDF)) {
-        # we need to add the dataid to the summary
+  if(!is.null(msdDF) || !is.null(alphaDF) || !is.null(deeDF)) {
         msdDF$dataid <- thisdataid
         alphaDF$dataid <- thisdataid
         deeDF$dataid <- thisdataid
       }
-
-      # jump distance calc with deltaT of 1
       deltaT <- 1
       jdObj <- calculateJD(dataList = tmObj, deltaT = l$deltaT, nPop = l$nPop, mode = l$mode, init = l$init, timeRes = l$timeRes, breaks = l$breaks)
       jdDF <- jdObj[[1]]
@@ -146,74 +142,88 @@ compareDatasets <- function(...) {
         timeRes <- jdObj[[2]]
         jdObj <- list(jdDF,timeRes)
       }
-      # track density with a radius of 1.5 units
       tdDF <- calculateTrackDensity(dataList = tmObj, radius = l$radius)
       if(!is.null(tdDF)) {
         tdDF$dataid <- thisdataid
       }
-      # fractal dimension
       fdDF <- calculateFD(dataList = tmObj)
       if(!is.null(fdDF)) {
         fdDF$dataid <- thisdataid
       }
-
-      # add to the big dataframes
-      if(!is.null(tmDF)) {
-        bigtm <- rbind(bigtm,tmDF)
-      }
-      if(!is.null(msdDF)) {
-        bigmsd <- rbind(bigmsd,msdDF)
-      }
-      if(!is.null(alphaDF)) {
-        bigalpha <- rbind(bigalpha,alphaDF)
-      }
-      if(!is.null(deeDF)) {
-        bigdee <- rbind(bigdee,deeDF)
-      }
-      if(!is.null(jdDF)) {
-        bigjd <- rbind(bigjd,jdDF)
-      }
-      if(!is.null(tdDF)) {
-        bigtd <- rbind(bigtd,tdDF)
-      }
-      if(!is.null(fdDF)) {
-        bigfd <- rbind(bigfd,fdDF)
-      }
-
-      # create the report for this dataset
-      fileName <- tools::file_path_sans_ext(basename(thisFilePath))
+      fileNameBase <- tools::file_path_sans_ext(basename(thisFilePath))
       both <- makeSummaryReport(tmList = tmObj, msdList = msdObj, jumpList = jdObj, tddf = tdDF, fddf = fdDF,
-                                titleStr = condFolderName, subStr = fileName, auto = TRUE, summary = FALSE,
+                                titleStr = condFolderName, subStr = fileNameBase, auto = TRUE, summary = FALSE,
                                 msdplot = l$msdplot)
       p <- both[[1]]
       destinationDir <- paste0("Output/Plots/", condFolderName)
       setupOutputPath(destinationDir)
       filePath <- paste0(destinationDir, "/report_",as.character(j),".pdf")
       ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
-      # retrieve other data
       df_report <- both[[2]]
       df_report$condition <- condFolderName
       df_report$dataid <- thisdataid
-      if(i == 1 & j == 1) {
-        megareport <- df_report
-      } else if(!exists("megareport")) {
-        megareport <- df_report
-      } else {
-        megareport <- rbind(megareport,df_report)
-      }
+      return(list(tmDF = tmDF, calibrationDF = calibrationDF, msdDF = msdDF, alphaDF = alphaDF, deeDF = deeDF, jdDF = jdDF, jdParams = timeRes, tdDF = tdDF, fdDF = fdDF, df_report = df_report))
+    }
+
+    # Use mclapply for parallel processing (on non-Windows)
+    results <- if (.Platform$OS.type == "windows") {
+      lapply(seq_along(allTrackMateFiles), process_file)
+    } else {
+      mclapply(seq_along(allTrackMateFiles), process_file, mc.cores = detectCores())
+    }
+
+    for (res in results) {
+      if (is.null(res)) next
+      if (!is.null(res$tmDF)) bigtm_list[[length(bigtm_list)+1]] <- res$tmDF
+      if (!is.null(res$calibrationDF)) bigcalibration_list[[length(bigcalibration_list)+1]] <- res$calibrationDF
+      if (!is.null(res$msdDF)) bigmsd_list[[length(bigmsd_list)+1]] <- res$msdDF
+      if (!is.null(res$alphaDF)) bigalpha_list[[length(bigalpha_list)+1]] <- res$alphaDF
+      if (!is.null(res$deeDF)) bigdee_list[[length(bigdee_list)+1]] <- res$deeDF
+      if (!is.null(res$jdDF)) bigjd_list[[length(bigjd_list)+1]] <- res$jdDF
+      if (!is.null(res$jdParams)) bigjdParams_list[[length(bigjdParams_list)+1]] <- res$jdParams
+      if (!is.null(res$tdDF)) bigtd_list[[length(bigtd_list)+1]] <- res$tdDF
+      if (!is.null(res$fdDF)) bigfd_list[[length(bigfd_list)+1]] <- res$fdDF
+      if (!is.null(res$df_report)) bigreport_list[[length(bigreport_list)+1]] <- res$df_report
     }
-    bigtmObj <- list(bigtm,calibrationDF)
-    bigmsdObj <- list(bigmsd,bigalpha,bigdee)
-    bigjdObj <- list(bigjd,timeRes)
-    # now we have our combined dataset we can make a summary
-    # note we use the timeRes of the final dataset; so it is suitable for only when all files have the same calibration
-    summaryObj <- makeSummaryReport(tmList = bigtmObj, msdList = bigmsdObj, jumpList = bigjdObj, tddf = bigtd, fddf = bigfd,
-                                    titleStr = condFolderName, subStr = "Summary", auto = TRUE, summary = TRUE,
-                                    msdplot = l$msdplot)
-    p <- summaryObj[[1]]
-    destinationDir <- paste0("Output/Plots/", condFolderName)
-    filePath <- paste0(destinationDir, "/combined.pdf")
-    ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
+
+  bigtm <- if (length(bigtm_list) > 0) do.call(rbind, bigtm_list) else data.frame()
+  bigmsd <- if (length(bigmsd_list) > 0) do.call(rbind, bigmsd_list) else data.frame()
+  bigalpha <- if (length(bigalpha_list) > 0) do.call(rbind, bigalpha_list) else data.frame()
+  bigdee <- if (length(bigdee_list) > 0) do.call(rbind, bigdee_list) else data.frame()
+  bigjd <- if (length(bigjd_list) > 0) do.call(rbind, bigjd_list) else data.frame()
+  bigtd <- if (length(bigtd_list) > 0) do.call(rbind, bigtd_list) else data.frame()
+  bigfd <- if (length(bigfd_list) > 0) do.call(rbind, bigfd_list) else data.frame()
+  bigreport <- if (length(bigreport_list) > 0) do.call(rbind, bigreport_list) else data.frame()
+
+  # Fix: collect all valid parameter lists from results and use the first one for summary
+    # Extract units_vec from first valid result
+if (is.null(units_vec)) {
+  for (res in results) {
+    if (!is.null(res) && !is.null(res$calibrationDF)) {
+      units_vec <- res$calibrationDF$unit[1:2]
+      break
+    }
+  }
+}
+
+# Extract summary_params from first valid result
+all_param_lists <- lapply(results, function(res) {
+  if (!is.null(res) && !is.null(res$jdParams)) res$jdParams else NULL
+})
+valid_param_lists <- Filter(Negate(is.null), all_param_lists)
+if (is.null(summary_params) || length(summary_params) == 0) {
+  summary_params <- if (length(valid_param_lists) > 0) valid_param_lists[[1]] else list()
+}
+  bigtmObj <- list(bigtm,calibrationDF)
+  bigmsdObj <- list(bigmsd,bigalpha,bigdee)
+  bigjdObj <- list(bigjd, summary_params)
+  summaryObj <- makeSummaryReport(tmList = bigtmObj, msdList = bigmsdObj, jumpList = bigjdObj, tddf = bigtd, fddf = bigfd,
+                  titleStr = condFolderName, subStr = "Summary", auto = TRUE, summary = TRUE,
+                  msdplot = l$msdplot)
+  p <- summaryObj[[1]]
+  destinationDir <- paste0("Output/Plots/", condFolderName)
+  filePath <- paste0(destinationDir, "/combined.pdf")
+  ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
     # save data as csv
     destinationDir <- paste0("Output/Data/", condFolderName)
     setupOutputPath(destinationDir)
@@ -230,20 +240,22 @@ compareDatasets <- function(...) {
       summarise(cumdist = max(cumulative_distance), cumtime = max(track_duration), intensity = max(mean_intensity))
     bigspeed$speed <- bigspeed$cumdist / bigspeed$cumtime
     bigspeed$condition <- condFolderName
-    if(i == 1 | !exists("megamsd")) {
+    if (is.null(megamsd)) {
       megamsd <- msdSummary
       megaalpha <- bigalpha
       megadee <- bigdee
       megatd <- bigtd
       megaspeed <- bigspeed
       megafd <- bigfd
+      megareport <- bigreport
     } else {
-      megamsd <- rbind(megamsd,msdSummary)
-      megaalpha <- rbind(megaalpha,bigalpha)
-      megadee <- rbind(megadee,bigdee)
-      megatd <- rbind(megatd,bigtd)
-      megaspeed <- rbind(megaspeed,bigspeed)
-      megafd <- rbind(megafd,bigfd)
+      megamsd <- rbind(megamsd, msdSummary)
+      megaalpha <- rbind(megaalpha, bigalpha)
+      megadee <- rbind(megadee, bigdee)
+      megatd <- rbind(megatd, bigtd)
+      megaspeed <- rbind(megaspeed, bigspeed)
+      megafd <- rbind(megafd, bigfd)
+      megareport <- rbind(megareport, bigreport)
     }
   }
 
@@ -259,7 +271,7 @@ compareDatasets <- function(...) {
   write.csv(megatrace, paste0(destinationDir, "/allTraceData.csv"), row.names = FALSE)
 
   # generate the comparison plots and save
-  p <- makeComparison(df = megareport, msddf = megamsd, units = units, msdplot = l$msdplot)
+  p <- makeComparison(df = megareport, msddf = megamsd, units = units_vec, msdplot = l$msdplot)
   destinationDir <- "Output/Plots/"
   filePath <- paste0(destinationDir, "/comparison.pdf")
   ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
diff --git a/R/fittingJD.R b/R/fittingJD.R
index c2ed8bf..012e39c 100644
--- a/R/fittingJD.R
+++ b/R/fittingJD.R
@@ -41,6 +41,10 @@ fittingJD <- function(jumpList) {
   timeRes <- params$timeRes
   breaks <- params$breaks
 
+  if(is.null(nPop)) {
+    return(NULL)
+  }
+
   if(nPop < 1 | nPop > 3) {
     return(NULL)
   }
diff --git a/R/makeComparison.R b/R/makeComparison.R
index e642deb..85c1050 100644
--- a/R/makeComparison.R
+++ b/R/makeComparison.R
@@ -66,7 +66,7 @@ makeComparison <- function (df, msddf, units = c("um","s"), msdplot = "linlin",
     geom_sina(alpha = 0.5, stroke = 0) +
     ylim(c(0,NA)) +
     guides(x =  guide_axis(angle = 90)) +
-    labs(x = "", y = substitute(paste("Diffusion coefficient (",mm^2,"/",nn,")"), list(mm = units[1], nn = units [2]))) +
+    labs(x = "", y = substitute(paste("Diffusion coefficient (",mm^2,"/",nn,")"), list(mm = units[1], nn = units[2]))) +
     theme_classic() +
     theme(legend.position = "none")
 

From 662355bf0a016e4a7644a0f14251c4da0375884f Mon Sep 17 00:00:00 2001
From: Stephen Royle <s.j.royle@warwick.ac.uk>
Date: Tue, 16 Sep 2025 16:14:39 +0100
Subject: [PATCH 3/4] v.0.3.12 cleanup if many/longname datasets are used

---
 R/compareDatasets.R   |  9 ++++++++-
 R/makeComparison.R    | 20 +++++++++++++++++++-
 man/makeComparison.Rd |  7 ++++++-
 3 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/R/compareDatasets.R b/R/compareDatasets.R
index a4d8527..61d52ff 100644
--- a/R/compareDatasets.R
+++ b/R/compareDatasets.R
@@ -274,5 +274,12 @@ if (is.null(summary_params) || length(summary_params) == 0) {
   p <- makeComparison(df = megareport, msddf = megamsd, units = units_vec, msdplot = l$msdplot)
   destinationDir <- "Output/Plots/"
   filePath <- paste0(destinationDir, "/comparison.pdf")
-  ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
+  # if there are many conditions, increase width of plot
+  if(length(condFolderNames) < 3) {
+    ggsave(filePath, plot = p, width = 19, height = 19, units = "cm")
+  } else if (length(condFolderNames) > 6) {
+    ggsave(filePath, plot = p, width = 35, height = 19, units = "cm")
+  } else{
+    ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
+  }
 }
diff --git a/R/makeComparison.R b/R/makeComparison.R
index 85c1050..b06d2a1 100644
--- a/R/makeComparison.R
+++ b/R/makeComparison.R
@@ -1,7 +1,11 @@
 #' Make Comparison Plots
 #'
 #' A series of ggplots to compare between conditions.
-#' Called from `compareDatasets()` this function generates plots using summary data of datasets, per condition.
+#' Called from `compareDatasets()`, this function generates plots using summary
+#' data of datasets, per condition.
+#'
+#' Note that if one dataset has a long name the names will be wrapped. Wrapping
+#' not graceful. If you have very long names consider renaming them before running.
 #'
 #' @param df data frame called megareport
 #' @param msddf data frame of msd averages per dataset
@@ -18,6 +22,20 @@ makeComparison <- function (df, msddf, units = c("um","s"), msdplot = "linlin",
   options(warn = -1)
 
   symlim <- findLog2YAxisLimits(df$alpha)
+  # find longest condition name
+  maxchar <- max(nchar(as.character(df$condition)))
+  lablength <- 15
+  if(maxchar > lablength) {
+    wrapit <- TRUE
+  } else {
+    wrapit <- FALSE
+  }
+  if(wrapit) {
+    # this is a hack but it deals with labels with no spaces (which is what we want)
+    # insert a newline every 20 characters
+    expr <- paste0("(.{",lablength,"})")
+    df$condition <- gsub(expr, "\\1\n", df$condition)
+  }
 
   # plot alpha comparison
   p_alpha <- ggplot(data = df, aes(x = condition, y = alpha, colour = condition)) +
diff --git a/man/makeComparison.Rd b/man/makeComparison.Rd
index 2fb31c1..f29e45c 100644
--- a/man/makeComparison.Rd
+++ b/man/makeComparison.Rd
@@ -31,5 +31,10 @@ patchwork ggplot
 }
 \description{
 A series of ggplots to compare between conditions.
-Called from `compareDatasets()` this function generates plots using summary data of datasets, per condition.
+Called from `compareDatasets()`, this function generates plots using summary
+data of datasets, per condition.
+}
+\details{
+Note that if one dataset has a long name the names will be wrapped. Wrapping
+not graceful. If you have very long names consider renaming them before running.
 }

From c8d5e7d1034e5e0e9f9d718d9457236bc91f3f30 Mon Sep 17 00:00:00 2001
From: Stephen Royle <s.j.royle@warwick.ac.uk>
Date: Wed, 17 Sep 2025 16:53:10 +0100
Subject: [PATCH 4/4] Fix: crash on large files, limit cores to n - 2, fix
 units on condition-level summary

---
 R/compareDatasets.R | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/R/compareDatasets.R b/R/compareDatasets.R
index 61d52ff..9ec4893 100644
--- a/R/compareDatasets.R
+++ b/R/compareDatasets.R
@@ -169,7 +169,7 @@ compareDatasets <- function(...) {
     results <- if (.Platform$OS.type == "windows") {
       lapply(seq_along(allTrackMateFiles), process_file)
     } else {
-      mclapply(seq_along(allTrackMateFiles), process_file, mc.cores = detectCores())
+      mclapply(seq_along(allTrackMateFiles), process_file, mc.cores = detectCores() - 2)
     }
 
     for (res in results) {
@@ -205,6 +205,9 @@ if (is.null(units_vec)) {
     }
   }
 }
+# package into a data frame that can be used by makeSummaryReport() to get the units
+dummyCalibrationDF <- data.frame(value = c(1, 1),
+                             unit = units_vec)
 
 # Extract summary_params from first valid result
 all_param_lists <- lapply(results, function(res) {
@@ -214,8 +217,8 @@ valid_param_lists <- Filter(Negate(is.null), all_param_lists)
 if (is.null(summary_params) || length(summary_params) == 0) {
   summary_params <- if (length(valid_param_lists) > 0) valid_param_lists[[1]] else list()
 }
-  bigtmObj <- list(bigtm,calibrationDF)
-  bigmsdObj <- list(bigmsd,bigalpha,bigdee)
+  bigtmObj <- list(bigtm, dummyCalibrationDF)
+  bigmsdObj <- list(bigmsd, bigalpha, bigdee)
   bigjdObj <- list(bigjd, summary_params)
   summaryObj <- makeSummaryReport(tmList = bigtmObj, msdList = bigmsdObj, jumpList = bigjdObj, tddf = bigtd, fddf = bigfd,
                   titleStr = condFolderName, subStr = "Summary", auto = TRUE, summary = TRUE,
@@ -279,7 +282,7 @@ if (is.null(summary_params) || length(summary_params) == 0) {
     ggsave(filePath, plot = p, width = 19, height = 19, units = "cm")
   } else if (length(condFolderNames) > 6) {
     ggsave(filePath, plot = p, width = 35, height = 19, units = "cm")
-  } else{
+  } else {
     ggsave(filePath, plot = p, width = 25, height = 19, units = "cm")
   }
 }