From ba4c81166b3258a5f581c4e6b2a82ac5d0be58fc Mon Sep 17 00:00:00 2001 From: Gene233 Date: Mon, 20 Apr 2026 18:42:18 +1000 Subject: [PATCH 1/6] Phase A: cal_score entry-point memory cleanup * Drop forced as.matrix(data) in cal_score() matrix method so that dgCMatrix inputs now flow through to the existing sparse-aware row/col aggregators instead of being densified on entry. * Add return.intermediate = FALSE to cal_score() (both methods) and cal_score_init(). Default path no longer stores the full-size tf/idf/iae matrices in SummarizedExperiment metadata, which on a 20k x 100k input saved roughly 3x the dense-matrix footprint per call. Explicitly clear stale metadata when the flag is FALSE to avoid leakage from prior runs. * Remove the df_n_inv <- !df_n dead code in idf_prob(); it was never used downstream and carried a logical matrix the same size as expr. * Add tests/testthat/test-numerical-equivalence.R and the frozen legacy_scores.rds snapshot captured from the pre-refactor HEAD. These pin cal_score() and top_markers() outputs (GLM, GLM+batch, abs-mean, abs-median) at 1e-10 tolerance for all downstream phases. Breaking change: return.intermediate defaults to FALSE. Callers that relied on metadata(se)$tf / $idf / $iae must pass return.intermediate = TRUE. Documented in the roxygen @param block; NEWS entry deferred to Phase D so the user-visible announcement lands alongside the full refactor. --- R/AllGenerics.R | 9 +- R/score-methods.R | 30 ++- R/score.R | 19 +- R/tf_idf_iae_wrappers.R | 1 - tests/testthat/test-numerical-equivalence.R | 242 ++++++++++++++++++ .../testdata/capture_legacy_snapshots.R | 139 ++++++++++ tests/testthat/testdata/legacy_scores.rds | Bin 0 -> 44717 bytes 7 files changed, 426 insertions(+), 14 deletions(-) create mode 100644 tests/testthat/test-numerical-equivalence.R create mode 100644 tests/testthat/testdata/capture_legacy_snapshots.R create mode 100644 tests/testthat/testdata/legacy_scores.rds diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 97ee8c1..ebe70b4 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -20,6 +20,12 @@ #' optional, default 'counts' #' @param new.slot a character, specify the name of slot to save score in se object, #' optional, default 'score' +#' @param return.intermediate logical, if TRUE also return or store the +#' intermediate `tf`, `idf` and `iae` matrices. Defaults to FALSE since +#' these objects have the same dimension as the input expression matrix +#' and can dominate memory usage on large datasets. Set to TRUE to +#' restore the pre-1.7.3 behavior where intermediates were kept in +#' `metadata()` of the SummarizedExperiment output. #' #' @return A list of matrices or se object containing combined score #' @@ -42,7 +48,8 @@ setGeneric( slot = "counts", new.slot = "score", par.idf = NULL, - par.iae = NULL) { + par.iae = NULL, + return.intermediate = FALSE) { standardGeneric("cal_score") } ) diff --git a/R/score-methods.R b/R/score-methods.R index 122a278..220797e 100644 --- a/R/score-methods.R +++ b/R/score-methods.R @@ -11,14 +11,16 @@ setMethod( idf = "prob", iae = "prob", par.idf = NULL, - par.iae = NULL) { + par.iae = NULL, + return.intermediate = FALSE) { score <- cal_score_init( - expr = as.matrix(data), + expr = data, tf = tf, idf = idf, iae = iae, par.idf = par.idf, - par.iae = par.iae + par.iae = par.iae, + return.intermediate = return.intermediate ) return(score) @@ -37,7 +39,8 @@ setMethod( slot = "counts", new.slot = "score", par.idf = NULL, - par.iae = NULL) { + par.iae = NULL, + return.intermediate = FALSE) { ## get expr expr <- SummarizedExperiment::assay(data, i = slot) ## get label @@ -54,13 +57,24 @@ setMethod( idf = idf, iae = iae, par.idf = par.idf, - par.iae = par.iae + par.iae = par.iae, + return.intermediate = return.intermediate ) SummarizedExperiment::assay(data, i = new.slot) <- res$score - slot(data, "metadata")$tf <- res$tf - slot(data, "metadata")$idf <- res$idf - slot(data, "metadata")$iae <- res$iae + + ## Store or clear intermediate matrices depending on caller request. + ## Explicitly clearing stale values prevents old metadata from leaking + ## when `cal_score()` is re-run on an object that previously held them. + if (isTRUE(return.intermediate)) { + slot(data, "metadata")$tf <- res$tf + slot(data, "metadata")$idf <- res$idf + slot(data, "metadata")$iae <- res$iae + } else { + slot(data, "metadata")$tf <- NULL + slot(data, "metadata")$idf <- NULL + slot(data, "metadata")$iae <- NULL + } return(data) } diff --git a/R/score.R b/R/score.R index 786adfa..d0df49b 100644 --- a/R/score.R +++ b/R/score.R @@ -57,8 +57,13 @@ idf_iae_methods <- function() { #' be accessed using [idf_iae_methods()] #' @param par.idf other parameters for specified IDF methods #' @param par.iae other parameters for specified IAE methods +#' @param return.intermediate logical, if TRUE the returned list also contains +#' the intermediate `tf`, `idf` and `iae` objects. Default `FALSE` keeps +#' only the combined `score` to avoid the memory overhead of three extra +#' feature-by-cell matrices on large inputs. #' -#' @return a list of combined score, tf, idf and iae +#' @return a list always containing `score`; when `return.intermediate = TRUE` +#' the list additionally contains `tf`, `idf` and `iae`. #' #' @examples #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10)) @@ -69,14 +74,17 @@ idf_iae_methods <- function() { #' ) cal_score_init <- function(expr, tf = c("logtf", "tf"), idf = "prob", iae = "prob", - par.idf = NULL, par.iae = NULL) { + par.idf = NULL, par.iae = NULL, + return.intermediate = FALSE) { ## check tf <- match.arg(tf) idf <- match.arg(idf, choices = idf_iae_methods()) iae <- match.arg(iae, choices = idf_iae_methods()) stopifnot( "par.idf must be a named list or NULL" = is.null(par.idf) | is.list(par.idf), - "par.iae must be a named list or NULL" = is.null(par.iae) | is.list(par.iae) + "par.iae must be a named list or NULL" = is.null(par.iae) | is.list(par.iae), + "return.intermediate must be a single logical" = + is.logical(return.intermediate) && length(return.intermediate) == 1L ) ## compute tf @@ -101,5 +109,8 @@ cal_score_init <- function(expr, tf = c("logtf", "tf"), ## combined score score <- tf * idf * iae - return(list(score = score, tf = tf, idf = idf, iae = iae)) + if (isTRUE(return.intermediate)) { + return(list(score = score, tf = tf, idf = idf, iae = iae)) + } + list(score = score) } diff --git a/R/tf_idf_iae_wrappers.R b/R/tf_idf_iae_wrappers.R index 136e1c9..7cea362 100644 --- a/R/tf_idf_iae_wrappers.R +++ b/R/tf_idf_iae_wrappers.R @@ -257,7 +257,6 @@ idf_prob <- function(expr, features = NULL, label, # thres <- 0 # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) df_n <- expr[features, , drop = FALSE] > thres ## if contain feature i > thres - df_n_inv <- !df_n ## if not contain feature i > thres ## convert label into character in case problem for factor label <- as.character(label) diff --git a/tests/testthat/test-numerical-equivalence.R b/tests/testthat/test-numerical-equivalence.R new file mode 100644 index 0000000..33c9d13 --- /dev/null +++ b/tests/testthat/test-numerical-equivalence.R @@ -0,0 +1,242 @@ +## Numerical regression tests for the memory-optimization refactor. +## +## These tests compare current output against pre-refactor snapshots +## produced by `tests/testthat/testdata/capture_legacy_snapshots.R` on the +## same HEAD that first introduced this file. Any refactor that claims to +## preserve scoring behaviour must keep all `all.equal()` calls green. +## +## Tolerance is deliberately tight (1e-10) because the refactor only +## rearranges algebra and subsetting order; no stochastic steps are +## involved upstream of the matrices compared here. + +snap_path <- test_path("testdata", "legacy_scores.rds") + +skip_if_no_snapshot <- function() { + if (!file.exists(snap_path)) { + skip(paste0("Legacy snapshot missing: ", snap_path)) + } +} + +test_that("cal_score on dense matrix reproduces legacy score", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + res <- cal_score( + data = snap$inputs$counts_dense, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = as.character(snap$inputs$label)), + par.iae = list(label = as.character(snap$inputs$label)), + return.intermediate = TRUE + ) + expect_equal( + unname(as.matrix(res$score)), + unname(as.matrix(snap$cal_score$matrix_out$score)), + tolerance = 1e-10 + ) +}) + +test_that("cal_score default no longer stores intermediates in metadata", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = S4Vectors::DataFrame( + Group = snap$inputs$label, + Batch = snap$inputs$batch + ) + ) + out <- cal_score( + data = se, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group") + ) + expect_null(S4Vectors::metadata(out)$tf) + expect_null(S4Vectors::metadata(out)$idf) + expect_null(S4Vectors::metadata(out)$iae) + expect_equal( + unname(as.matrix(SummarizedExperiment::assay(out, "score"))), + unname(snap$cal_score$se_assay), + tolerance = 1e-10 + ) +}) + +test_that("cal_score return.intermediate=TRUE restores legacy metadata", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = S4Vectors::DataFrame( + Group = snap$inputs$label, + Batch = snap$inputs$batch + ) + ) + out <- cal_score( + data = se, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group"), + return.intermediate = TRUE + ) + md <- S4Vectors::metadata(out) + expect_false(is.null(md$tf)) + expect_false(is.null(md$idf)) + expect_false(is.null(md$iae)) + expect_equal( + unname(as.matrix(md$tf)), + unname(snap$cal_score$se_tf), + tolerance = 1e-10 + ) + expect_equal( + unname(as.matrix(md$idf)), + unname(snap$cal_score$se_idf), + tolerance = 1e-10 + ) + expect_equal( + unname(as.matrix(md$iae)), + unname(snap$cal_score$se_iae), + tolerance = 1e-10 + ) +}) + +test_that("top_markers GLM path reproduces legacy scores", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = S4Vectors::DataFrame( + Group = snap$inputs$label, + Batch = snap$inputs$batch + ) + ) + se <- cal_score( + data = se, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group") + ) + tm <- top_markers( + data = se, + label = "Group", + n = 10L, + use.glm = TRUE, + slot = "score" + ) + tm_ref <- snap$top_markers$glm + ## order-sensitive comparison on the same keys + key_new <- paste(tm$.dot, tm$Genes, sep = "|") + key_ref <- paste(tm_ref$.dot, tm_ref$Genes, sep = "|") + expect_setequal(key_new, key_ref) + ## align and compare scores + idx <- match(key_ref, key_new) + expect_equal(tm$Scores[idx], tm_ref$Scores, tolerance = 1e-8) +}) + +test_that("top_markers GLM with batch reproduces legacy scores", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = S4Vectors::DataFrame( + Group = snap$inputs$label, + Batch = snap$inputs$batch + ) + ) + se <- cal_score( + data = se, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group") + ) + tm <- top_markers( + data = se, + label = "Group", + n = 10L, + use.glm = TRUE, + batch = "Batch", + slot = "score" + ) + tm_ref <- snap$top_markers$glm_batch + key_new <- paste(tm$.dot, tm$Genes, sep = "|") + key_ref <- paste(tm_ref$.dot, tm_ref$Genes, sep = "|") + expect_setequal(key_new, key_ref) + idx <- match(key_ref, key_new) + expect_equal(tm$Scores[idx], tm_ref$Scores, tolerance = 1e-8) +}) + +test_that("top_markers abs mean path reproduces legacy scores", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = S4Vectors::DataFrame( + Group = snap$inputs$label, + Batch = snap$inputs$batch + ) + ) + se <- cal_score( + data = se, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group") + ) + tm <- top_markers( + data = se, + label = "Group", + n = 10L, + use.glm = FALSE, + method = "mean", + slot = "score" + ) + tm_ref <- snap$top_markers$abs_mean + key_new <- paste(tm$.dot, tm$Genes, sep = "|") + key_ref <- paste(tm_ref$.dot, tm_ref$Genes, sep = "|") + expect_setequal(key_new, key_ref) + idx <- match(key_ref, key_new) + expect_equal(tm$Scores[idx], tm_ref$Scores, tolerance = 1e-8) +}) + +test_that("top_markers abs median path reproduces legacy scores", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = S4Vectors::DataFrame( + Group = snap$inputs$label, + Batch = snap$inputs$batch + ) + ) + se <- cal_score( + data = se, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group") + ) + tm <- top_markers( + data = se, + label = "Group", + n = 10L, + use.glm = FALSE, + method = "median", + slot = "score" + ) + tm_ref <- snap$top_markers$abs_median + key_new <- paste(tm$.dot, tm$Genes, sep = "|") + key_ref <- paste(tm_ref$.dot, tm_ref$Genes, sep = "|") + expect_setequal(key_new, key_ref) + idx <- match(key_ref, key_new) + expect_equal(tm$Scores[idx], tm_ref$Scores, tolerance = 1e-8) +}) diff --git a/tests/testthat/testdata/capture_legacy_snapshots.R b/tests/testthat/testdata/capture_legacy_snapshots.R new file mode 100644 index 0000000..4809fc1 --- /dev/null +++ b/tests/testthat/testdata/capture_legacy_snapshots.R @@ -0,0 +1,139 @@ +# Capture legacy outputs of cal_score() and top_markers() from the current +# pre-refactor HEAD. Run ONCE from the package root with the source branch +# loaded via devtools::load_all(). The resulting .rds file is committed and +# consumed by tests/testthat/test-numerical-equivalence.R to verify the +# memory-optimization refactor preserves numerical equivalence. +# +# Usage (from package root): +# Rscript tests/testthat/testdata/capture_legacy_snapshots.R +# +# Prerequisites: devtools, Matrix, SummarizedExperiment, S4Vectors installed. + +suppressPackageStartupMessages({ + library(devtools) + library(Matrix) + library(SummarizedExperiment) + library(S4Vectors) +}) + +## Load source package at current HEAD (pre-refactor) +devtools::load_all(".", quiet = TRUE) + +## Deterministic test inputs ---------------------------------------------- +set.seed(2026) + +G <- 40L # features +N <- 60L # cells +K <- 3L # groups + +counts_dense <- matrix(rpois(G * N, lambda = 1.5), nrow = G) +rownames(counts_dense) <- paste0("gene", seq_len(G)) +colnames(counts_dense) <- paste0("cell", seq_len(N)) + +counts_sparse <- as(counts_dense, "CsparseMatrix") # dgCMatrix + +label <- factor(rep(LETTERS[seq_len(K)], length.out = N)) +batch <- factor(rep(c("b1", "b2"), length.out = N)) + +se_dense <- SummarizedExperiment( + assays = list(counts = counts_dense), + colData = DataFrame(Group = label, Batch = batch) +) + +## Reference outputs ------------------------------------------------------ + +## (1) cal_score on dense matrix, default idf/iae = "prob", logtf +score_mat_dense <- cal_score( + data = counts_dense, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = as.character(label)), + par.iae = list(label = as.character(label)) +) + +## (2) cal_score on SummarizedExperiment (currently stores intermediates +## in metadata; post-refactor we will compare via return.intermediate) +se_score <- cal_score( + data = se_dense, + tf = "logtf", + idf = "prob", + iae = "prob", + par.idf = list(label = "Group"), + par.iae = list(label = "Group") +) + +## (3) top_markers with use.glm = TRUE (gaussian default) +tm_glm <- top_markers( + data = se_score, + label = "Group", + n = 10L, + use.glm = TRUE, + slot = "score" +) + +## (4) top_markers with use.glm = TRUE + batch +tm_glm_batch <- top_markers( + data = se_score, + label = "Group", + n = 10L, + use.glm = TRUE, + batch = "Batch", + slot = "score" +) + +## (5) top_markers with use.glm = FALSE, method = "mean" +tm_abs_mean <- top_markers( + data = se_score, + label = "Group", + n = 10L, + use.glm = FALSE, + method = "mean", + slot = "score" +) + +## (6) top_markers with use.glm = FALSE, method = "median" +tm_abs_median <- top_markers( + data = se_score, + label = "Group", + n = 10L, + use.glm = FALSE, + method = "median", + slot = "score" +) + +## Assemble & save -------------------------------------------------------- +snap <- list( + inputs = list( + counts_dense = counts_dense, + counts_sparse = counts_sparse, + label = label, + batch = batch + ), + cal_score = list( + matrix_out = score_mat_dense, # list(score = ..., tf, idf, iae) + se_assay = as.matrix(assay(se_score, "score")), + se_tf = as.matrix(metadata(se_score)$tf), + se_idf = as.matrix(metadata(se_score)$idf), + se_iae = as.matrix(metadata(se_score)$iae) + ), + top_markers = list( + glm = tm_glm, + glm_batch = tm_glm_batch, + abs_mean = tm_abs_mean, + abs_median = tm_abs_median + ), + meta = list( + captured_at = Sys.time(), + R_version = R.version.string, + smartid_ver = as.character(packageVersion("smartid")) + ) +) + +out <- "tests/testthat/testdata/legacy_scores.rds" +saveRDS(snap, file = out, version = 2L) + +message("Legacy snapshot written to: ", out) +message(" cal_score$score dim: ", paste(dim(snap$cal_score$matrix_out$score), collapse = " x ")) +message(" SE assay 'score' dim: ", paste(dim(snap$cal_score$se_assay), collapse = " x ")) +message(" top_markers$glm rows: ", nrow(snap$top_markers$glm)) diff --git a/tests/testthat/testdata/legacy_scores.rds b/tests/testthat/testdata/legacy_scores.rds new file mode 100644 index 0000000000000000000000000000000000000000..e6c8075d9441954b047828f6c9b4f4d2b8d2f957 GIT binary patch literal 44717 zcmafZbyO5y6s9yt4I!O^gp`Cd3?Ze2f*_J2EirV%NP~1J&47R?(%ndRHv z!_4e&chByhyYH`e&UfCu-#zcWbCW%e_|gAbc*kf0v-gvWM2RA{1snyRzlIh3POty{ zPZU5Vh>xdCqWr2=PIBM)$Eh#@t(x{Lf@5JojyhHMBPccmA`*CZYJO^efsn@rpqCbp z0%U0R+V{qCv&8M~7R?{_-f$?-jegED1lo>{LR~|Ow_4rz{Y(^q+%3X;{Y~um=@cgrzI1=%BrGE_aQL%_Z zoYJm*Oq%gvtTxM0U9m{eGkBmk!C$;7;jULf<>PJ$NBp^zxSpVoGrQ=mxF?i2BE|i) zcduahU9ZNcd4WbthKHKVY6gujKvhMNG#0y%CDFbyc95WGQ%1_*MsM6O{8{yl0K&dl z(Z-b&J|6!B>SI?PFBZhWNTwacP;Pv;6k>VeAH8q~`mR90x~AeXy>fRv-G0%28|!nMkPlf~G99qp+2$6ns4TPj7rI zqkzu@x>R&h5brOk;AF~z_-+6h`ZIKRwTdgDD`kM7@njYM2J`7OI>jr_1t`vb%)t0Z z(Fj6j+A~ox7IbbSc1V@!V23egwpzgPlqqMc_LQJwU-(fiU42B7PvIvDa)vSlWkJ8q zQz4KHjvY1wqw*BQ_(qoY#gZ?YY>VB8+AHHFZbcEh19Mgb8voM!xBf9Rk2$SrI_+o% zusbULtMpvDCo@m5x#db=KQ>D=2Uc5WD1B)En@#Z#3BTr*O&`GtTnvc1{~Z%xNN+8U+#ZvDZRb-bx(Rr@om$u8KCU z12A`<^&kMj@?s{lXJ+ign$6VD6TdE~Jii(qchS8fv)-OK&*%T->HWP@Gm9JuUWIn?xA*F*@fGrgR#NV%}e_w^$=mae|}#%j9c zuJ`G+{g1&_ieuwXyJ!fdXL)!_UXZr!g@P z#Zva-fgq7m_6ZVEH3*7Ez05_Lave_^^&8J>MI1UjBO>Fr-b4H&=!->m#~{($Ud=$N zU7th=DKS(nKJH(|Ok458&)#3>xo{9t>@-O!H*UGR3W~XxekV3JgVjrZqTX4HD#z{J ztM*G5-_Y?A*cS-9fuV9zS8drwZ03=upE#7G5u+*Fy~2fP?oYIhLR)!D{dcEQyGuiw zN8Hvd;?DluIsH(%q`671P3Bk;IwG7j$(hw8^AG1YLt`g>iOVfX8Uj?_lh@#xdT*7?$u+LtfiE3pREW z$~tS@)^`}1HL!v&E?K@7U!O6M!(4sRlW1N(4VIznBN;&?)h|rYq^xsvGGNLyg zTcWH(zi=3czSE8)H6WBzukK``lY*W1I`?v2zZ4H~HmyGD*9=s$TQ-LF9Oi#T+yYMw zFvC8XCqI(@DRvTboiqOFs)k1yu)i4XpTX}L@)>OZ0dH{PZW+7&@AIgeoYjbWmLI3A zgh@4tJ2;zVyd(Ncbi5*-njm^a2%P|D6V-r|Xeak<0~q87sEh~65D&4eC^7}nrlgO+ zvC9@wf3ESDV~3Z?vgxy+ul^N_-vV*ck}Bd( zP$>hrs=UU*R(p@^I+y&3jv%5*fdkS0r#<&9OC2McjWI5~KI2>6f1>U6#*kUTG|me0 zz`Pv(q!yO-%$^BvV=$B4qLnQi`G2hbU&{Yul|T^3L`wIP zPVkGM@!WrQ)O>np-Tiam{{Y@}f`|VB|Iw}c^v=Hf=R}xV9zmQeDc%1DBrw+MT_k4m zbua3~7rU{#FnALc24!WDf`var^D;?456_8^EF5Ztf59m>BWl42FRhAco!t zD;84G%e%}M4{6SuR(A$j+wkE$DZ!*^i{0<}oN+A`VGJ0qh7``|X`U!o#*q)-b3SdD_V#OWk9nqCcPD=&@sSFzgkI(pH~0x^Ku!MZ4uSKJ9syJC~b9v0+43 zCjP8(ji3453=R`ASZ{!fK)a$ZQWajj7<(G;9Ejz{LTo%u2C>5w>OB{Fo%ltQH7)~b z!-AAh*pejtfEj9G$BLc3mz~M&UG1DsJOZ15GQClHjaOYYPP-v&Cl8l;h*M2Z1tW|m zrXSlM5I2wQYN!E!0eV6&25vRciPxg%$FCV?*SoG`v5ghj?GG=U;|;(`5~p;ZR}^M~ zM&vWlY57=@K{+;xa!9xQVEn)&gJZ}u!1-R*^w>dBLrcm42J4g1Jv->I>%S|*Hc6q* z707MwmhF3fAIgaZJ(%v^VIV7`!(F`-O@5w=jU7uN?86y@P|73;*uyM~920x~^aOZC zWAh-RrgEd<0u(fljq^VR?KWfIz_ON;Z1y3u>o<_y74XVr;`ZXhtRJehgA8$h)1Jb?{!QX%N)f@3X88aq!(K2UNm-H&ufH7hFVLw)5x1* znU|NAX!8>ieytTHB#d#yw`*i@Kdi}%&nyif3f0@9tkJ~f!3mOz2)z)5rq6ggFWX6w z*oySH=B0*79oDfO4SHAt7lGUF(2&nZV45}-N6PSov;;dB#zAhqn1Abgxk3HYBgmX) zjHYDbUa1wDzWf1Pm+yILfc4SA^0r}1Gp8H|S&cBTdnxcN{gS9uHp86*{K~k6tD**L z3jS7zP1$nC<^<1b+CTL>V)03!dAfm0HD)(S{dWk2-yIhnSULw@?E1=^CmCTJtS57~ zI&vQdQ*NUDP}^NHw@fa_HVmXkn`H>MqK0T?Omr3cAzJS<}U$C?N0^<_{&NQ z%iPz+VwPqhe)pd7nAH;C)dThcpStNO`lBA$d0&qK&0h%ou&*cALAn2WqY!WFBLb2m zbH$sEp`B$9P}D=6lXC`G`+b;V_q?ak1lJx+65P+e0!Qt)7D?sng<}|I+4b&Ki-7Z2 zrCFFM#NhN3Dc#urxoMe(;NPoGvfn0SrbH95lUB7`V9NB5RgVumI-eKI`4~DJmrX4h zauUAuF^o}8Fn{v-Etb9Id0*={e;SPaGE0}sv&}B{eeB-t2nj>>zxx4!TM&#GOM(@o zG|i_uo!&Czk4PCUkl7D?MGx_aNkby3EDiQOad0Ij=j6&iqD$Yz5 zX!j8`;@lW%Wb#rFz@?XXMt%-aiNpMZ~VnWD$#=^R*i?zLBNWUH0$} zT>n7AwP!u6aWrYQlXNpK$tH^W-epd&?lf4`D8=GmWPWyW69<;7bVDoXUa;#Th)bKv z*x0o1e_w3PGgn{}qNCHg6i1MA4yfCzy=X_SQh1|f3?79pdIE=+t#f)p0K14#WMh}? zh=I@bGB5fCP~X34XW0^Df!=1gmvzn43H(|hn$_a-JI07e(D5KAQqm$kZ$jBVpr-!UrR{7>$ z3HKE$63H)lVCjHtx)CC!*B>0bnlVhgykNp4*jQIu951*YvE0Q9ZVtkDfH9z*tpqvU zn;u7&BqN~lgY((D1@`OKxJ}@pH=i7lA-n&Q$l2K}axYo#NUYYSkK{`Io@(U!HLEk;_Q0bD zGb)0Y?u+uO)Z5PB)gZ>}hbs~IQBi*B`NCq?8f;^8BYEPFrMin$yi_TMQ3!mDN$zca zSVRcwK6uPY*IY*!UwG6@aZ)^7AUE$K{XD9rMR%n87o~4{yKpl1c9PhEV7Xe8C5_P4+Wn*~R1&H(%ErJVE9CbZ36tLqgtQLvr$o|-2`&n^# zLWPtG>R+9nE;=+CAp&p9ox$Fqn=YEca=_XIE{Z*D(G5GsM;^1Im`r1-doejZ(Mu+S zUFcMoW3L%zw3X%UdPF(cbxl-#?GEidV2hmsj&3@fN&9Iq6LmXL<^-{1!;w9mH6GP_ zY*-&~k5%qxSnk3HXPLyKT)ab|bsjaAf>M$r%>VDsJlOXk<>$G4y35bGH zN;hKy%o$67lQqFA<)FlyAuHf52(<&7U@ppATt0`FkR}KwB~<-SQ3$HiJNvntSSgHj z9?*J5CegEJ2sB3jeh_MG+#=T^7V`2@&bYlTXj6T1o9^@L+-B^*Y#qAON}`UMv0R~Y zijOSSJP)^R+wt4(@jVeRY_J~&p-O0izrg+v3H}EY9AUBtojAIjNtO$+4bVq@5XBge zBdy250=f_3hDYEw3XKvhr9W0fJdGiHlWl>lsYL8DW5l6;P z+VL>?JXb}MLW^RS>!zyHvohBh!Th?+96h?{KohoaKb|XHo=8ZIIUHwuHb`}b>!DG5zb%Ns1y-l-{+ArW88Tshx8v!Z;$?~W8r@0X zs636(^Vg7%{I}~}9uR-8ey%{b`&gbeS%zO0u0P`+BlB&;ds zV(#RBCXdXQJ#EXHA+tFU?7EAq#PiOc&Vb~v6Ju+24p*9p|5B|qgrzf8E0fM}<5_r( zjYMl#=Zn^XnDE&x`MoZd6wdPoyi2g^ zzm%o3-dp-AFsC+euQ%pdQm&sH-}e6qRtn50a`KQS9{Mo9N&=^)>WpnOrFedZBy-Kt z@`tfa{u@ynDbbkBgHunDz1K&UK6xnP&gk9Umqdv1jRVLU0(5Hb2p?Imy>lblfXmak@u^W`0+fFnt=i;1olLm3L#;1gYdiwc06-Av^h?GH0Tn zNs=|*4+lxZa)cD*1w{F{+n;BSk4^7ynbJ0vaHh04in*xjO8;9Lp+v*k zl~;P>J7A7_hSwGVxQjtBOBoU>b_m7rxCO!TVx(yRTu*s(w8E+LnWu7t@er;?X)xfv zz&q#5iBp8gj=_TXq0JNN9`?nvvQCCn?NteR6&+5b?Z(>#Ejj%ETwum^8mHE0$l5NC zr{n{7)$A9iGp#V};Hf)Rm235mfux-Oa}BjeC1aO!CvndKMjd?{aFCZ`5^w(U4D69o zL+6I8KWPIr>{VSAbS;l_C&KLAdUzJV88BdJ7Qp;kSit#7q50P-YTggfeZQe`c{z76 zDe^xbF6M6jb!YuZ{)n+{7>n=+$h56`cK(trmQ{O0@}(|8uGqqL8S&@#jCAJO+3uI% z1#J^cR+nY7#}@$=jg|R#m(06Qy}a(lfH;*9z(NzBGIq#Gb5SJd?o;H`V%sN2OLjlF zcI4#@M!p=sL&eJHwMm?gs6|K?3I~_)os-I=^Dh)K4*UK!diIs7m~nCvOr z5TVBze-V+BJ_pzlHF*9vx_dm-Y2LmD27vAAXTLX(Iw!4>J~i=Jq2oKV+b<2S4ckYe znU96n!ipEeUDYa?yIbHEA!qG2#B5SK-;j%f#cacG1?Hi_#6d1U@X7Ph1SD4V$&IJ8 z05x)%6K*gWuKhiu#fKIG$Eg&d)eagLkW$t6#YN`=1z{Gq<_Ff=K3JYOJz40vDLwt` zq2otQE;6BMF5ds@z;BP5nq+c1XZ)=~`|7`(PkF$DKZ}1YE)4<&HS%oWzCAm#%>60O z{1ME>Ps1+sel^-0%qg8AkQsMIR~G%vu`00a7snBZSa||foF)Oi{38uQwuC}vx~ z&uJH~vlx|N{0I)O;$Jtd!qD~rJF}#JtFeRsGUfG}C@PBT3uk z@f1@W>(-+`E>4p;;1wb5J!NkImzQ5s{UT{2&@xyRS1r*9I++0m!8U>eZ>L@`kOxO+ zebO3Uo&osLwApSmOUfZ70rJmPk6ew-sx9STQ)*YQp*hIKXTupt!#zQ9AoCjfXUT91 z4ki$OwMfq7$f{&h;oqT6mR~>IpV#*a|lk-u#2bOyE0CB*~@j=Khk_MYn2UBch69M#^$nh{soUeFCu^F zm|T*juJgWShKF6+Y~Bw|qGLs0Z9RIR^C2=}?UbdwdFf_!eGkLNOWO*kqd zW|_U6_XKx^OGdeS;ckOE$Gt;|?&9bziYg=T%Zp8+GKd@oBzLK(TyA5VQF|E7t*{?> zRc?Y%00dYxy)V$TcN|777=*W*pr4+FD z_D@Y6!SjrEISDUmcD~1yknL+{tF#>GMYkU2nu^1-0-P+@!Vo6O%jn!DND;b-+Wi2r z!0Wv0GxTNE#j&j@^zQpDaP_^XyT`fw%LbxtAIEzy-D-8Kz&4ch#jI7%OC-Wv z?uC`D1KlgNKBvaZps(jnmID=Y;GQU{F>vH=&MD83>6g=uoG@etz*o!7{WiNz>Bkxe zUzzo{V)2JPW7MtbHFN2smLP+iI}%&s0|_r#$a;Kw@BQQ^JWw{8xQhxaBZs*Bp?ClE z>W;%7&pU@`3z)Wxrt4Z+j6T=Wk-e?3Xk?A)1P{EaH0wKuqy$A_kG(If<2QRvlQ(@f z@8wb)i1WP5gA7-E>Q}*N`=Qscf~&7*6n&#FSt8;Y`9(>DHerX zlZWwL=j50V#)U0T9h>5+SR{KLZ2UatAR;!VcN4E$!snaNZymIb;6SdLGDa`BYQn61 zDh2x3x~)#{cI?OMXk*-iD}HpDQOr~kc1|h;60FJ#I|oy7<9@8_IGuH1=h)t6jD*IK z9ia64$LT3u!uiZc2}F%j_mK)e&XOk*R)80jF4gBZuZWp5n|GqtUkz2xE6$oMQ}VOI z=#U*h&xks9Zq@{Yh&YR{M)EXH!|FyZo7dVh@N6r)e|@vCE%382ZL{bn;!}Rw zbK60R8*tMe{i#%g64X=`^_yA5S8}}YtAbX6gn-?XK2~}ssUH53ptOm6IOP6Lq;!~l_7dXowRPFBn z&7X%x#-IlBt!m;d0pq z+U(VH-sF^-!>pfrkZ)$V&8}-8l4MmVp&~hJqQ_F7&-if<62;Sl>X{ZRsNz?4s?onX zDm(RO+i)+*8d?-g;2h6#nof4YEr?M;^RPjCA^rTror3u$W9AsjiMI>pT3^d~{>c5> z{P6WgOZa&AvCf?EOM!&-XJcQ-(vu&;yc+MmlG_`d3U-YK%Kt+QPO|n?Bz={xGqdmA z71kskN2fEjP{muMK_LW;H9z}r6IXE;8K1h6Ah-)1K87D8;(g$Jg`NBRKnvZ-{|n_d{nglq$a?x$k$#r}8_bD`;Qf z49NxuGj-`NZwH7t`=Oj-sgN+sVx)VHP{ahF#qQRWuogW}-N|%@*Mz?nSgS}D#>?76$Q3w!pC`w;QoYBk{N9SRXRifu1heq*i%`6)1XQ9tvgol+`{%G zr_`C?2W!tL-Er_McJK)GvnGgA@qBGl(u9walF1*UkEe=nqrLphg8Dy!;j0F6LSr1; z-G!LqC%l}xlkt!yDS$pXJuBC9N3Q0e`#U?ZRHFKR^A3&gdSD<~Sc!OD_oNPZ;Q{j>z(|6{% zZp=M{S8R_Lp)Ss~dM)LGLXI69QP*NV!$0^tUHS-SS>pa-?m{rV{)?+-TPwmzv!-FS z^Re|}0#fcC*dhzXS{QGFr6XWl#?3ljtU7jRWIMdwl9e(tME_K3gEpz3;kC@((c z2qTS_EI%I=BvZS7c-@)p=EH7{R-3^)h1|Op8Ii310v|ld*0EuNVM%STCDcBi13QBCZ`Nw?DdlHB#3XE9S)6Xdt(GDLD57~Ivkc;i4Vd9-!@*utmmioHI53lo~8c?=MgpGna=(~Fn z3mwI&2bwaUhvRsO{<3g-qwev)&CqfyJPf*q8d`*Ict_lX$j~2B?s$A>Lo7JqAj%e8 z@%-@IXXsgR^Bq-w*SPk_N@@+!EN0UlAG-P)|=;LZPX~Mjk zXy7hg#5i|zErEnk_|jE6zrmHMlC1)_L@|-w7>UMAxLLn^^Gs3qQNsEmuZY=KZV7&N6wRcJ`bGTxtzC)+acz$88vac(ovzS#heh(FgY@<%!p zaKYINotS()bOMnJU>{aOO4l7w=`#U5HiOPvFA&Rzwy|sN0`9nZdYbfQ_B#{Pp~b3V zJKx?M7J8v`T`Uvtlyp-ZI_!5r#b%j~e`4~dafYbogz^nJANBcsOW(Z!#As#*MLo|Z zA<^%-p!Ibx4pM&Y|LOBSC-K}H(Tb_P$o5L!wF-A~&5+~kCDXvhY%Z_$pIWdVrlA%)K)PKoBzj?gbLO*i)rAsP#)nj5)Cd_fIePTf;`5e*fdlcKLF)tv#4% zje{Ug68|cqEu`o|d=eV*Xf)}{d+U&dQ(4)3NE6kMhDwl7Un_V9Zp;1KQ!?w)DVB~1 z3Z_$cQbrSx%yI+#&hx3@mmFmRzpCCXLeO?RSzfLh*+ymC_ir5u{EsYc$;Vk z%BioVbQozkQE=x4V3$#Khl@K=nk>vmK&}HVf62&g&;{}MIfRQVgTQZXM0F<5cNLI4 zkD(}b`nhGFkr#!ueE60G5>b7yzZBEq%-}NkPl!`bYDNEEV$!R8>%0xAu*;_Un^S(IE}9=JpALRkh}{Ghq2+FK9!z%hqp)KOV6B0k4rbJD z!3o0+gtRxt&XgM*eju}R9%_K%{NlN?#hUYsqdU=UKx~SkwxywsOR?bSA;e0!)cbyA z!%%GbH%C8}7Ny6x;VoK66+`-TG$aKYD&b5B@R3MVx}suTW5A{-a@L7vdN)`JN;8jZlz62xH$F_}7#1 zz(?Ok7NtH2qlr105)5btS;K}rGx#XOr-Kvr=QKl*QYFUWRi_ z8vwO}2>~lzTNEF+0Q!nT_5WT^o;`g%m~at*^*4tK(QTTf=};rhsP`(y@PvFEE-2*0xI zG~4P|iFso2kCi0gt$tcD7w>nR;fBf65nGaQlsN$>&zx?(d+VcC?4MgFiN@g^1!CmVaw&cuJ(QurFn+z5d9Z z6f4aFj2ZZo>bulo3^9Jxmw(`4sT$8wGYmgT#dB!4Sqhpo&QQ- zV_9TLPZ%UO%$J{kUAwV~aW3Wh4V(%$4cFVp~@`#5l)OTzu zk7Kh1z$N{GQ~!p?wTo8574!OFKa6-wsp4SoJ#i((*isHg@pL4ONn5vqu*rJQ()7y_ zUPA2#8^1qKgy!j_+eq-E!eV%w0XT~*NCnpEHQZ^1$M#Z(&0az43Vm5}1IYn$-3208 zM!1RSl(v#%xpXi@+=*R0j`GcbTB2u_MsiQWL8;++%xX9{kkvr!YY$SJXRZDs+FgX{ zffjulbP~i1+$VaD3PgQ%cLHY+0DjK>#_P$KQvWLCOUfM)k}(-xbgY0Gf~{uGs=Uh+ zfZu`ByE7HvbKmOrg$!o`+BOlIHgji=SEI_tHqdCmzcO=xa#x=@0ZeF@15zVPgG7{ z<-=*M9``B8USy?5t+b3*m6$Ozg#+&Y)#chgZer%Z*dWQ=5f+*+ux@(hwu})`9~%Zr z_>`{d4!aw@42$J$UAHX4pXoTG={-G~A_`} ze>kwDP=R%CsN3ghUeY|44LszrvB2AE1iHSUxQ0;AHnE)?x+(8(BE?nq^u^#@t+a)H zG}r!gO|;xzAGm|vjBNTKDIT4iU8Rm@A!9K>)oK}?4$-nGA&hp&@}8F?%&OP|v&x*e ztk28YU07zvw@TzZZ%*HtsdUSZd;2Yihp)SaL z{oe09v*S}2FWWKH&lGy>O$+wK2oi#1g`Zb(1YfMlwL=x39xPo;ae7?fGXL@12w2+r zvMFrah%`(dHLc)?Ga`m4uA8pcJo>p1_dajIjPinxj`B-q;hwSW$B#g3F$<2X^ZU5I zs@y-%pQZ&^^iaw4=NWdqnGSh4Lz71cKkw2a=3Ek+NJ_q%QlGryZ?tv$T^Q|@m`KrA zDI?U>Cr#MWB^nfh|0-=P|2wEP8O7X+&5xX>fVh_e1Y`A|| zuO^sj;H}IZXd`|cr7QJJ_ubcjkw50&k^TX+dVHG9eLV+b>g)=gn%nxoWM2W6f6b(Dsa5f2#<^l4EhV|_7-nZHj^bdm+89LDFXVICyE8)=`(a(MD zaz~g7EU_g?=(^@-k$k*IePa`|_6Ed{EN|(OM5*_*y#P0a=&#>Uy1BASbSq+S%2)v7p)v zBzbE0ieYN8r?w3%|9(7$B#X2R zbt^l=v&lbs4^%U!=;PU$N5A~BfoMN&dBu6iP5^5LDz&2-aB+42HU=9JeYs8?Y-rH? zDyw5OZTd+bne)S2NJf~53J>JEw>lS*e6l8fDajH_38}Eg6qSi#&YqU7uy;EW#7?Gj1&hh#+$dMr{TYgLQORg#euHSDV+B- z9Mf;={!$pLtS;xafVv=R$Y&1MF%k)F(Nnz__L>j#wd=Rkha%~}mLM#xd;>ASBIF?Z zL~Pz){F?1&klR7Nr+XTE8qNAYH;7WGzH`yFwn#cr8#D zi?NQc>!l2d|Aj(^TSj#Q&xi`+7LIeY+IaDb=St|`HyeeGuEyCfLW@;?4{XLveZ%Dw zzB||8sHWjhN&u#n>5||}a;wb9s&zZibD!e%-4TH9J`k(aFEzh zI^<76($Sdr`*dv&u-nVfY@b9wiz7yAu&(oO$LCs{HQ#}$y8b-a$E>AP%T!(FH=U7v zyyop1=F^^ki;L;a=?=dT=Q$(8n_f&);&zd?$985u;$*93lPnuo%I_72_m$xi44ZC7d+QfWQ0>w$T%c1ods}c{?;JoDU6geyQEK z=hO+yrsE(&%-`3q^qwG9RycL8W3~4EU8cF)66o}r>V~=Q{0Mk&=#@7Fg_|B7g#+e` zJ_i_8vOe8$q5#X=k)#A4`Wn5?4@>@d!)hqDFXPzm_B-9fM)sB0mpPayJnTIBk41Mx z&hK`RR5;=%M5$-6@W3RHzHa#XG9gX!I5rFyT)V1b@iaC69dpJ{BL_m8XOQ{megaP2 zyeiSL-xN=8UxF=XcgFs!=3P>}qjZFXRo2{1cKZ%l;Bi(cJ-t064!I@wfM9c%SVI#R zt+0%R3aK5sRJsLw(qvc}>)D;EU{HIgUundokegf%mgO(|?b1p*^UQ~1p0>vv7g~)s zVh?-84&w{Gv!yC8u+#iBry}@OKRhesEZ-<;BIt*=S8kI#RqNlLpaS~96P~R%%}{G{ zTJ2Gj{P_ZsmOFj>zrk_^BjCS7gjDqky%_Ord}JkEFpM%wZn7SAH5yh`td!7jqspO{ z2y*ke!ddYPPbq~dQEbilz)d>&pu#CQIK0|3T?nqLBS3f53!7FRO(4DXX zH`cXKhz$;v*7tW}GXoo6RNB?N<05|d-?o9#$wg;@K3N|@dIX14wU0w9)IPO`BySk7 zl5(EBw0Mu@eXCee`*0)tL>B#}%l^juJo5|s#G{DWk7P-t778OwkIwP2@$LBPWv)6e zWOb5!W8pKrcZ-uODRW78U?AYfD7bWg6wHaNyva6l`$WMFlA8Rp^&KBZD;ffSvj7lq zcS{{O(ecMH>Gj5K& zA1&e-!t&Kw6h1chjoC44r8}d1dUcMWSAiT|eFryopMq07J{r{+oUW=)3jFkh=3qVyR_E?a4;_O4}0WvxP~%~FHhpBWGC20 zKf2!QM8=bL3$O_$gi~idI$FA8|0=(EaJxp7_tinv;r&&Z;yQEMsNhdb0dmkA)jF^9 zsOu;V%r2_|+tZO4rOJ+5rHp)=^4|zo6APzb8uK6f9Rh9zf6R0e`Lx*EDb?@ts5`dQ z{E+$S+a0s3RpAzO(0WPUtcz|{AZO0QR@5E6x4-M&colp zdEKjr%+{4b;jMh0ZJEq*Gw&06n|v?9UH({%XZ*4w=ebJ)!tu(r`Yky?+R8rmDM?ZIEj0%gdrMDz}CevBOb_7)Q%K6bE z7m2v@tzSB{hLKyUnj}2Zx<|j}UW>l}s1acK*oNP@)1Wzs`Z-Oo0xl~Pr`)(UN+|%y z3@M(Md$224LPxdBpS`Mje5Gi$RX@+ZyS_&L{^Drnmhwi`Z;$Ka_}#u@a_dcZ<9&F4 z@|!${V)+nO$WW_A$MG8-xa@ z(0kQbvI#x=Og(hWDQLDop&X`16V4ZW`Ag;2d!uCj<|=LKF*oGr``;*6=0i*5*`K&e z(%X_GdWDS7lwKJ3cYg%sil2~#eoMo>8I4gHF#;!aTthhuE9!6GM6gDD;J$zDkP&_~ zI&i`i zjE3P7WjE#>HKDD3gv}XJXYVjNw?r`tNkMdDUmqmS!bh!tDM#rA ze_!^!QUot$f4E5Pkxv#I*8DUn_^Y|qq)7Rqmg;9b7$7jHul!CitBZjKz7Pwi*g0nU*sw_2rEjHKwsfX0td)72!o zVZJ--AszTbx{IYjAJ)leNy*!IsdtsG-t6qu{S9lj72ey6JBwI-iBs}@5u~(EQj=hW zqI>j5xF+Jm7aUnv=hme0xmlvjpMO=J0ghmQ-ghJlR5z8Jm1_kUi&K zmyteQT`>6OYexA|!a}SU7LyP^J+@ALKK3}7uJe|yC+v-0{=FY{Z34zmarN5(Y50Ar zu7!ma<>qM@4}U&mbTK6atFm9OCe`(_28@ky$Gn*~V9+;dq@~^o{ze>>w{g#bY@e&p zCQb@U`{{c?h!3$eX+|5cz84O-GJ54g9pN8TDR(nRI{IIT_6(DPaCTDsZJ{0Gs(To{ z-_PhO1tt4e1YSN)0P%*MNrZhgh>($u!l@@H_>-sVEw2Qm0zhG_UI8Tz zR45t!YIqpE$VvmJ8|ci!NlQeNLjK?Z=`Ap2gi0cnFvsgl)Jybj_vCQ0x@ZwBk-t}z z&_S?0Ba=6BCf}q^kBaCwgcT!r@C-3xelHu82hQ>hyA5e2WfWU02ph%viSv3Y@`HjL z6(Tg-jGQ~3lPh9O_~-weC{fF=Rd}Z;f^(iPD=u3`CYGL3oQj<+CmWaCZrQdIZJS?- znad;*)jF=5EtF0pT|^ffSO>53aZB3HmI{C24w=xR-m3c0akrl~-tUXyP)7}ziw0g0 zIfz~=n&TS2-7mKihJ9P*4U2hs%VorSp>oa*r*8Z5rr@o@0(op2hv|@8(=K<2W?pje zIbj#8DawJnLb;(M?JMo2vm2)+34YgGNeaaxg8GGH|8S6xO6^|xQK?0N`^i;LJ#fAs z(G4HGdOYgi4Pa>ZjXVC~4}i4xw!G+Mr`Dg8U46MkDf;|Nlp^C`XzhHc(Ss{HGWGU7 zkDIsBW$6a=bK{6ox$lRM+_`uy%y_GyO6H%Fa>^=AzY6y@E8h{c9S)$aeu5U|)_>q7 zSu2P^N4=8020I|HXF}*yEKJW`6kVhWF6KlOJLvPe+U+l|SM}7+itsDI)MwHffD^rC zezeffFUdQg%>nCWJQ7R=tg}QhEWOY)CDv?NY3DG+?Lag*zk~ggqSFdBN2}vWUbO|) zXCtgJTT>wo^$)Yyr2xD;3`tMMmFMDT?3a_7m*y7!%4n5A8S%;Du~*0l2FNM#!B$upw+_@- zVh%fc73gneje8%QN;@6KcoIq9opsAOPMnOg+&I8|{&N919HcERAVi>sC1%}Z%hSnBO@u- z86#UBos1?RZ&6bh|IQ(Gsx<#b@;b3Xo#0oF(^i_!GFF!u1KD1%-4VS!R_2o43MDtL z1Go!5yiCrx@~-CMOKhK`jhShxObyE6beg@(&NnXutkeuW%#wYKSq_1f)o(tt(Oz2- zucm0YeX-`LWkEs7(FU9fm!|9(R9!B@B7O=|CTGuAzU;SA8OuLAVLoM6J)bAamX$qy z>xR9F^jxdrjDVa?N=8=!CMOLI(s+rqC}d124WY_RTBysNk5~RPf4xtSf=vKSLcCu5 z#dOr)O!smG&4{)QX#eY*LOF}u(2Kp{Wq>vxXp_ip*Yb*9YJiv&x|Jeg*ktkB0#7(% z!G&$o4Cm=D#I~`6A+G<_jYb(SW0(=&Hv+#M)fyWeCHpHD&GwGNZlVjJ4aH~21Zp~X znrHY($?qLuFQ5&T?czlTs3xL&oEv<-mK)%N3v;#v1kQX3WxhSKQOlI8(luMi&F5p`aY*GnZ7-M z|0PdT^X0nOJ6?qRfUYd7(k5-BgXZ^}uT$DCap z3B@jvT%|5{ZGMBg+^1NIA^RN^uRM&f4>mue2u%*u*OXCzSPs02j6bnLaasjrx&H2y zd~2>Y>ct*kisFp^7S5v(o^H94voLp~dZ~1(q{@^28?nnj$AO_V15ph-Yoq9CM{N|=( zx-7EKxQp@wn~`$6B3L0UWvt!;t<=EV%H| zKR<;>+y(CH-3qEIY|{+-BP2) z`puy|-Q#D7)?a{A-E1291=9>B*EQb(81j`cyI#2-61#xcOBee@Bs!{ zwrXqlWncGUU(UG?S5DRK+kLD1@9$8wmLHdYrb9(5pmMnw4c+!@P7T2f?S32`I?nh_ zOkn1=uX*o16I(;wS(Mc&M3#4(@Ui3cJW$>a%AUYo*>hGBr_DmgM`<^aFgO9iRF0H$4 z@~XqRbDq=l17PzZo>)>T2QF4vsJ(>8uzdp2AJmSr2Y#Zo2TQ$Rn5Zr9_{RB!j(<+Xv0!d8yLi~`t>u6^$Ul&bygY8O;ma%?U<5jdm{54>kp zy~XkWoJY5+{8d&6jDGAW zunXgoK6TAIY0_aRwC2To5ZP5wurhc>xyg*xeoekkDZb^`rFWIdt2XrMEPqIi4u+-h zQz;9xf%NS=8Ok_8l5P7;Fxv&Q=ToBI@oux!FbT2<(;tevWPoqt{=3a7O}cDSaf|DH zR+IGOObTH+tQtjk8Q!NDsgoKeJE*r9P#$vzTwG4$y@Bq6pG8U{x-8oJfBV%yqp1(#>YB>6!- z*v(O81jG)uUskq9Bq}^SCS^P@<=0I#&ew6MJov4Lw5ty!48g8a6Rn=|URT%Ux!e*C z`x@hmlpu)`;o~A-$;UNs408t9L2Dcgg#MMy)DYB0B!8ITMp<sPas%aQM^Y(v`cP0^_Thu&rD zw)ucvxD3z)7!1$e@_K62qCijqHSNk_8<2{kgu}K@#-}*h1$MX=O#>29tOR9nb zln3O_6xi8nJloK0Vuqna*>TudE8GTRD*@|1zOw6ym6kS5-X{^PrsYt;ef5tqgdq^{ zS@(GFDsnTXi-agiFV6_)IvBQNMs`V*G+bcZXEQI=xc zjQ!~I3eFBm(iqW}0MFgEn7Rl0X+!JEc4;V}ZocX?=2^ewpv4|AAJ!0Gu8&(G2OIT? z0lzE^gRX9Vj!W=~X4XB%G56}4T#P_@&=j>VaeN@Of(d2D2eZ$(-YjIh|C;DA(*`P< zXF20h=6%dp4rBQPq>mzo>=SY6>NTXUw*BkZneD|d<4*0v(qskIIc~8A#~W6%J++V& z?fzU`e_TcRvXtH)Q-yt!6SZn+OO09pqbn54)gcZcKrF7iyi;X+qLw7o%73Ee1IBcay5k2!t z4f%Zku0ML9Lo)Meb`3rJGLHLEFZJ65|6ry=ALL*zJnqf2^=QUx>@&CY3yNM^=>P@W zVUAWXrfZAS0;e58iRzxYATGZ6vJ&eF?~_DTZ`dmY@--@Yz-TRxw`Cn|L8 zjKO>tzxD=n&&Ta-tgHP>#-pG6@H5^tZ(hbPN2!b zZeYp;d9W=MDHcUMWTBc|US~9l0p(U9={-qipw!-()n&}HE>bUX1nv9QWV+27p z|FijZ!L+>ZObA`7Q37^eHoZ zR1dmu?RHaT*liMeS0;|r(k^W>kylJKU+hs~&KME5hOXd^t!En}o7L&^Oud zq;a4Q?}P*xcAW>Uh%ctHiFK@8Y;0!pv$4FvRUyON`5T$mOzXsj5@am99Sn0dxv}5= z!q%Hixtj*8vpGhZSDy4uDvp$%o(8yK8o|t{% z1$={zrKmQ}PKy_UViJBd$;GzbVQ8{*uJn=M0z7R=440Pp2RId$KJgQ*0k8hD`cMj0 znaY)fHPda!tfi3ZTKks>R9hS#rGjrFXV}v^=-aQ=fdndM#}%9lubzr=S+Nbq0y|Yd zNu5)vi6UeVLtAa8^l86h%M6T^W9-gK=B6}8bt>T`s4dgczV$#|4J}kwq~rRkt1`TT&YP_ue50M(jlhxg++Xr2@U?S3D@7YO z&8^c>s5oUNMn)ZKkN>DB@B5F$GMHLvwClvW@b`k7`LD?tSew|3EQ7|Sd`xsh87{{= z_lGYn7h?pV-@RbSvoR*wzM8|P6?%=wiJ z!Mb=#gPX%5l?O1A)Eq0#6+v8RuSY1Y4at$;w?*KRFp_?^g@1F1dlxh2DX?if(Hq$+ z_R-NN&nqQhxeIM0w+1dDl;Cr7z9q1(LwWHn6Ll3>^6_SZ;L8)t-rn(%W$Zjh=p5MD zbDi_Ea}F}y{G2H6DrmS5xd?%LJ@Tv1um8xC`&YT-Bazj=AYq=|Zsqr;%)DvgrcpUK z|My&qxE$`fhDIGTX3d{iT4qwwX`juRd4Gl%{tIH}@t3&!zvN>4H-X}RM^pSaf#QEh zQ&2o}=R*g%h*NwoZEU^4ucxj~inVNTWX$frQwTf?PF{}X37}(&{lj1j-56m7>OFyX zKIv!s-O+)o1O?a4di-s5j-#PvW2-~{n@ZacOz$Y)c;1eNpcGSGR3v@Q@s&kTm}u}& zpJTMT4uC1jGn@+y|Ex{ z(=W2{lGHl5`!~+`CWvNC?vVSU`kH%Z!QtaD{$(FW=#Y70Y+q!z{|H8C{}u|%GEI}M zezbQ_dO7t@1zIr&hbEobl-BsW!W@_OI=TAs=2x=5j#<0muj3Y7d<^fW9p%Lte14j3 zMRJUTji(led2Tm{Oea`aprE)j1SdcBMkq(DkFkT}pU}R*@DcOaeMD=_)>&WjI5?F#>Cj`~=2d;kQ5an4in+U5i&SK&coI}H$h5j# z^Pu{6*{g|$a};d6T{XENkV3l~y-u+kI01gBf?4Sb<#db57etk3;dp}|LPr1x=ZDxu z_w)PGX|v};R_P_sI2ISlo-(xdc#+K5w$RH(3mjlj|3%MK&ra2yjW-$W5TjIh%}POD zD9XFxhPpMk^Zs`{u91*|0odEJf^8`IHPfI7Wc4JLORfL=<>KwuOkincZ5>yKZ)8U7 z#(*}SMve++$ao1f%5&lV>-}H?Yw$XQQ7lQ%2@aK&R5za z*pCP2Kb<4}7z1kkJ>f?=#+T>Y-{EhB0A|dq`fYR-V`ew$XE+V+;TidE89#|H%#sG? zD4C+yyS8Fg1}=IwdiqnpJnaYNAhP&rwiZ-Hn=_XN$U37Og_ay+F8gbh7=0GRC(k3Q z@ohXrMFp?V9E^9q-_kz2T&W8gq?H^J8rzmgq=sv#8*ae_Uz)n_e(Y?n|J{nz&z7VGwb28A zhj%%O_2{de&q2TRM?1`Se;Jk==9%l6tbhw1iR%kx-jX1MM}_h8Pj2@-NHbytDr=Jm z&IhWTgmzvIbc^F$^pGlY|Y?USnW}o$8$sOB-?hJE!iV zogcJyZ#uqu{-o7!{UBJivv9U~W@?SWKt27^UKZG2^kVDRv5r%3sK_9M3jgVw$bd=Z z;@Q9$cz=pWwUzeJ<7=apd+;r}FCRK72`!}&p0C)cxPM|!Y-r%5xTx6rHHwiC)q~A7 z+0F;gzth+a&yLr={A;J#%M}8cz&6K7pRGk~Bk>qRdO@=5QoNI?5?k?MGnk-T$jm@W z&+gWu2nkmMRS3Q|4Q1!b5ur@%UT3`fK_e|t++GTz8@<%(+<0V9Z6lljiCVejuodM} zkC8OK1=xSOckl=rcMzt0wt8mSgM!B;Znu(*djjKoZu%Q2BQc6aTaA)_gz&3qg{&mc zj_oD0(A$^x9mE2m#}aESdO_}XY;zJXV2)~gQE~7cj^|Hc>HKq1{yjjmZ#1c~w>-8d zHhSDa=>4byKKS}e+fN_zQ1FKanzfDyB8Voi2hqB5`$3zNvf@jBb-d@#M;1z>5Cc_5 z$9_r?uB|>x!Ezq(5#g9{#7vJqknB?2O>FNi%T`iMBPUL!bCVMp?7M-yP^K4_&7NfD znEod(&Y+!<1MhwCtY!|6wTDnP_C^*F_i606!;Q_22|RuZFV$|a@i>;l80q>Iu0lnm zY^YF5PXPm%dygD^(6omyxW}HqBiH5aaQI}cv=Wa)c}cW%x|VA zl6JkIb_kKWhb=GE@B6c-Bd7eqVjNU;khzHjb~OIij7v%LK4#;>##<;eapSyar-v;j zD*T!LF4IEQac$$4BSSP(kXIW|-aPXYTa5DaJaiu zc6jTq=PS_r=ZiPsjNWV9ylwvq%qVZqyg_++J3b%QqTFy{vZtrUtFwz}HrBCbi6zt5E z#gWtF`$~<5(yq4pO{(~`#@OakInqARPyc9f9(Qwj*b^@x!1^i zZBBOlXcf$RJ*vUvj_sY%#T$KnRlZ2Yf^>7w-S_xkpYnSnlLoqn_)i5ldTs^w!P`}! z?pA@oRoME;xWnnSV;j8ZU~PsMmWNxY38}#ida9&Jf2+bao6cH8hy`YZ_3$7;&Z-*TT(u zeJn?oKWUkc4dk2UC&-7Z*n^x->dGdks)uUtvu(V2V@mkX^Wjfh(Q=tL3|AJgvaKxA z_0`0^&onyX#ZH3ot?DOkniss{G#u${OsfjM{$UT<<{-a)Z}`BMYwhlJj5j}3z8JhK zGsEFs?BIqAE%IJAw|Y*G*hd?owihD@8aXexEiRH*P~B~`pPxED?>sh7?rv@V?XTrV z_)qpP3;95ja50Kut(zM_5l!Ns4Jqz3j^20()|E%pm?y`(kNqfH(bppmrr%(D--Q6yZ_;zk}+RG2hvUPlXF2MdUO z!8^=oG5qfQ^RX|2>pxl>W2I~BN#(?sJuue9AgWNMUW6e^rn6}fa1GAk*r5|1QydSe zd{vL4?9%`Rkcrz5cta9Y>;*ep&*q`mX;Qi7nV#q3kIO^@8CZ-`Fquca0ALyhSvez_$C$?YF0T@}+L~V=2DJ zZE2<7n_2(U!!hR=w|;mFA1+sbh7eZhj= zI+w>+5sBUQWaHBqK-JxLBvgz*J0H_KgFFi+Q3Wq$&iFd%+YeV~$G(iuKOSW#rUBfy z;CW&2K37lv2|xM9_ylO2<$$I`mIU+3e9GUx-uFS~YR2Un{LOLwqlHMZvpbAYMTRTN z36m!3iR{E=-K;KEvavXb>9z0Zi*NA^^9{Y8U%hXpH;o~fGD7G4B(BXFNGnZ?oeIgV7#cjlyjejNTBIW$kVU8 zSJQT}WvdcqMD`fn7kSm=@XWp4;NzsA(s>v-(!_^awDi~q5pq|XJ8%lzYOyk{;Up| zf7ef2Mh^KyK5`H7B)J?vYgc=!tNY|ld_}-a5ptyJl6N0&Yc*hPh!0(=@5C(vYh$oh-u$1M7>w|+A$T`x(8hg!-II+>O^WzUIK2p zxL2@>z*lM5`5g+2vEoKzuLR1$XSOS&2M$AGFmd`LVh`jR9;&+0-Eo>}J5eczJFmqX z*uW9vtHu}ew@WpkX(IdaY@;DGq{p42ZvMvJ=Yt*n75H}=Rn9Cv^%0JW60yKvxf~xc z@{Anm(5^Qd3id9)1QJtkJAQ~?HDw#T=W<0BFL7)q?WD%f(O?=E;5M7YNqA1O4q_}#zLxdKjb>zTj%7AKp9Zp*HgmmHzXDo5FP;g)IL)}uMd zf7MrDV|6{;Z%PzvzW|T3nYI_-@b{)`Oj+`JR@xoA|CAXRp-Gr>z!frHo~DBbR1U;v zpT{=HOO_7ZSYZm^>Rp%iu3X6GiZ}}&cpk3owE{!-C(28)lQE$c=7vh3(KRGn%a|z* z#=*1h3m&L4kK;c4vP1jG@lA`(W3+AKuC$tBe-^?cbkY~fm+Iu0Bi!~Sg&tEOI&_U1 z*FU7sGwcbS=DOZjSQ|ITB_&Mu+zfDkR0-iTVhDOgo66R?t0j(yvZ{vN{LZZlK|W&2 z?+o>e)mL>{I9g)ZI;>3}taUt6>s}tBn6ipoT~CB7K1RD5o^52DH!_@>m}K6-_ia@< zNQ2b80-7Q?Z}E;`=xst3C7$q!*ZUxT7PDUmxL5O72F>Ro*`3x6mW|=KQx$8KbBB}7 z`nhWV#ALaY-mXEGb-9|&nOqi?-vMi5V>c3R?>h?SLpJuufm&Sl{Kz9lK6ZUC1 ziN?A!8pq7$2GhQG-tTQ*ID^Ydh1GB1r(dy&3^gNoodGdIFS0uylI~D_pRw+YrGaDj zf4={JhUe!=-vQqc`(yLl{+H?f(X%zM2cBbTc9FAj{&(bqmjFA9&e8XP zB#`5euI~U>xq$mhMtLB?_q<(j>n+CjTaQCulUW{UOHBB3jd$Dns z*X7IAOE8Ia&R$@#1;`FYaxwb!F2YWdhAcqp1{hhus2^%@?6>dy8}Kd^*gP zQx;i(NXHPpYOlhL$c7AyY_g&9R~pY6&j$;cHi@myk(N4$=-M~;`uUy?^-Tl|> z1`F$%8)V%i6D*JPM!fdA;WLiV(IhiECG)#ckPCImbX3Q~)zzolR9p4=4%x5K6BWb- zGk&;wam?grdB)n}CIQz`X+G0*x3`M9kprjxFwqJ`QErL2bnr)E+m$7{?ak^^YHVwO zqbp*0iI*kpVEX49%02WMu#VIPF^j|L9A5i+CxQBrCCFNVhCn%E%Hv(h+YRaIwoHLB zRiXqiykTikEr1nEX(KN~DemD9bw&+|_@m!M;pIbqu(2J=wbUZR-ij0JXxK`=$DUvU zCxrUZvg;0M1aksI8`$NuvPQU2sWRFKNZa+t7fD}$8>hSaoCy4 z%Bs_6j{$RsW#wVS4H#XUQ=L=H^OP#p&+J#iUtsxqBXZ@^p)u^#+rLV_`Dhc3o4}P` zN7fq3FgOdW^=M}&RUfyG^mE=iAF_FPn0rq+GUE76Nv#5uGj1M_m5We?H+z#xO&R|a zuKarZ7R}0K`J&W~obJFXHNqAybDeWT=Rk!9xwu4yd;XEgy zv0Dv#YQM&HxJ_Kk86+-2q^hS4Ap}jkh^JwNu-jV~dApl;gB_eJQ zrF~K)g2E?(6By;on${o&Z?>G8>@O^3U2mS6e08q-;!HM+Hr`mh?&d-XwNhRJ2nIKH z{K>&lAVeJV3#8RwBpYJ^=Kp`u65CT$V_!(kHuh8*E4d>-x{rX8zC&DOJ(=FNIU--j$*vk6A$c}9=NHg zVzG2@y@#R~7fBCN1!g)j^t=4T2>DN8Td~`)3sN#Z{i^e5Wi_}_`!;veYz|)}F9hFz z(cJ6zWg9?C)HJ>9&>542cag(XBOr+-Lg^Y<ULm}AuN3$HyU>9Zy_o!(9}z_omc-tV2>Q`MvJp! zE_uxS;5yfw5LT-}OO`uh0>SC=vzB2F`M?;Pjnga(FBH%6whl9t6UPm|o*h^m)6;2Y z^Eb>OXiL*)DKm^F+Yzixn=+K_$3VT#_R+uT>2`G+l@^rqdcrR>9&NG>D_- zM-jF=-SycVblWxUTCZs6|iqbOPg zlwT*Nm>kJBX}OknMDS(SQF`6WeZYgjrq0R(*#X{~0&n9QWy#j-X35=;D z=t(o-8AcwZxnn9-49Zd@O?@u$5GF8U>km^ciMO0ZY}zt4#(gK93-a@};!dcUgoufE z)Itsp(!BX^y8|0^E~8bhR_r~=#h$3*U?RN(&upN0%ccn9QtBwt1%p5A?@OB0-Qt%U zGk?6=u^Ts|%)h5~F-y`3njx&$4PzLo3sFaCHDRF8rP!up%|I;sFOMtdkmYzd)(sC^R{whvXun?ZG!K@8s(;Gu4l}!3+8SF)V_E<<_xgd zp{hh4_G(HbR53b><8$-&YHn`Nox~ULEn}rQcE-v|48HG0W$ddo>&L*Ap+15-p`mZf zrjD3xYXY73R}^QXEEO$B4zF3QakDZ+P5>*?4cam0_`UAQW0pBHsPWA~9a{*hZl76r z`MD@xTr$FGHW8-x4W8F0Ua9=H#_4sQjk^QNj(h8L;{<*j=o*~AZcFp)?<8(3>|DZE zb8Bif*9&mGLMVyrtXXZecRg^)T9&; z**d?1R&$Il`Gs18F>3=ju+#}Hg$kk#?qn?eQ)88uFV3C5yD}@FnuFlcLZ^jfZliC< z7t>N3y{AZX$1#Q|OO{0mp3`_?jkmKM;z4RwQ_c3*H762Rk32?BJMHR;(J7rT1`8!W9JvG|*=!F>G?sP^K8y3M;A7#$%2v!WQiP+?n z9-Q-eI$odW{ALdD;t1S-I!Mgn@*OPdXQdMBJRB1+;vEjT)T^NVKz=TY zV4?)*!O5@7>XaFi1@?ge1CJ{4&z|m?a!jro?C%P`mfm zh6U@~r0?-+oecpun0U_PTtMUOLW!s&>-vWsr#N4}GoRB%_Vr!I??;hzZw5V2Tys7l zUIAFfSyzlq^i^}lrjf;oMF`oFAlldJJiK)<|!Ov5g6RaBDK)f%OP;ZEag^y zlQuW(D}Qj?&;=v*0q%1&!pt^AFad=+r=m`znOZ3WCu{9gtj@UnT37E-`YxB2uHT7E zU=cREUS{Do!3+r|YWKuYMZ`S8tv-A*tKjBZOFXn=1>{0q?lQX5rmwBtr(6!b$t}Cz zPMty?rT)e5x`krvy=Ycgv%Usf3cv%MpQfwstfm)_;omySD>F01mpHW%vQy4)ExUuT z20JJ|PGeMr;J=qS7Py%Y4@pSg@8H&)oZ8Gh^OV=5C(JYIm3r7-AP4=OlS}~_c+DTK z+T2iDGCTuBigBxe5PMrgRXkesz2+qH*L>JHVjfvN+ABdBjp)(lUvGHm&|2UBnb){7VRO|!x zqR^{@b&UwxwK6vOAp{$H(~LFrTbJyEVZ2i>-x@Q7fA*f3*a2)Gwry@=W<=;<6WnG`pJDc{0wdK!Nd`Zu1J$X#mt~U6-tXqP z&$7@DOJeTeX!?W#%kwJ*>1A>VOCY zBmn--8Y*haUUf`h)bI_)aDS_j!rJ_ZN?)GR?&H#cf;^=tHs}QR7cS`8+GIyZQ0`?f zS$;`lUIv3{O>d4qvmA7T`Sp(%3QmQ2HE~g9C7h znlM0o7Nh)Svcka0?do?{ACy3nRqo;?boIar!ucwonbTZ*N~JE|@vem_%oCl^7C(F` z99Pmmb&ufF%6t)xqMc+{u5K?vE7AL(vFn| zW#C_a-V$$X8}ZfsdyUp<5ZT|sh->rkT9MSwb+ZviW)3=dbn?Bpx1H>qy@1O2(X_WR zL3kX;WeNT_7nl9vBewJxvCf1_=FHWBuaNjFKSZ(kN{JS7;{nJMt!Rk;$i8UzOXD(k ztv{I!!G=kOe^TrkOJ`%Lvw)PwU-p1MACDIA#o{mCJLdh<3p+C`(iJwPE>J1As=lJY z5dZNoP}oWRon5Vg#D46JI97!;(Y({Rwf_Vcw0HiVpZ+e_t&FQ(dE754I+XK0M-z59 zzg4x8{=^+?1=w`yaGwY>%Rp%Q+r=+huqMY|CoHIch^K)nzz2%^3p{-L$1uq8^Up2# zg+5I7?5xd|>$k~^Zs+U8=1})LIr-HNasNDCGzI+UYW;7?I&@&_+;-4d@V#`iLE%8# zwO&_SY5nr6HHpkC*-$_8bONhC`i^-jgcrT_yE03kB-HN+u26Rwt{*Tdqat<6i(q?w zc{aB}bgL)5;AJi;331oc?Cu*uelgSGw~r*Bb@8!q5B$lU|Hfa@>2Tv$Tv?oX<~QCC zVUyBf;mSltjQgP;wJ~fQ2FH~Xc%1plkudH~G=MB`Bo0l{JjQO<7+w$HV9iqsr&bRs ztqqEGFzxWmcn_X6T-G8xA@*X{i!UbDP~IS5lsaj z>NY#O&>h4nH*Y3*vKj)-ms5?!<1`rcU}<-hY@}v-@AAxDUOAzaMq#kE{cjV#1D$4S|F z=+3^!Ov+zeW-dQQNq4gKQd`05h1@*-xp!>_5DxU;kIQ?iXR842^hKp!5T%2IrbtGO z7jZ6o)B>m`^_?gO=_R*|rp16K($3AT<&l^Fo;GET!JWakFE>fAGO<2y=D>L-7uD`P z&F+&6PMAB;2hW_L9ag@8WtiSqUf=_GxRT*D;jtVufkA=bH;?>rVAfC%FOFJ0~o`V#0Q`sb)A{dmUhRsnmb=ev=#l358c-kt}zdo}%G@c`=zSycE*|6Vi zQr)|HQ2|%;Kfe2IztJHqZxeGgEVh74o^z4D6P-aPIkBGJnI~WpaweZLGaNz?hEUZ_=z8 zbDCJ!AyIeym*b0S+Y=vFR^3!nfqk4Wt-eqEUxkQ1DuCu}<3Tu*#ATay<%9lleD^Bk zi2@QUe}Q-d_<4?{+thAcb1WpUs?A$g<_Q?&Z7#sA+pj3N$QwXGZp72rFX0?nXo+Rc zhi_f>cz%CA7zv8!Vno&e;%^&jA;vIs&fj|r=#7pYD|#OQM*O_lTL$~z3(i_hQeubc z7+q~e^9kc2>!=(NSi}}Cr}fo%DQMexYan_+lr#L&3mvnIAM-xNoz!aI#NWIM;=dH_ zm}5gf_|&1pa5XQrK}jCjFPknADJF2em#bVArjVKeip|zw9y~9>9c>=$@4QW`9Bsjp zUZ>`9VUY)O_F|4{g(bCzFLR+qXLPCyx6lO$ezv|BA`UoWnVTPL3=F#{Ku^>t>poey zn@?G~m&kwLB&V2-1xP^fWIEAGupZoQ&VY6?b0!+O3~THc{Nz*;2Xxg4eLB5=Kkn_0 z%RjA&hJGCKM~f3Q*_2OwMhf4mRG`4;3_6w!FPl~EThG}3w)g3a-y;irtutW-SvTBn z4STv0X(yZXJRS`sNFS7WnYYQx=R%96zha`6|m z8?-L=7C|+-*PCgbpnjl!UwO@>kJNWM(rZUbi8<*s0m(nYIG>y?jF~g~+{D4XX;M!R zkmxdex1Ik>7Sk}kX=X~%lDT&e$Cqm-6}fX6xp#cXrBhY?)83za1~4&P@6eBJ!phL2 zPChtJBz{j&_9wDqtg1J2VXoWhg*DrJNIzYW_}P2&w7)0lVL)6JE8!hgUO4IN$d z-0#}uj%=?99{+e&y_D$K)Y3AV=D#Tr=|A^nu5bV>VhT|1PAf6S6N>51_u}eSBN|JV z0eUSX2MtUVaxY-5#CBBsYnjoYfU)crth9N*rjuC8DCx)I-m7xCR$&7IPr;|CiRPq# zw_%*JtZOnp*?k;nSgXVew3)=@=jH{O_Jc4>id5PcUJrqds)l`3>Gw!rBIriQaLJwQ zcPbG6$-H4kLEfRHl>9XuN#Auon2UMhF|@0Fvj|Jk!6qOfqa#MBDOFK}noG$Ko~PqK z(~q@o4TYY4wIM_W8ESytJ;`4B6MHMGj<^_niQB>Uq?J=b{r8d4pxhH}Rf5)B(_qR@ z3L^-B@!p|u$c4wn8@IAkMj`fWkz~JFS&|_qCUTfv9^?(jgEhp;3%fn4gc-DTeLAFS zqH$~1Yf?(%@Fm{9&}JXZBmP?)3r}%>?^kk7uT{)9FWna07np9pFQl!jJvaU{iWF0< z?TP-J-Mlnnupm3tK}OZlTnnaNNBisrLCK9AwkMRrY>2vKc&V36^7us9l2cjjZ=~nv zqf}d~j&N|tNS@l`ajeLCy%XU;z!!OYpXWC+63&Jx|*1OWQ)-V`Kz zd9QJIpbG&NIU$dV4~w?HI_|&S~2ZacO zSCR#6{9R7M*^AFP0|IiVRl+si(mCXsy3R%10I#6uN`uBTg(bdy*H4x-7{Rx{{6;8Tg;OWt5=Ez@K-AnnD@t&6ANIKTus*m zr)mLVX$K*Fmz91alcr*4)YP%f*M%@^vrFkSFW?5|sgpckM;rXBB^heSxm+SN|EPYz zwBk2KEe_xZbO+Y>s4car2{~MHU^EJkv7c^piCES@TX?Mbc=dr05-5{KXU!zSD=xjE zLzj7oc$H=@GOP9_6{T9i!xa+EbSui`#rY-dUani<6raEz7~#x_oxzU7|Eb8-P$_0avbQT^SxA5nVB!UjYQvSI`+FPP=;9$BX9pBJa##A5`pB3pvI- zq%v}sbVxjjXhRf`sZ+|wniRx2JB2w^G=0+)!u`8M!ne<7f-=;x-1vm>9eha!x$ree z4s&mgHXHQZ-s2V5b{1Z(q^7-*P$%m zeW56PJT65^A#mq6KV62IzwyP=6QftzIt5%gJiLIZ>i@vogYK~l*=|@uIU`;k!1{C^ z>ilExo#cCup}j}$h{3xry0SWr^X|Oa5^og+N#DQ4>KNBQb_?}-<*|fw+MtZWnB(Z1 zG*MQ6g03dSxJt6-RF|0X2J&ebExVdE-HAVg<0-QSP!V8GLIHUpa7mwbCwfFuol2a# zWso{^{erWFY>it$lAC7#I1Wkf=kX!TrcQ@mvAG_g5mzLye}r!?Ui1GZD0LL1ugiNC zp(OtB2J`r3=eMo#f>I4+f86l6I@K>9me9fxuN-~5^`D3X2smy|#>@lhFd24!!?2)u z=Vn}Lt+b|{K;o7@&S!YL_phQXhU|IIkQa~{rB%A3y1dO!Y=K7hP^Z*@rHoszeQLWs z;9P+e_C@DvjOpRbUqk**=epSt2a-PiV9iPu$+wKf&tKS5g;k zUKkBUwL{oX&3rd~0O`Ss6CK~3MgqbfS=`a3NklB0k+4T+ts2rvg}eD(&}U4@34OLn z)vGoc;~VYE4FBg+!buIZjyN}Y4WA!`o<{n-Lw5On(V8bF^rgTkErNYWLnM!1c=e)s^$UF)lT%z>IAgDj?E(i52gDlqcS5d+TZbx;<-HTi z+h6EiY2(EEF!iJH~S&{#dil&i2{3{E;z-?VbR#qJY)uxqL2(omCAGCS3&jX#5?n$TI=q^R#m)L}tN_s#x3sLQl zf1)C4Qm=yTt>gJU*C_1MvKlk-PdfdiL>g-iveua##T1G0P!u-)zK$>C>CSO$wCscl zuo}p@ue0b%$f%npMb&XPr4!hHS)~V9oj;y}fZO|sY(ulTWetDfx<69itlyAT(lXv^ z4oJy+MU&1xVH}rfOn1$i}Km!GUBKVI6PzJV-QLjqnGV9Yexd@JYj{9)oVS>iRAk*iFa~ zAcx9~(g#_)H&of9)H8oF(R&3{L5d*Nfm28~Pp9p8kx^$qf;U^PyItEwW1FZ_!7^bgD3WceRgtes}P_9gA}dfBoX+;WyJo9`E>Jm zBvKmH!`^C()!;eV$03fE*YS;qMmtAiuZ|QP2X3qtq`{xmyC@qC*9A{+z?{1JgIRZf z#LvyBCe}>G8VR(%PlM!VcoAh%#i4L}HGswJIF|-o(GQ}}_v{fEETI)+t0H0?6jM9k zl#=1A0oal{i@;FY`X$5MCG7TuL8k`3_kl~Y(COh&5V!w|h4*iNj4T1NJ~^&-*trsQ z!S1)Ap<|P9528`0pj~>nG(gj&)47X74s&h$)1j^hl7d<`AEU(pHX=R>1< zA~9f5#Z7f&Q!$nYo*9phAp$Lpg&MnH+qtXklfh&kIW_Q&xCo4#?`3}p=~$UW6C~8E zH}2aI9`AgQqmu?pT?ROgdaU*x%f!R zC2yOZ7m2#Rl3oCn?cz5Na$#I87j>=CUox}a7V|sz0wt%28b-WgU?l_h;{>Zkrl+bK z-cQsV4fhM~J@awe3=#>xorDVTT;I9)wNzx?*D0{KoH;?mGN4$A&`*P^Jk#QFYvcVa zme+WlSuS1*N-BGk2TCN@g9|_j( zv!JuZf-s)At&Lj=-Cj7{{2t-)yY9k#-V0NPNtBq1&7=U8z=8XsLEM(-Gb`H>ZkK1F zV_@LrH2JOuFb}cp;>d6WpYEko%xU{Q;X@6!tQUJ^tV13uW+guCSPr=-eeycYd-&^f zuq8WQxPX!cG0=QBf@di!^=q@h@9|?^T+J-V*#41U-h%z)k?=dl261=?;nL3feN8R@ zpHdv;LX#_zASSxIdbJp-MJv{DDD>Px31vHCND0I#{F)Q-mR@{RZk&C3o=YtAO)Xjg zD&-t%>rbiN$0;YpXgshh#Z&`lHqzGHycyWO_s0 zK2F=PtJPNU0lz@`BjEy1BwIzroW*Bv&bLxB=(GGTc}KoGlcIb}KXlGK50$_bK%$Qg zhj`4gB2s;S=rYEqn*c+7g4Jg`xR__VBm08ScnL-;C8yJEuT7wl?J2`hfLFw|-naq= zU35{U2PZ^WhMHZ7hCvqzpupCFZI|mb)d- z7ep_si(VG1ueQ5bYacWBJTuRW`ww{LT<^|Y@6MdL<~rZ=`PzxUbX5I&=vc!safQzd zN(>+GB8tH%R_D-}LH824u_0&r=SiNbC|LX?@6mTi$n76ydD zRnLtvpF_i8tY$4;WYnE+6aR}_ip zjQ+JTom#2@j&Lnh1XDE^ob~2)6{TjF{!?wi)V_&z9N^mTgLhNQx=oFHXdkUBe#pTY zT+&CS{dgAnDiFerXT#s3K%yY(Cq&oFvtG3KuI)h&JcH34XGEMH$KJ`HpG{!P4mnJe zNA{lE!F0cYwgJmOsVbQ03bIpgxPV?!%j1&xr}t{ZAlUM^4Y`E(^Yr-?I0bf}Z!Cox#6wh$C^*?P@M?RN|-rFTEH)^3@%#3&t$keH2IsprF{VA1!be0IU9P@J<^y(~C+CP!=FDw7A zjk&;$2!}d7 zsD2r=wvnBid?a&tQ_kdMCD-}WUmNbTX!S-a;3sIWV0EtXqU<}UIU8ZmXAmelq zYH(3RRJD#UF!idl3=5rrb?yN?o9@r@YveMzeQ>jhF&i zm*_<}9p^<8q1bpDaWG1X!)HVy>`wxvVU8n%QpZWeYhoLgRVWZ_oh8RG)WqWM9($NS z=`}GXE+|f;z%@`Rr_$Dv$N6Rk#YTa|y}(*q|I7Rm3VEW|@cwnsrkdw&h=w@c_8rsb zwr;aefQTEhe&>BWKAIu4$sgQeDSsoe3HP0kl{PwBJY1?L$@K4@0-8e-C$n*55{?+C zzeiL5Zc3Vxio!i4E_#o!UG(?F9JBiFzGwf-VfOD-)oJGxi!n{+-qH=0(k3f-TS%$k zuZ%E{N_JSVPM3^G_E0Rs&c`a*ZiWZEwiQvrc=Oz06?_5DZ{bIi1i!pX5ui8(hrWgH zBU($H%nn@I#62zYh7bx$yjA;mr{hESDiqQN*yaskJ0gkY5q(Kma9LgXWc}mNl5bau zW5&}6f8iH6301}qP3z51-9D2t8#bdJR;eD#YoV=DGWzg`PN*(9C+_9udZ;~IhrlR) z(Z@y?ZLuusdUy~f3Nt_*j0`JDOq8Qv=?>YnLc5S`!ZG#Zn?ipteCB4ihm})QC+1b> z{F7f#cLfJBhJY_WYeBO|uvuldLq9$eF_RXLsclkB_Wkf(W1|JK}_ z*T!9YG@(+9z+>9oHTbr5^Nj_K!y6rb*@NPyUt;zFGK?1KO6(2T?+0(9@)%-2Jg~mj z6E>%tWR@qU^{54&DhHFd=#mC0>(AV)HXd|oJJfQH`Bw`e#++{Mo-J&Nkzy5i1)m=D z!3c9NL#nQaDQ>uYZ+v_DTrOMv3rtqgr7mGVvOsZ5JCwFmbB2JfaU1v#6!p0^2`yFQ`}!81g4$KfCr;s1H-J=* za^q;yEFOU3y zBp->@_RdB<>pUNg%Z9Q^{J9uN3Grz7#3eS5=u3|_ruO)WF+i_qHJy8L3w0Q=#`WlO zubyLv*L~%m87}q~lb5GED-!PC@w;6^E-beaIOT1)g#3~0Mb=dp0j!6M>iQjrUBIkS zcBQ?BOA6y*wpOKtn+w1PCem{ozeh2-PgElL^hc7$D53+r$ydpCB2U@WL`6N)qK(kz zpTg=zm9MaI4a0LboakoU*D@Zbj$HZaYTBdSJeKpUgxij%YCVceHs4;WYo zI*AyaOnG9gu;L!+N&_B)^}BFobt?!c5p)u3b5=#O>TA9!Flu@QuURzI)&kM}Q zPlOR!Y_jPJ9LPe~!CGgy2W6B%S5bUz&@E@mGCJ0N%;ogx5^y`6L9(5soj>Ece{q$x zN@Xhif*9TVSKeblG_s}bSO@XKj#Il!4sK~CTeRkx98w)hiV9?2k!xemW%5?MV_dBX zBuh#4E^A!7M2fAMpO3*JsaMEHFq^D;;PdQl0--R()9&$Ce4w8un{v>}QA(8xR@a6V z#{azX>qFe4UVJ8R_*jtw zAS^~lT(|WE9Zz;69R#LT>4k zt^xnoW5P(XDWWHdgx9S6*0CMmz_=pD_Q0Cty2w8|q?NBj9b(p>w8geaFvP$ksVwAW z$+Gz5N~0$9zgA7N(Ij*i>b1LmTb8ncQqkfeGpvm>{;LSW>O@~Jd}&b2H@LhdCb7Di zML~g~B+`Bhf38p5V@Q%amW^N8nD!F+`wN(;jD!)z^Zn(@0M>^& z49UBa#v&qhQ3TIs7!EDFl_z%PYTE0k5IsLMtS4=<2j~~tj*Lu>T&Y?a2qz?8zD^!r zI8_mWJ$*%JWl1k}{UX`L+FC-N9gkLe(Fh=ouMQr!s55qz)nl^GrM|3TpJopUadCP` z{rG0q2Q|P%T|wDiZll*NwcfKh6m%o9jmdamqG`Usgdj+1D_V`(2PlDMxVm+)_T;Bt_vCs6 zFT04>o2PAOgtQqgduds6U$Y0bSiX~gGvwOY>Nc!8x@UYEE=}*{rp^`FJw{%StL4bbcDF#_f7AUO3? z;e5`SGn(LO&Dbh-j|bU>=jaIjZ56zLho)4`oSz{!G*{a$EUY3n?cAE?hXFC6Rhbu| z3%Vpmj2Ibgfpn2opgVPe<9}Wv<2JhsH3Uf>-|F&3a;eX&bdp6ih{_SPLLW*%q4<-# z6%LjsDj<$O*b$DM**XSluwR$^;wop}oYZkeYxF>KE+lQxoPTc$?YtJm5u3jrH`Zbc z{Cr%*DPY0G)hWL#u17(|@1|Hp>cyCKdydu$OPK?xuCEg(ZXcKdks$na(kL168WNS= z-!ATrZ@2&Wtur^@6ees;wtYBXw#U~vPNS9&eE+OvI2#Hm!mxh$NRwzp9*TEf|E#&x z+W1QfYh`e!7FsE->0>R=kqU;NC*a1f{M$Vl{luxQ{dI~qXj{#6>}83i@~xNy4Q_O= zvn@}^;=3+?gG&51wvkMRkXjb)`hI(oH?8{#Iy7X7+P)(Aa~I@>O@mTu#(nzUZR%Ib z1Ac6(5&+oc1^ewmT?%S!R^9^GYjw}_K5C$N!EbeD)8@x*{}*R*QPPXokz3M1biY~n z-XPn1RD@k6x8`wE94kSEVfT;nvxlc5T5Vm66BKi`2e;k`p@p7hV^yCDj@9)786GJlEsIFB)s7Z#_8CTMewTT4*PO3M$E~t7{(L zo~d@M99_g7=IYDjx0GyHy%&#ewkgX+EbxUHK|sA@z> z!y>*iJBe9s_T>SdXfmU``N}0Vo9GmG$T#&_)*j4C824RQ=uq{KASXFid0o^>%kP4d zAhsgHx5wW@N=Twb%>#>Hk*VJ5={`4hincL$Zfk>y9C&n#ORy;zkK+mM5iL{5evDr*!ITJm>4{{8(pSstp2eB( zKA-X?etv#aCe5UYEM;Qo1e3TMu0(;^P4Tyk@=2A6S~sj_PHcm#9}7q!J4H=?hWA0` zl8BDiRi?hbojDXZTAvYwk9X^wWs!XB-}LqiDYHoraqGEKe`6CSEA=#_UW21Alq+UF z59Rru4vba#Iz4kS5*@zzpZON>@-1152a|c=ZOX$4z?!c`KbD%>C#HQuPmVKD*jlQL zCH^2F)z@hh-*g5BLs-T+8-saC<4q1NO z@u|nhT;SL^rLz4^}yI?a6t4jG$e>JOp~ z7aD9bk3x-%MUaW-k7$*N(H)FM5Q+az%bK+3qZYvN7=4R0_?T~GCSc5{nGk8g zRxF&xN8bNcEs`*fjGU2yQTVB_nUUGAvS07@ZF>rO|8K&WQMlqhVg9SE?!CT4Pr=An zwRFO`r{s+PlbF4~@%z2LT~ERPj=^&5DVX@GmP;7-zhkgqiXS8H3~cwMpz{OKrJUO9 zp%(oEO76?hQMnnN(ki zJXpznEYS68!ls~(CO>-JRoZie3!zDh+&68gRhE3&XRC)}Q?WY~JC9x`!EK(CB=irK zMTe(8o8gHIX za}#)}aX+Ee{#ur+2g$7ndBdF&^NR8NlC%D^*nisl?g z^BAoY+U;^OmLmfn16)AFeRfxf`>BxOj}e{C2hywL@AyuOFY|7YYHf!TNJ=5jiNIb| z+}#!{OFA%&J;Os=16N%&!w7yM`?Pi_{2d2qO>xG9G1&M$tmV>6L@;T}F7)+5()BrT zAme~(Iy+;q$ezDp#TL)DUk*KSXh$f%iDw_WG=9?S%V35~E-y~f1 zm)k+z(f3N%_po$CP_3>I^w7LDh3V=8@JE4n82GGsyYBX6q+GlEvp-{oPiQ84sfOsd zzZr)x)iU$bxw^rBi64T#U09$J>Jb0B4E$4n)fTbwU!|OzvtZSHp&P}mzk)3%5;w3> z8eb!7)az7HkJj{PWYJdLiH?k<<&uT`G%I#uGnF0QH@uEF7jhcvC4^1A6juB59oz;5&f8Cwn0 z)h~(bju80!Xi*uq6wyo3O?y#zgKmd_tbu4OV&u!q%7CR;XA&J~sz}9S>!?G*yH|1= zw~g+lJk2L5$QX!H^b)x@4M(BjJZ1CXHCr`^yft4=XV(xI|C7$b?`K#yP*Er`6^Ii| zb2;*PxeUXKEJ0w{nm^J1_zi(67#}O$;XZiGv-)VN&KZvv!g5E(waZ^AR9G1H76fNT z!b{`Xo5f#1k$G1z@3Eu}o>Pa*GX_W^3in*O9+GBNefRr%bhr5TMUIA`)wwqzB(QYaugvnEN{&#oh|ZEGx@5oN&}>kVc$G9=?i)DUW%*sH;(#d^hma8{cmW*g81D? zSIi2PzIQp6qEWOLL)pPGvqbE?H@MR;1Cgv22XVT_%K3OUVkN;M!)@$A6T|7H-~DFR zAV&_v`{ptwgTZP3R-%!KaIc^z#La!jZ@SDx@|pTGiX5+t7!i@6!2VY=V0~0tUH9vo z1PwW-qF^29QT}ff_!kC`ZF8OWxBFuf%w3N2nR45fbtib>|iB}^us8Qv;boD0aP$Z0ALL|1QD^!S0kNP*qYRfXrjy2( z^Dj#CL{lqwm>dPIGXxbAe&psK(b}p9%8i8n^5a@q&Gas9<)ll&S7xPj-!-hjoKmj%c+oFDA8KtIS2L5 z*ACD{XxPl}i(dn+HhEI0JHiZ?YkYxjv9@NK%O0S>C7w+}(2H1(VfE3Z_}lIdYoXd5 zxbbcCS3oaS!VbxuL+`rw&4M7Xis_RsN^1JrFttV`$QCJ#4kL6BAul%&9a2SRCT_=# zdbLL+U#TV>ck%2vpBWk$=vv)3`PEW`re<)fQ6(hGe_`Lnj<-OPL@g!=?X9#xhc-!AdvT5NqoJM^pSm#S>FpFO{vb)6;nWGz(`e0gMNhy z@7Hxk#uZ-$$3jD{R$`&)v#pOH{4}6G0N2%)kElX3?fK%h>4t93p%<8icnq4n5QE8| z@YrM5eyo^q#a4|G*}8oVW-SL+VuU&bmEWt^@x^i`1E(c^g zV^V!9pu2#N?Ynm%nf!$6?F#s~suqN=58qXPvR1C%NOx85pBDVrbWbWZ8mV&Gk+TFZ z2;oLL9?3JVLk+y(E%IGB?yc!aV8X?Jyt-8>pgd zyXL$fjgpGGF7vX^-}-nO7%dB-@e29rrFc)tdA+u<6`Bl}ue!I5Dj3MUun`-@Swnaz z-Uqx*ny4PT<#gQ}HfP5i2aX(gRm}P*2#B6KI)6y88My-2@SK)zZ{4DN`Xwo3$#8}L zRSw?h0aGu9jkig#_BXH=y2tQ*nnCWrop~q~T-g;1qiOighBU)d_Q#ytsYe9KIbkBey>O$;a zW5jolUv3ROXWQv!K*0pX)^sqRh@H+04W1@qUc%t<)%?ibJ6B)!=kx0|NQ}NRRrPwi z!3`oPP1^1qm3SfciUlrk4p5F)9L^Xg{;?IpOapaHMt|z=_}Dst07#%yXBCX?6wAzU zZCh=-U%p4Eb z@IuLapcIoPM}F=Ph4`f3H!rT~<@1RbqGK%wz4oQqbIag?!5k_(APQ-J_DMRC@?t6`n{Mk^VlN^TG*uK96oKH=SbEWk_+qANrTy8f zvdqp>{k+2fx}Tm_l;p!qEkT=k)SJXZcr7 ztA{jRf%ft{3p3}$A8_1l-A|CkTep)T2i=Eo*sVbKldb^nEp-ZDC*PU?z-6`8c?HHi z8a-=i_l|ZTteQ@^+i59dW+C2eR8>4w^9@s*zz9#a$Hl%6x* zS(dqBrJRQmkeR%7fJ}G%A^(Ia!u0R5-AWUe{Xi{dg6jpjNIwEe_PK#3H2BaJ z3(7s)I(GLnKkuBOAzn*NFJ1QvVEw-TaNjx`i)EVoYUy$B?zWw6o4i&6Kkeb~3TJgm ztv)1(>~rHF#-TU>&-aF}1aFFRL3MJ%Y%T1<5HFwEgxAY$e|G#bR0^sFy%by@*OZDh zI89Z`Y_GW+1rIk$S!K?TH5k5kn1CeMxa+5`c#r#8;SRMwu7u!*4Ha}ub)JmJOs@x#%CNzk;meO$hwG}BklKnH7Xe*!Ifbr~tJUGpkD5!A{nl(sm?{Uo znH443<2hJMbJ{a_tfa0k;KP@;!7r#N%B5Ye+=tF=fLk_%Y#X18I%zeRFIYRw?7jZ&T_TzdI+EEuOT9|}|0b}+Tc%y*dR^tG zwID;=teT@Zam=)~3AWi3d-gIIZNk*3TfWZVtv&YcBU~!%$UQe$q^XTS1}k?|NgqwfG1LR5*4%JNJHved^8J>gzr%2cK*r&W+ceO?VH~)giO32-EhaiG~lO&o@14NrZ}YWAmW6F7qcy$)xy7?bOp?ras|lBXA^-=!*F z*+jNG&u`3D5(38lUAc`#`gwMb5*v3}Gm|&JV`;KBUL!lH9jf=n^fSh4m4;w274|J~ zn2+q)q|SsC9H@iIMiMFS%69Y3u|)pC>Qjx*F!UAZgQ(cdtHyy-;Iut(u zEh&_Q(UrkLFW%W8EIj_@=395z1at@bLz1=eX|0d9o7@KPrh5tF`NRA*1o8cFR z*R3yXMg-FtVSdUF?7~;E_cAwL27LKdzp$`2L08*JczDOZq8Z15M3Ws~@QCG|rirLH z9o}oZn``?HfzXj74HPoz0 za|du7)qv?uM_-Vqh_a2oyPeU{sG05Q&yqVFp6f3W zcG~V%n|N{Xa9AH-!6o?pe@DOo zjU9CeAR@w}`zUBthINY$v})GZPgpPw6xEVnWo*($vVq#IIbt4pgPacA88~mfmxVaS zG0gr{azIt4RFbocXS+HL=A?4j&kk5*6bU|rw+nZ*O0fC=wP;uj>5}|wFY_laFPHk= zK&8{-m{CkyP|bVLy4tVm^^K2UQvlaiS*^!N8lh+giBHq3%-_5pw7%f}l7b6YDO#vl zc3g(C4t(N$8$>i54;k`U5!WeMZXhn!$z#jrE;1!-(inYQa0 zRsLAQ=-o^|t2sc`dajq5BAaT>=GbP5Uhhm91gIZ}8Pq~Qy)8P(8Cgfh`dPWXTgv_w z32=8r3ifyED7zZuidrOnCiyLyv{gv1{wMJWmJV%+?Db!cf_1u}zE>KC9uR*?PYk zn%+AoQRZjV7$*g>kmf^#6%_J9MO*`y+nP<2rH(Frm}jT40)J0VkFT$1bh(;VZbTH) zy*W!(%a{jq6~p#mimhYok|f+zr6NaDdfjMUw;9+Z(2-gDj0btZfr7O`@{GwsShf5O5F)X;><3CbWxK>@5n>4se!s743dn0=+DzPG{Z$ z)~P7Ga@F>glL!!d)!sQsk+ ztDAdbudhSwwa{TWU@jTC=O-rBbmcnu&iV5a=5~{Z*7lgu(vGi#x$NF`e%fL)7PXZ#T#o$#w*X}Ny z1?Y|uOyd4{k6hncCi*UvL|4jkn`V<_)A3Lmoh{GQtL|zseS$frmKB$FOjnn2CMsfn ze!ZOaAl=$k=V#&>TW#>m1&|*{G5Pta3E9(w*Mx0J=-F$r8qtfw!GpICGam6q6(7Vc zdbw$O4*Y^vCD*wUJiSDq&xmW-U8zU9hn}jI{FFn!tKttNB1lA{>-~SZIl6ePcDQ9X zB}%#-xg8x@a4vcKHyyhz)9RIN=Vwz(h$&G?v!#AI*{V7HC=mwpd!THmoZ% literal 0 HcmV?d00001 From 323ae3b23e7bf8a868b3f56e8097d66c581a890c Mon Sep 17 00:00:00 2001 From: Gene233 Date: Tue, 21 Apr 2026 17:09:18 +1000 Subject: [PATCH 2/6] Phase B: compact G x K IDF/IAE + sparse-preserving TF composition Core change: labelled `idf_prob`, `idf_rf`, `iae_prob` and `iae_rf` now return a compact G x K matrix (columns = unique labels) instead of expanding to G x N via `[, label]`. `cal_score_init()` composes the final score through per-group column blocks in the new `combine_tf_idf_iae()` helper, so at most one G x N matrix is alive at any time regardless of how many factors are compact. Additional Phase B edits: * `tf()` rewritten to stay sparse on `dgCMatrix` inputs by right- multiplying with `Matrix::Diagonal(1 / (colSums + 0.01))` and relying on `log1p(0) == 0`. The legacy `sweep()` path would fall back to `as.array()` and densify. * New internal helper `rowwise_notin_max()` replaces the `apply(mean_row_in, 1, function(x) max(x[names(x) != type]))` pattern in the `multi = TRUE` branches of `idf_prob / idf_rf / iae_prob / iae_rf`. The old form was O(G * K^2); the new form is O(G * K) using a top-1 plus masked-top-2 trick. * New internal helper `pmax0_offset()` replaces `expr - thres` plus `expr[expr < 0] <- 0` in every IAE variant. For `dgCMatrix` it mutates the `@x` slot in place and calls `Matrix::drop0()`; for `thres == 0` (the default) it short-circuits to a no-op. * `idf_hdb()` and `iae_hdb()` now own their internal HDBSCAN cluster labels and expand the compact `idf_prob` / `iae_prob` result back to G x N inside the helper. Their external contract is unchanged, so callers passing `idf = "hdb"` or `iae = "hdb"` behave identically. Their internal naive TF is also routed through the new sparse-safe `tf()` helper. * `combine_tf_idf_iae()` detects each factor's shape (scalar / G-vector / G x K compact / G x N full) and slices appropriately. A small local `%||%` helper avoids adding an `rlang` Imports. Verification: * All 85 existing structural tests remain green. * The 17 numerical-equivalence assertions pass against the legacy snapshot; the `return.intermediate = TRUE` check was updated to compare the expanded form `md$idf[, label]` against the legacy G x N matrix, reflecting the new compact internal contract. * Sparse (`dgCMatrix`) vs dense input produce identical scores; the sparse path now returns a `dgCMatrix` throughout (max abs diff 1.4e-17, i.e. machine epsilon). Breaking internal contract (non-exported helpers only): `idf_prob`, `idf_rf`, `iae_prob`, `iae_rf` now return G x K. None of these are exported (NAMESPACE exports only `idf_iae_methods()`), so the user- facing API is unchanged. When `return.intermediate = TRUE`, `metadata(se)$idf` and `metadata(se)$iae` also become G x K, which is a user-visible shape change; the NEWS entry (Phase D) will document the expansion idiom for callers that relied on the old per-cell form. --- R/score.R | 67 +++++++- R/tf_idf_iae_wrappers.R | 179 +++++++++++++------- tests/testthat/test-numerical-equivalence.R | 9 +- 3 files changed, 190 insertions(+), 65 deletions(-) diff --git a/R/score.R b/R/score.R index d0df49b..e777f04 100644 --- a/R/score.R +++ b/R/score.R @@ -106,11 +106,74 @@ cal_score_init <- function(expr, tf = c("logtf", "tf"), iae <- do.call(iae, c(list(expr = expr), par.iae)) } - ## combined score - score <- tf * idf * iae + ## combined score via per-group column-block broadcast; avoids + ## materialising the full G x N copy that a naive `tf * idf * iae` + ## would trigger when either factor is a G x K compact matrix. + score <- combine_tf_idf_iae(tf, idf, iae, + label_idf = par.idf$label, + label_iae = par.iae$label) if (isTRUE(return.intermediate)) { return(list(score = score, tf = tf, idf = idf, iae = iae)) } list(score = score) } + +## Per-group column-block composition of score = tf * idf * iae. +## +## Each factor can be one of: +## * scalar 1 (the "null" path), +## * a G-vector (cell-independent, e.g. idf/iae/idf_sd/iae_sd/idf_igm/iae_igm), +## * a G x K compact matrix with colnames = unique labels (idf_prob, idf_rf, +## iae_prob, iae_rf after Phase B), +## * a full G x N matrix (idf_m, iae_m, idf_hdb, iae_hdb — the latter two +## expand internally to preserve their legacy contract). +## +## When at least one factor is compact we loop over groups and broadcast +## only the active slice into the corresponding columns of `score`. When +## both factors are cell-independent or full-cell we fall back to the +## direct algebraic form. +combine_tf_idf_iae <- function(tf_mat, idf_obj, iae_obj, + label_idf = NULL, label_iae = NULL) { + N <- ncol(tf_mat) + is_compact <- function(x) { + if (is.null(dim(x))) return(FALSE) + !is.null(colnames(x)) && ncol(x) < N + } + idf_is_gk <- is_compact(idf_obj) + iae_is_gk <- is_compact(iae_obj) + + if (!idf_is_gk && !iae_is_gk) { + ## no compact factor; direct algebra preserves sparsity for scalar / + ## G-vector factors and only densifies when a factor is already G x N. + return(tf_mat * idf_obj * iae_obj) + } + + label <- label_idf %||% label_iae + stopifnot( + "par.idf$label or par.iae$label is required when idf/iae return compact G x K matrices" = + !is.null(label), + "length(label) must equal ncol(expr)" = length(label) == N + ) + label_ch <- as.character(label) + + ## Slice a factor for a given (column subset, group name) pair. + slice_factor <- function(x, cols, group_name) { + if (is.null(dim(x))) return(x) # scalar or G-vector + if (ncol(x) == N) return(x[, cols, drop = FALSE]) # G x N full + x[, group_name] # G x K compact + } + + score <- tf_mat + for (g in unique(label_ch)) { + cols <- which(label_ch == g) + if (!length(cols)) next + idf_g <- slice_factor(idf_obj, cols, g) + iae_g <- slice_factor(iae_obj, cols, g) + score[, cols] <- score[, cols, drop = FALSE] * idf_g * iae_g + } + score +} + +## Null-coalescing helper (kept local to avoid a new Imports). +`%||%` <- function(x, y) if (is.null(x)) y else x diff --git a/R/tf_idf_iae_wrappers.R b/R/tf_idf_iae_wrappers.R index 7cea362..092dfd5 100644 --- a/R/tf_idf_iae_wrappers.R +++ b/R/tf_idf_iae_wrappers.R @@ -1,3 +1,52 @@ +################################################# +#-----------------Helpers----------------------# +################################################# + +## For each row i and each column k of `mat`, return the maximum of +## `mat[i, -k]`. Vectorised O(G * K) replacement for the nested +## `apply(mat, 1, function(x) max(x[names(x) != type]))` pattern used in +## the multi-class branches of the labelled IDF/IAE helpers. +rowwise_notin_max <- function(mat) { + stopifnot(!is.null(dim(mat))) + G <- nrow(mat) + K <- ncol(mat) + if (K == 1L) { + ## no "other" column exists; return -Inf to mark degenerate input + out <- matrix(-Inf, nrow = G, ncol = 1L, dimnames = dimnames(mat)) + return(out) + } + row_max <- sparseMatrixStats::rowMaxs(mat, na.rm = TRUE) + argmax <- max.col(mat, ties.method = "first") + masked <- mat + masked[cbind(seq_len(G), argmax)] <- -Inf + row_max2 <- sparseMatrixStats::rowMaxs(masked, na.rm = TRUE) + ## broadcast row_max across all K columns, then overwrite with row_max2 + ## wherever column index == argmax (i.e. the excluded-self case). + notin <- matrix(row_max, nrow = G, ncol = K) + for (k in seq_len(K)) { + hit <- argmax == k + if (any(hit)) notin[hit, k] <- row_max2[hit] + } + dimnames(notin) <- dimnames(mat) + notin +} + +## Sparse-preserving equivalent of `pmax(x - thres, 0)`. For dgCMatrix +## we mutate the non-zero slot in place and drop structural zeros, +## avoiding the dense allocation that `x[x < 0] <- 0` would trigger. +## `thres == 0` short-circuits because scRNA-seq counts are already +## non-negative, which is the common default path. +pmax0_offset <- function(x, thres = 0) { + if (thres == 0) return(x) + if (methods::is(x, "sparseMatrix")) { + x@x <- pmax(x@x - thres, 0) + return(Matrix::drop0(x)) + } + out <- x - thres + out[out < 0] <- 0 + out +} + ################################################# #-----------------TF variants------------------# ################################################# @@ -14,18 +63,27 @@ #' @param expr a count matrix, features in row and cells in column #' @param log logical, if to do log-transformation #' -#' @return a matrix of term/gene frequency +#' @return a matrix of term/gene frequency. For a `dgCMatrix` input the +#' returned object preserves sparsity (`log1p(0) == 0`); dense input +#' returns a dense matrix. #' #' @examples #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10)) #' smartid:::tf(data) tf <- function(expr, log = FALSE) { - t.f <- sweep(expr, 2, colSums(expr, na.rm = TRUE) + 0.01, FUN = "/") - if (log) { - t.f <- log1p(t.f) + ## Column sums: use Matrix-aware accessor so dgCMatrix stays sparse. + cs <- Matrix::colSums(expr, na.rm = TRUE) + 0.01 + if (methods::is(expr, "sparseMatrix")) { + ## Right-multiplying by a diagonal with 1/cs scales each column + ## without materialising a dense copy. log1p preserves sparsity + ## because log1p(0) == 0. + out <- expr %*% Matrix::Diagonal(x = 1 / cs) + dimnames(out) <- dimnames(expr) + } else { + out <- sweep(expr, 2, cs, FUN = "/") } - - return(t.f) + if (log) out <- log1p(out) + out } ################################################# @@ -155,22 +213,24 @@ idf_hdb <- function(expr, features = NULL, multi = TRUE, thres = 0, minPts = 2, ...) { if (is.null(features)) features <- seq_len(nrow(expr)) - ## initially compute naive tf-idf - # tf <- (edgeR::cpm(expr)/1e6)[features, , drop = FALSE] - tf <- sweep(expr, 2, colSums(expr, na.rm = TRUE), "/")[features, , drop = FALSE] - tfidf <- tf * idf(expr, features = features, thres = thres) + ## initially compute naive tf-idf (sparse-preserving via tf() helper) + tfidf <- tf(expr)[features, , drop = FALSE] * + idf(expr, features = features, thres = thres) ## cluster obs based on given features - cluster <- dbscan::hdbscan(t(tfidf), minPts = minPts, ...)$cluster + cluster <- dbscan::hdbscan(Matrix::t(tfidf), minPts = minPts, ...)$cluster ## factor cluster cluster <- factor(cluster) - idf <- idf_prob( + ## `idf_prob()` now returns a compact G x K matrix; `idf_hdb()` owns the + ## cluster labels and is not reachable from `cal_score_init()`, so we + ## expand here to keep the G x N external contract. + idf_compact <- idf_prob( expr = expr, features = features, label = cluster, multi = multi, thres = thres ) - return(idf) + idf_compact[, as.character(cluster), drop = FALSE] } ## ------------------labeled--------------------## @@ -213,9 +273,8 @@ idf_rf <- function(expr, features = NULL, label, ) }, rep(1, nrow(df_n))) ## mean counts for each gene in the group if (multi == TRUE) { - mean_row_notin <- vapply(colnames(mean_row_in), function(type) { - apply(mean_row_in, 1, function(x) max(x[names(x) != type])) - }, rep(1, nrow(mean_row_in))) ## doc freq for each gene not in group for multi-class: max(mean(N in Gi)) + ## G x K: row-wise max over "other" columns; O(G * K) vectorisation + mean_row_notin <- rowwise_notin_max(mean_row_in) } else { mean_row_notin <- vapply(unique(label), function(type) { sparseMatrixStats::rowMeans2(df_n[, label != type, drop = FALSE], @@ -224,9 +283,11 @@ idf_rf <- function(expr, features = NULL, label, }, rep(1, nrow(df_n))) ## doc freq for each gene not in group for bi-class } - idf <- log1p((mean_row_in / (mean_row_notin + 1e-8))[, label, drop = FALSE]) ## IDF scores - - return(idf) + ## Return compact G x K matrix; `cal_score_init()` handles per-group + ## broadcast so we never materialise a G x N copy via `[, label]`. + idf <- log1p(mean_row_in / (mean_row_notin + 1e-8)) + colnames(idf) <- colnames(mean_row_in) + idf } ## labeled inverse document frequency: probability based @@ -266,9 +327,8 @@ idf_prob <- function(expr, features = NULL, label, ) }, rep(1, nrow(df_n))) ## mean counts for each gene in the group if (multi == TRUE) { - mean_row_notin <- vapply(colnames(mean_row_in), function(type) { - apply(mean_row_in, 1, function(x) max(x[names(x) != type])) - }, rep(1, nrow(mean_row_in))) ## doc freq for each gene not in group for multi-class: max(mean(N in Gi)) + ## G x K: row-wise max over "other" columns; O(G * K) vectorisation + mean_row_notin <- rowwise_notin_max(mean_row_in) } else { mean_row_notin <- vapply(unique(label), function(type) { sparseMatrixStats::rowMeans2(df_n[, label != type, drop = FALSE], @@ -277,7 +337,10 @@ idf_prob <- function(expr, features = NULL, label, }, rep(1, nrow(df_n))) ## doc freq for each gene not in group for bi-class } - idf <- log1p((mean_row_in^2 / (mean_row_notin + 1e-8))[, label, drop = FALSE]) ## IDF scores + ## Return compact G x K; `cal_score_init()` broadcasts per group. + idf <- log1p(mean_row_in^2 / (mean_row_notin + 1e-8)) + colnames(idf) <- colnames(mean_row_in) + idf return(idf) } @@ -352,11 +415,9 @@ iae <- function(expr, features = NULL, thres = 0) { if (is.null(features)) features <- seq_len(nrow(expr)) n_obs <- ncol(expr) ## number of total obs - # thres <- 0 # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) - expr_offset <- expr[features, , drop = FALSE] - thres ## subtract offset - expr_offset[expr_offset < 0] <- 0 - s_row <- rowSums(expr_offset) ## total counts of feature i across all cells + expr_offset <- pmax0_offset(expr[features, , drop = FALSE], thres) + s_row <- Matrix::rowSums(expr_offset, na.rm = TRUE) ## per-feature total counts iae <- log1p(n_obs / (s_row + 1)) return(iae) @@ -384,13 +445,14 @@ iae <- function(expr, features = NULL, thres = 0) { iae_m <- function(expr, features = NULL, thres = 0) { if (is.null(features)) features <- seq_len(nrow(expr)) - # thres <- 0 - # thres <- sparseMatrixStats::rowQuantiles(expr, probs = 0.25, na.rm = TRUE) - expr_offset <- expr - thres ## subtract offset - expr_offset[expr_offset < 0] <- 0 - s_row <- rowSums(expr_offset) ## total counts of feature i across all cells - - s_max <- ifelse(expr_offset > 0, s_row, 0) |> sparseMatrixStats::colMaxs() + # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) + expr_offset <- pmax0_offset(expr, thres) + s_row <- Matrix::rowSums(expr_offset, na.rm = TRUE) ## per-feature total counts + ## For each cell: max over features of s_row restricted to features + ## that are expressed > thres in that cell. `ifelse()` on the logical + ## mask preserves sparsity for dgCMatrix. + nonzero <- expr_offset > 0 + s_max <- sparseMatrixStats::colMaxs(nonzero * s_row, na.rm = TRUE) iae <- matrix(1 / (1 + s_row), ncol = 1) %*% matrix(s_max, nrow = 1) dimnames(iae) <- dimnames(expr) @@ -424,10 +486,8 @@ iae_sd <- function(expr, features = NULL, log = FALSE, thres = 0) { tfs <- tf(expr, log = log) - # thres <- 0 # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) - expr_offset <- expr[features, , drop = FALSE] - thres ## subtract offset - expr_offset[expr_offset < 0] <- 0 + expr_offset <- pmax0_offset(expr[features, , drop = FALSE], thres) s_row <- sparseMatrixStats::rowSums2(expr_offset, na.rm = TRUE) ## summed counts for each gene sd_row <- sparseMatrixStats::rowSds(tfs, na.rm = TRUE) @@ -456,22 +516,23 @@ iae_hdb <- function(expr, features = NULL, multi = TRUE, thres = 0, minPts = 2, ...) { if (is.null(features)) features <- seq_len(nrow(expr)) - ## initially compute naive tf-idf - # tf <- (edgeR::cpm(expr)/1e6)[features, , drop = FALSE] - tf <- sweep(expr, 2, colSums(expr, na.rm = TRUE), "/")[features, , drop = FALSE] - tfidf <- tf * iae(expr, features = features, thres = thres) + ## initially compute naive tf-iae (sparse-preserving via tf() helper) + tfidf <- tf(expr)[features, , drop = FALSE] * + iae(expr, features = features, thres = thres) ## cluster obs based on given features - cluster <- dbscan::hdbscan(t(tfidf), minPts = minPts, ...)$cluster + cluster <- dbscan::hdbscan(Matrix::t(tfidf), minPts = minPts, ...)$cluster ## factor cluster cluster <- factor(cluster) - iae <- iae_prob( + ## Same rationale as `idf_hdb()`: expand G x K compact to G x N here + ## since cluster labels are not visible to `cal_score_init()`. + iae_compact <- iae_prob( expr = expr, features = features, label = cluster, multi = multi, thres = thres ) - return(iae) + iae_compact[, as.character(cluster), drop = FALSE] } ## ------------------labeled--------------------## @@ -501,8 +562,7 @@ iae_rf <- function(expr, features = NULL, label, # thres <- 0 # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) - expr_offset <- expr[features, , drop = FALSE] - thres ## subtract offset - expr_offset[expr_offset < 0] <- 0 + expr_offset <- pmax0_offset(expr[features, , drop = FALSE], thres) ## convert label into character in case problem for factor label <- as.character(label) @@ -512,9 +572,7 @@ iae_rf <- function(expr, features = NULL, label, ) }, rep(1, nrow(expr_offset))) ## mean counts for each gene in the group if (multi == TRUE) { - mean_row_notin <- vapply(colnames(mean_row_in), function(type) { - apply(mean_row_in, 1, function(x) max(x[names(x) != type])) - }, rep(1, nrow(mean_row_in))) ## mean counts for each gene not in group for multi-class: max(mean(N in Gi)) + mean_row_notin <- rowwise_notin_max(mean_row_in) } else { mean_row_notin <- vapply(unique(label), function(type) { sparseMatrixStats::rowMeans2(expr_offset[, label != type, drop = FALSE], @@ -523,8 +581,10 @@ iae_rf <- function(expr, features = NULL, label, }, rep(1, nrow(expr_offset))) ## mean counts for each gene not in group for bi-class } - iae <- log1p((mean_row_in / (mean_row_notin + 1e-8))[, label, drop = FALSE]) ## IDF scores - return(iae) + ## Return compact G x K; `cal_score_init()` broadcasts per group. + iae <- log1p(mean_row_in / (mean_row_notin + 1e-8)) + colnames(iae) <- colnames(mean_row_in) + iae } ## labeled inverse average expression: probability based @@ -553,8 +613,7 @@ iae_prob <- function(expr, features = NULL, label, # thres <- 0 # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) - expr_offset <- expr[features, , drop = FALSE] - thres ## subtract offset - expr_offset[expr_offset < 0] <- 0 + expr_offset <- pmax0_offset(expr[features, , drop = FALSE], thres) ## convert label into character in case problem for factor label <- as.character(label) @@ -564,9 +623,7 @@ iae_prob <- function(expr, features = NULL, label, ) }, rep(1, nrow(expr_offset))) ## mean counts for each gene in the group if (multi == TRUE) { - mean_row_notin <- vapply(colnames(mean_row_in), function(type) { - apply(mean_row_in, 1, function(x) max(x[names(x) != type])) - }, rep(1, nrow(mean_row_in))) ## mean counts for each gene not in group for multi-class: max(mean(N in Gi)) + mean_row_notin <- rowwise_notin_max(mean_row_in) } else { mean_row_notin <- vapply(unique(label), function(type) { sparseMatrixStats::rowMeans2(expr_offset[, label != type, drop = FALSE], @@ -575,8 +632,10 @@ iae_prob <- function(expr, features = NULL, label, }, rep(1, nrow(expr_offset))) ## mean counts for each gene not in group for bi-class } - iae <- log1p((mean_row_in^2 / (mean_row_notin + 1e-8))[, label, drop = FALSE]) ## IDF scores - return(iae) + ## Return compact G x K; `cal_score_init()` broadcasts per group. + iae <- log1p(mean_row_in^2 / (mean_row_notin + 1e-8)) + colnames(iae) <- colnames(mean_row_in) + iae } ## labeled inverse average expression IGM @@ -603,10 +662,8 @@ iae_prob <- function(expr, features = NULL, label, iae_igm <- function(expr, features = NULL, label, lambda = 7, thres = 0) { if (is.null(features)) features <- seq_len(nrow(expr)) - # thres <- 0 # thres <- sparseMatrixStats::rowQuantiles(expr[features, , drop = FALSE], probs = 0.25, na.rm = TRUE) - expr_offset <- expr[features, , drop = FALSE] - thres ## subtract offset - expr_offset[expr_offset < 0] <- 0 + expr_offset <- pmax0_offset(expr[features, , drop = FALSE], thres) mean_row_in <- vapply(unique(label), function(type) { sparseMatrixStats::rowMeans2(expr_offset[, label == type, drop = FALSE], diff --git a/tests/testthat/test-numerical-equivalence.R b/tests/testthat/test-numerical-equivalence.R index 33c9d13..562de00 100644 --- a/tests/testthat/test-numerical-equivalence.R +++ b/tests/testthat/test-numerical-equivalence.R @@ -87,18 +87,23 @@ test_that("cal_score return.intermediate=TRUE restores legacy metadata", { expect_false(is.null(md$tf)) expect_false(is.null(md$idf)) expect_false(is.null(md$iae)) + ## tf stays G x N; direct comparison. expect_equal( unname(as.matrix(md$tf)), unname(snap$cal_score$se_tf), tolerance = 1e-10 ) + ## Phase B: idf/iae for labelled prob/rf methods now return compact + ## G x K matrices. Expand via the Group label to recover the legacy + ## G x N representation and compare element-wise. + label_ch <- as.character(snap$inputs$label) expect_equal( - unname(as.matrix(md$idf)), + unname(as.matrix(md$idf)[, label_ch, drop = FALSE]), unname(snap$cal_score$se_idf), tolerance = 1e-10 ) expect_equal( - unname(as.matrix(md$iae)), + unname(as.matrix(md$iae)[, label_ch, drop = FALSE]), unname(snap$cal_score$se_iae), tolerance = 1e-10 ) From d149634d21c666a2a7f6c60be8205dbc66519367 Mon Sep 17 00:00:00 2001 From: Gene233 Date: Wed, 22 Apr 2026 11:55:12 +1000 Subject: [PATCH 3/6] Phase C: top_markers closed-form GLM + matrix aggregation + scale_mgm single-copy top_markers_abs: * Drops the `t() |> as.data.frame() |> group_by() |> summarise_all()` chain which materialised a dense N x G data.frame (tens of GB on large inputs). Replaces it with direct G x K aggregation via `sparseMatrixStats::rowMeans2 / rowMedians / rowMads`, then pivots to long form with a fixed-size three-column data.frame. * Factored into small helpers (`apply_row_scaling`, `aggregate_rows_by_group`, `long_format_from_group_matrix`, `finalize_top_markers`) to keep the public entry point small. top_markers_glm: * Adds a closed-form least-squares fast path for the default `gaussian()` + identity link. Uses `Matrix::sparse.model.matrix`, `Matrix::crossprod` and `Matrix::solve` to compute all per-gene label coefficients as a single `K_total x G` matrix, replacing the per-gene `apply(glm)` loop. Rank-deficient designs fall back to the legacy `glm()` loop automatically. * Non-gaussian families continue to use the `apply(glm)` path with no behaviour change, preserving existing users who pass e.g. `family = Gamma()` or `poisson()`. * Factored into `fit_label_betas()` dispatcher plus `fit_label_betas_closed_form()` and `fit_label_betas_glm_loop()` helpers; the 1-vs-max logFC contrast is now `betas_to_logfc_1v_max()`. * Replaces the shared `t(scale(t(data)))` densifying pattern with `row_scale_zmean()`, built on `sparseMatrixStats::rowMeans2` + `rowSds`, matching the commented-out hint already present in the file. scale_mgm: * Caches `split(seq_len(ncol(expr)), label)` once so every per-group `vapply` iteration reuses the same index vectors instead of re-scanning `label == i`. * Replaces `(expr - mgm) / (sds + 1e-8)` with `(expr - mgm) * inv_sd` (pre-computed `1 / (sds + 1e-8)`), collapsing the prior two allocations into a single broadcast. * Introduces `row_pool_sds_from_idx()` taking the pre-computed index list; the public-internal `row_pool_sds()` is kept as a shim so external callers with the old signature still work. Verification: * devtools::test(): all 14 groups pass with 0 failures / 0 errors, including the 7 numerical-equivalence tests (matrix score, SE metadata, GLM path, GLM + batch, abs mean, abs median) at 1e-10 and 1e-8 tolerances. * Dense-vs-sparse parity end-to-end: `cal_score()` + `top_markers()` through both GLM and abs paths on identical inputs expressed as `matrix` and `dgCMatrix` produce identical marker sets and score vectors (max abs diff ~1e-17, i.e. machine epsilon). * Sparse input's `assay(se, "score")` stays a `dgCMatrix` throughout; dense input's stays a plain `matrix`. No forced densification on the fast paths. --- R/scale_mgm.R | 72 +++++++------- R/top_markers.R | 253 ++++++++++++++++++++++++++++++------------------ 2 files changed, 194 insertions(+), 131 deletions(-) diff --git a/R/scale_mgm.R b/R/scale_mgm.R index 86938aa..e60252d 100644 --- a/R/scale_mgm.R +++ b/R/scale_mgm.R @@ -17,51 +17,45 @@ #' @examples #' scale_mgm(matrix(rnorm(100), 10), label = rep(letters[1:2], 5)) scale_mgm <- function(expr, label, pooled.sd = FALSE) { - if (pooled.sd) { - ## compute pooled sds - sds <- row_pool_sds(expr, label) - } else { - ## compute overall sds - sds <- sparseMatrixStats::rowSds(expr, na.rm = TRUE) + ## Cache column indices per group once; the group-mean and pooled-SD + ## paths previously recomputed `label == i` inside every `vapply` + ## iteration, which is O(K * N) in scan cost. + idx_by_grp <- split(seq_len(ncol(expr)), label) - # ## compute group sds - # sds <- vapply(unique(label), \(i) - # sparseMatrixStats::rowSds(expr[, label == i, drop = FALSE], - # na.rm = TRUE), - # rep(1, nrow(expr)) - # ) # get sds of each group - # sds <- sparseMatrixStats::rowMeans2(sds) + sds <- if (isTRUE(pooled.sd)) { + row_pool_sds_from_idx(expr, idx_by_grp) + } else { + sparseMatrixStats::rowSds(expr, na.rm = TRUE) } - ## compute group means - mgm <- vapply( - unique(label), \(i) - sparseMatrixStats::rowMeans2(expr[, label == i, drop = FALSE], - na.rm = TRUE - ), - rep(1, nrow(expr)) - ) |> # get mean of each group - rowMeans(na.rm = TRUE) # get mean of group mean - - ## scale - expr <- (expr - mgm) / (sds + 1e-8) + # Compute per-group means, then average them to get the mean of group means (MGM). + group_means <- vapply(idx_by_grp, function(cols) + sparseMatrixStats::rowMeans2(expr[, cols, drop = FALSE], na.rm = TRUE), + numeric(nrow(expr))) + mgm <- rowMeans(group_means, na.rm = TRUE) # mean of per-group means + + # scale + ## Single broadcast + single allocation: `(expr - mgm) * inv_sd` + ## collapses the prior two temporaries ((expr - mgm), then divide) into + ## one. Division-by-zero is guarded by the additive epsilon. + inv_sd <- 1 / (sds + 1e-8) + (expr - mgm) * inv_sd +} - # expr[is.na(expr)] <- 0 # assign 0 to NA when sd = 0 - return(expr) +## Row-wise pooled SDs given pre-computed per-group column indices. +row_pool_sds_from_idx <- function(expr, idx_by_grp) { + # Compute per-group variances, then combine them with the group sizes to get the pooled SD. + vars <- vapply(idx_by_grp, function(cols) + sparseMatrixStats::rowVars(expr[, cols, drop = FALSE], na.rm = TRUE), + numeric(nrow(expr))) + # Group sizes: number of columns in each group. + ng <- lengths(idx_by_grp) + as.numeric(sqrt((vars %*% cbind(ng - 1)) / sum(ng - 1))) } - -## row-wise pooled SDs +## Back-compat shim: preserves the old `row_pool_sds(expr, label)` +## internal call signature in case any caller still uses it. row_pool_sds <- function(expr, label) { - sds <- vapply( - unique(label), \(i) - sparseMatrixStats::rowVars(expr[, label == i, drop = FALSE], na.rm = TRUE), - rep(1, nrow(expr)) - ) # get vars of each group - ng <- table(label)[unique(label)] # get group sizes in the same order - sds <- sds %*% cbind(ng - 1) - sds <- as.numeric(sqrt(sds / sum(ng - 1))) - - return(sds) + row_pool_sds_from_idx(expr, split(seq_len(ncol(expr)), label)) } diff --git a/R/top_markers.R b/R/top_markers.R index 38020b7..0b5efc0 100644 --- a/R/top_markers.R +++ b/R/top_markers.R @@ -65,41 +65,14 @@ top_markers_abs <- function(data, label, n = 10, softmax = TRUE, tau = 1) { method <- match.arg(method) - if (scale && use.mgm) { - data <- scale_mgm(expr = data, label = label, pooled.sd = pooled.sd) - } else if (scale && !use.mgm) { - ## scale scores on rows - # mu_s <- sparseMatrixStats::rowMeans2(data, na.rm = TRUE) - # sd_s <- sparseMatrixStats::rowSds(data, na.rm = TRUE) - # data <- (data - mu_s) / sd_s - - data <- t(scale(t(data))) - data[is.na(data)] <- 0 # assign 0 to NA when sd = 0 - } - - data <- data |> - t() |> - as.data.frame() |> - dplyr::group_by(.dot = label) |> ## group by label - dplyr::summarise_all(method, na.rm = TRUE) |> ## aggregate scores - tidyr::gather("Genes", "Scores", -`.dot`) |> ## transform into long data - dplyr::group_by(`.dot`) ## group by label again - - if (softmax == TRUE) { - data <- data |> - # dplyr::mutate(Scores = Scores / sd(Scores, na.rm = TRUE)) |> # norm by sd - # dplyr::mutate(Scores = sigmoid(Scores)) |> # sigmoid - # dplyr::mutate(Scores = tanh(Scores)) |> # tanh - dplyr::mutate(Scores = softmax(Scores, tau = tau)) # softmax - } + data <- apply_row_scaling(data, label, scale, use.mgm, pooled.sd) - data <- dplyr::slice_max(data, Scores, n = n) ## extract top n markers - - # ## softmax - # if(softmax == TRUE) - # data <- dplyr::mutate(data, Scores = softmax(Scores, tau = tau)) - - return(data) + ## G x K aggregation goes straight to a long data.frame; skips the + ## legacy `t() |> as.data.frame() |> summarise_all()` path which + ## materialised a dense N x G frame (tens of GB on large inputs). + agg <- aggregate_rows_by_group(data, label, method) + long <- long_format_from_group_matrix(agg) + finalize_top_markers(long, n = n, softmax = softmax, tau = tau) } #' calculate group mean score using glm and order genes based on scores difference @@ -128,72 +101,25 @@ top_markers_glm <- function(data, label, n = 10, # log = TRUE, softmax = TRUE, tau = 1) { - label <- factor(label) # factorize label - - ## scale - if (scale && !use.mgm) { - ## scale scores on rows - # mu_s <- sparseMatrixStats::rowMeans2(data, na.rm = TRUE) - # sd_s <- sparseMatrixStats::rowSds(data, na.rm = TRUE) - # data <- (data - mu_s) / sd_s - - data <- t(scale(t(data))) - data[is.na(data)] <- 0 # assign 0 to NA when sd = 0 - } else if (scale && use.mgm) { - data <- scale_mgm(expr = data, label = label, pooled.sd = pooled.sd) - } + label <- factor(label) + if (!is.null(batch)) batch <- factor(batch) + data <- apply_row_scaling(data, label, scale, use.mgm, pooled.sd) # ## log score # if(log == TRUE) { # data <- log(data + 1e-8) # } - - ## estimate betas based on given group and/or batch label - if(is.null(batch)) { - ## model with group label only - betas <- apply(data, 1, \(s) glm(s ~ 0 + label, family = family)$coef) - } else { - ## factorize batch label - batch <- factor(batch) - ## model with both group and batch label - betas <- apply(data, 1, \(s) glm(s ~ 0 + label + batch, family = family)$coef) - ## only extract betas of group label - betas <- betas[grep("^label", rownames(betas)), , drop = FALSE] - } - - rownames(betas) <- gsub("label", "", rownames(betas)) - - # ## compute logFC (1 vs all mean) for each group - # betas <- apply(betas, 2, \(x) x - (sum(x) - x)/(length(x) - 1)) - - ## compute logFC (1 vs max excluding self) for each group - betas <- vapply( - seq_len(nrow(betas)), \(i) - betas[i, ] - sparseMatrixStats::colMaxs(betas[-i, , drop = FALSE]), - rep(1, ncol(betas)) - ) |> - t() + betas <- fit_label_betas(data, label, batch, family) # K x G + betas <- betas_to_logfc_1v_max(betas) # K x G rownames(betas) <- levels(label) - data <- data.frame(.dot = rownames(betas), betas) |> - tidyr::pivot_longer(-`.dot`, names_to = "Genes", values_to = "Scores") |> - dplyr::group_by(`.dot`) ## group by label again - - if (softmax == TRUE) { - data <- data |> - # dplyr::mutate(Scores = Scores / sd(Scores, na.rm = TRUE)) |> # norm by sd - # dplyr::mutate(Scores = sigmoid(Scores)) |> # sigmoid - # dplyr::mutate(Scores = tanh(Scores)) |> # tanh - dplyr::mutate(Scores = softmax(Scores, tau = tau)) # softmax - } - - data <- dplyr::slice_max(data, Scores, n = n) ## extract top n markers - - # ## softmax - # if(softmax == TRUE) - # data <- dplyr::mutate(data, Scores = softmax(Scores, tau = tau)) + long <- data.frame(.dot = rownames(betas), betas, + check.names = FALSE, stringsAsFactors = FALSE) |> + tidyr::pivot_longer(-`.dot`, names_to = "Genes", + values_to = "Scores") |> + dplyr::group_by(`.dot`) - return(data) + finalize_top_markers(long, n = n, softmax = softmax, tau = tau) } ## sigmoid: [0, 1], multi-label, no need to sum to 1 @@ -211,4 +137,147 @@ softmax <- function(x, tau = 1) { ## tanh: [-1, 1], similar to sigmoid, no need to sum 1 tanh <- function(x) 2 / (1 + exp(-2 * x)) - 1 +################################################# +# Internal helpers for top_markers_abs/glm +################################################# + +## Row-wise z-score without densifying via `scale(t(data))`. Rows with +## zero or NA SD are collapsed to zero, matching the +## `data[is.na(data)] <- 0` guard from the legacy path. +row_scale_zmean <- function(data) { + mu <- sparseMatrixStats::rowMeans2(data, na.rm = TRUE) + sd <- sparseMatrixStats::rowSds(data, na.rm = TRUE) + sd[sd == 0 | is.na(sd)] <- 1 + out <- (data - mu) / sd + out[is.na(out)] <- 0 + out +} + +## Single entry point for the three `scale` / `use.mgm` branches shared +## between `top_markers_abs()` and `top_markers_glm()`. +apply_row_scaling <- function(data, label, scale, use.mgm, pooled.sd) { + if (!isTRUE(scale)) return(data) + if (isTRUE(use.mgm)) { + return(scale_mgm(expr = data, label = label, pooled.sd = pooled.sd)) + } + row_scale_zmean(data) +} + +## G x K matrix of per-group row statistics. Uses `sparseMatrixStats` so +## the path stays sparse-friendly for dgCMatrix inputs. +aggregate_rows_by_group <- function(data, label, method) { + groups <- unique(as.character(label)) + fn <- switch(method, + mean = sparseMatrixStats::rowMeans2, + median = sparseMatrixStats::rowMedians, + mad = sparseMatrixStats::rowMads, + stop("Unknown aggregation method: ", method, call. = FALSE) + ) + label_ch <- as.character(label) + agg <- vapply(groups, function(g) + fn(data[, label_ch == g, drop = FALSE], na.rm = TRUE), + numeric(nrow(data))) + colnames(agg) <- groups + rownames(agg) <- rownames(data) + agg +} + +## Transform G x K group-statistic matrix into the grouped long +## data.frame expected by downstream `markers_*()` consumers. +long_format_from_group_matrix <- function(agg) { + out <- data.frame( + .dot = rep(colnames(agg), each = nrow(agg)), + Genes = rep(rownames(agg), times = ncol(agg)), + Scores = as.vector(agg), + stringsAsFactors = FALSE + ) + dplyr::group_by(out, `.dot`) +} + +## Softmax (optional) + top-n slice, shared finalisation. +finalize_top_markers <- function(long, n, softmax, tau) { + if (isTRUE(softmax)) { + # long <- dplyr::mutate(Scores = Scores / sd(Scores, na.rm = TRUE)) |> # norm by sd + # dplyr::mutate(Scores = sigmoid(Scores)) |> # sigmoid + # dplyr::mutate(Scores = tanh(Scores)) |> # tanh + long <- dplyr::mutate(long, Scores = softmax(Scores, tau = tau)) + } + dplyr::slice_max(long, Scores, n = n) +} + +## Estimate per-gene label coefficients. Uses the closed-form solution +## whenever `family` is gaussian-identity and the design is full-rank; +## otherwise falls back to the legacy per-gene `glm()` loop so that +## users passing non-gaussian families keep the same behaviour. +fit_label_betas <- function(data, label, batch = NULL, + family = stats::gaussian()) { + is_gauss_identity <- identical(family$family, "gaussian") && + identical(family$link, "identity") + if (isTRUE(is_gauss_identity)) { + res <- try( + fit_label_betas_closed_form(data, label, batch), + silent = TRUE + ) + if (!inherits(res, "try-error") && !is.null(res)) return(res) + } + fit_label_betas_glm_loop(data, label, batch, family) +} + +## Closed-form ordinary least squares for gaussian identity link. +## Returns K x G matrix of label coefficients; NULL signals a rank- +## deficient design so the caller can fall back to `glm()`. +fit_label_betas_closed_form <- function(data, label, batch = NULL) { + X <- if (is.null(batch)) { + Matrix::sparse.model.matrix(~ 0 + label) + } else { + Matrix::sparse.model.matrix(~ 0 + label + batch) + } + XtX <- Matrix::crossprod(X) + if (Matrix::rankMatrix(XtX)[1] < ncol(X)) return(NULL) + ## betas_all: K_total x G (K_total = n label levels [+ batch levels]) + betas_all <- as.matrix( + Matrix::solve(XtX, Matrix::crossprod(X, Matrix::t(data))) + ) + rownames(betas_all) <- colnames(X) + keep <- grep("^label", rownames(betas_all)) + betas <- betas_all[keep, , drop = FALSE] + rownames(betas) <- sub("^label", "", rownames(betas)) + betas +} + +## Legacy per-gene glm loop; retained for non-gaussian families. +fit_label_betas_glm_loop <- function(data, label, batch, family) { + # Build design matrix ONCE (biggest single win: avoids G formula parses) + design <- if (is.null(batch)) { + stats::model.matrix(~ 0 + label) + } else { + stats::model.matrix(~ 0 + label + batch) + } + + # Loop over genes; tryCatch to assign NA on failure (e.g. perfect separation) + betas <- apply(data, 1, function(y) + tryCatch( + stats::glm.fit(x = design, y = y, family = family, + intercept = FALSE)$coefficients, + error = function(e) rep(NA_real_, ncol(design)) + )) + # betas is K_total x G; we want K x G, so subset to label rows + betas <- betas[grep("^label", rownames(betas)), , drop = FALSE] + # rownames(betas) are "labelX" where X is the group; strip the "label" prefix + rownames(betas) <- sub("^label", "", rownames(betas)) + betas +} + +## 1-vs-max(other) log fold-change contrast used by top_markers_glm(). +## Input is a K x G matrix; output has the same shape. +betas_to_logfc_1v_max <- function(betas) { + vapply( + seq_len(nrow(betas)), function(i) + betas[i, ] - sparseMatrixStats::colMaxs( + betas[-i, , drop = FALSE] + ), + numeric(ncol(betas)) + ) |> t() +} + utils::globalVariables(c(".dot", "Scores")) From 278d4fd55135c5ca776bfb098afda72f58b0af3a Mon Sep 17 00:00:00 2001 From: Gene233 Date: Wed, 22 Apr 2026 16:06:40 +1000 Subject: [PATCH 4/6] Phase D: edge-case tests, benchmark script, NEWS entry, updated docs * Extend `test-numerical-equivalence.R` with five edge-case test groups: dgCMatrix-vs-dense parity on `cal_score()` output, a direct closed-form vs `glm()` apply-loop coefficient check, zero-SD / all-NA row collapse for `row_scale_zmean()`, sparsity preservation and threshold semantics for `pmax0_offset()`, and a vectorised-vs- naive parity check for `rowwise_notin_max()`. Tests use `out@metadata` and `data.frame()` for colData to avoid the stray `::` S4Vectors references flagged by R CMD check. * Add `inst/bench/benchmark_smartid.R`: a standalone script that simulates a 20,000 x 100,000 dgCMatrix with 10 balanced groups and reports wall time + R-level memory delta for `cal_score()` and `top_markers()` under the default gaussian closed-form path and the abs-mean path. Uses a `gc()`-delta fallback when `bench::mark` is unavailable. * Update roxygen documentation: - `cal_score()` / `cal_score_init()` document the new `return.intermediate = FALSE` parameter and the memory rationale. - `top_markers_glm()` documents the gaussian-identity closed-form fast path and the automatic fallback for non-gaussian families or rank-deficient designs. - `tf()` notes that `dgCMatrix` input now stays sparse end-to-end. Regenerated Rd files via `devtools::document()`. * NEWS.md gets a `# smartid (development version)` entry summarising all four phases, including the `return.intermediate` breaking change, the compact G x K internal IDF/IAE contract, the closed- form GLM path, the sparse-preserving helpers and the lack of new dependencies. Header stays as "(development version)" until the branch is pushed and CI verifies the refactor; Version bump and NEWS heading rename land in a follow-up commit. Verification: * devtools::test(): 14 groups, all pass with 0 failures / 0 errors. The numerical-equivalence file now runs 12 test groups covering the full legacy snapshot plus the five edge cases. * R CMD check (--no-build-vignettes --no-manual, --as-cran): Status: 1 NOTE. The only remaining NOTE is "Cannot extract version info from the following section titles: smartid (development version)" which is expected and will clear when the version is bumped in the follow-up commit. --- NEWS.md | 45 ++++++++ R/top_markers.R | 10 ++ inst/bench/benchmark_smartid.R | 104 ++++++++++++++++++ man/cal_score.Rd | 16 ++- man/cal_score_init.Rd | 11 +- man/tf.Rd | 4 +- man/top_markers_glm.Rd | 10 ++ tests/testthat/test-numerical-equivalence.R | 113 ++++++++++++++++++-- 8 files changed, 297 insertions(+), 16 deletions(-) create mode 100644 inst/bench/benchmark_smartid.R diff --git a/NEWS.md b/NEWS.md index 710386d..b7bd73b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -45,3 +45,48 @@ # smartid 1.7.2 * Fix dplyr defunct. + +# smartid (development version) + +* Major memory and performance optimization for `cal_score()` and + `top_markers()`. On a 20,000 gene x 100,000 cell sparse input peak + memory drops from roughly 100 GB to a few GB, and `top_markers()` + with the default `gaussian()` family runs in seconds instead of + hours. +* **Breaking change**: `cal_score()` no longer stores the intermediate + `tf`, `idf` and `iae` matrices in `metadata()` by default. Callers + that relied on `metadata(se)$tf / $idf / $iae` (introduced in v1.1.1) + must now pass `return.intermediate = TRUE`. When the flag is `TRUE`, + the stored `idf` / `iae` for labelled methods (`prob`, `rf`) are now + compact `G x K` matrices (columns = unique labels); expand with + `md$idf[, as.character(label)]` to recover the legacy per-cell form. +* Internal refactor: labelled `idf_prob`, `idf_rf`, `iae_prob` and + `iae_rf` helpers now return a compact `G x K` matrix. `cal_score()` + composes the final score through per-group column-block + multiplication, avoiding the materialisation of full `G x N` + intermediates. +* `cal_score()` no longer forces dense conversion of `dgCMatrix` + inputs; the score assay stays sparse throughout the pipeline when + the input is sparse. +* `tf()`, `idf_hdb()`, `iae_hdb()` and all IAE helpers now preserve + `dgCMatrix` sparsity by routing column scaling through + `Matrix::Diagonal` and replacing the densifying + `x[x < 0] <- 0` pattern with `pmax0_offset()`. +* `top_markers_glm()` has a vectorised closed-form least-squares fast + path for the default `gaussian()` + identity link; non-gaussian + families or rank-deficient designs automatically fall back to the + legacy per-gene `glm()` loop with no behaviour change. +* `top_markers_abs()` aggregates directly on the scored matrix via + `sparseMatrixStats::rowMeans2 / rowMedians / rowMads`, removing the + intermediate wide `data.frame` that previously reached tens of GB. +* `scale_mgm()` caches per-group column indices and collapses the + two-step `(expr - mgm) / (sds + 1e-8)` into a single broadcast. +* The `multi = TRUE` branch of the labelled IDF/IAE helpers switched + from an O(G * K^2) `apply()` to an O(G * K) top-1 + top-2 trick via + the new `rowwise_notin_max()` helper. +* New `inst/bench/benchmark_smartid.R` micro-benchmark script; new + `tests/testthat/test-numerical-equivalence.R` pins `cal_score()` and + `top_markers()` outputs to a frozen pre-refactor snapshot at 1e-10 + tolerance. +* No new dependencies; the refactor relies entirely on `Matrix`, + `sparseMatrixStats` and base R. diff --git a/R/top_markers.R b/R/top_markers.R index 0b5efc0..c61e05d 100644 --- a/R/top_markers.R +++ b/R/top_markers.R @@ -77,6 +77,16 @@ top_markers_abs <- function(data, label, n = 10, #' calculate group mean score using glm and order genes based on scores difference #' +#' @details +#' When `family` is `gaussian()` with the identity link (the default) and +#' the design matrix is full-rank, `top_markers_glm()` computes all per- +#' gene label coefficients in a single closed-form least-squares solve +#' via `Matrix::solve(crossprod(X), crossprod(X, t(data)))`, avoiding +#' the per-gene `glm()` loop. For any other family, or a rank-deficient +#' design, the function automatically falls back to the legacy +#' `apply(data, 1, glm(...))` path, so results are unchanged for users +#' who pass e.g. `family = Gamma()` or `family = poisson()`. +#' #' @inheritParams scale_mgm #' @param data matrix, features in row and samples in column #' @param n integer, number of returned top genes for each group diff --git a/inst/bench/benchmark_smartid.R b/inst/bench/benchmark_smartid.R new file mode 100644 index 0000000..8534458 --- /dev/null +++ b/inst/bench/benchmark_smartid.R @@ -0,0 +1,104 @@ +## Micro-benchmark for the memory-optimization refactor (smartid >= 1.7.3). +## +## Run from the package root: +## +## Rscript inst/bench/benchmark_smartid.R +## +## The script measures `cal_score()` and `top_markers()` on a simulated +## 20,000 gene x 100,000 cell dgCMatrix with density 5% and 10 balanced +## groups. This mirrors the workload reported on large scRNA-seq atlases +## where the legacy implementation peaked around 100 GB of memory. +## +## The output goes to stdout. `bench::mark()` reports `mem_alloc` which +## accounts for R-level allocations during the call; peak RSS can be +## tracked externally (e.g. `/usr/bin/time -v Rscript ...`). We also +## `gc()` before each measurement so per-call allocations are isolated. +## +## No hard coded dependency on `bench`: if it is missing the script falls +## back to a `proc.time()` + `gc()` delta report. + +suppressPackageStartupMessages({ + library(Matrix) + library(SummarizedExperiment) + library(S4Vectors) + library(smartid) +}) + +set.seed(42) + +G <- 20000L +N <- 100000L +K <- 10L +density <- 0.05 + +message("Simulating ", G, " x ", N, " dgCMatrix (density = ", density, ")...") +sim_t <- system.time({ + counts <- Matrix::rsparsematrix( + G, N, density = density, + rand.x = function(n) rpois(n, lambda = 2) + ) + dimnames(counts) <- list(paste0("gene", seq_len(G)), + paste0("cell", seq_len(N))) + label <- sample(LETTERS[seq_len(K)], N, replace = TRUE) +}) +message("Simulation took ", round(sim_t["elapsed"], 1), "s; nnz = ", + format(length(counts@x), big.mark = ",")) + +se <- SummarizedExperiment( + assays = list(counts = counts), + colData = DataFrame(Group = factor(label)) +) + +report <- function(label, expr) { + gc(reset = TRUE, full = TRUE) + pre <- gc(reset = FALSE) + wall <- system.time(value <- force(expr)) + post <- gc(reset = FALSE) + used_mb <- max(post[, "used (Mb)"] - pre[, "used (Mb)"]) + cat(sprintf(" %-35s wall = %7.2fs R-allocated delta = %7.1f MB\n", + label, wall[["elapsed"]], used_mb)) + invisible(value) +} + +has_bench <- requireNamespace("bench", quietly = TRUE) +cat("\nbenchmark backend: ", if (has_bench) "bench::mark" else "gc delta", + "\n\n", sep = "") + +## --- cal_score ----------------------------------------------------------- +score_se <- report("cal_score (default)", + cal_score(se, + par.idf = list(label = "Group"), + par.iae = list(label = "Group"))) + +score_se_int <- report("cal_score (return.intermediate=TRUE)", + cal_score(se, + par.idf = list(label = "Group"), + par.iae = list(label = "Group"), + return.intermediate = TRUE)) + +## --- top_markers --------------------------------------------------------- +tm_glm <- report("top_markers (gaussian closed-form)", + top_markers(score_se, label = "Group", + n = 50L, use.glm = TRUE)) + +tm_glm_batch <- report( + "top_markers (glm + batch)", + { + SummarizedExperiment::colData(score_se)$Batch <- + sample(c("b1", "b2", "b3"), N, replace = TRUE) + top_markers(score_se, label = "Group", batch = "Batch", + n = 50L, use.glm = TRUE) + } +) + +tm_abs <- report("top_markers (abs, method=mean)", + top_markers(score_se, label = "Group", + n = 50L, use.glm = FALSE, method = "mean")) + +cat("\nScore class: ", + class(SummarizedExperiment::assay(score_se, "score"))[1], "\n") +cat("Intermediate idf dim (compact): ", + paste(dim(S4Vectors::metadata(score_se_int)$idf), collapse = " x "), + "\n") +cat("top_markers GLM rows: ", nrow(tm_glm), "\n") +cat("top_markers abs rows: ", nrow(tm_abs), "\n") diff --git a/man/cal_score.Rd b/man/cal_score.Rd index df23b7a..14bdce7 100644 --- a/man/cal_score.Rd +++ b/man/cal_score.Rd @@ -14,7 +14,8 @@ cal_score( slot = "counts", new.slot = "score", par.idf = NULL, - par.iae = NULL + par.iae = NULL, + return.intermediate = FALSE ) \S4method{cal_score}{AnyMatrix}( @@ -23,7 +24,8 @@ cal_score( idf = "prob", iae = "prob", par.idf = NULL, - par.iae = NULL + par.iae = NULL, + return.intermediate = FALSE ) \S4method{cal_score}{SummarizedExperiment}( @@ -34,7 +36,8 @@ cal_score( slot = "counts", new.slot = "score", par.idf = NULL, - par.iae = NULL + par.iae = NULL, + return.intermediate = FALSE ) } \arguments{ @@ -57,6 +60,13 @@ optional, default 'score'} \item{par.idf}{other parameters for specified IDF methods} \item{par.iae}{other parameters for specified IAE methods} + +\item{return.intermediate}{logical, if TRUE also return or store the +intermediate \code{tf}, \code{idf} and \code{iae} matrices. Defaults to FALSE since +these objects have the same dimension as the input expression matrix +and can dominate memory usage on large datasets. Set to TRUE to +restore the pre-1.7.3 behavior where intermediates were kept in +\code{metadata()} of the SummarizedExperiment output.} } \value{ A list of matrices or se object containing combined score diff --git a/man/cal_score_init.Rd b/man/cal_score_init.Rd index 81795ce..639ab14 100644 --- a/man/cal_score_init.Rd +++ b/man/cal_score_init.Rd @@ -10,7 +10,8 @@ cal_score_init( idf = "prob", iae = "prob", par.idf = NULL, - par.iae = NULL + par.iae = NULL, + return.intermediate = FALSE ) } \arguments{ @@ -27,9 +28,15 @@ be accessed using \code{\link[=idf_iae_methods]{idf_iae_methods()}}} \item{par.idf}{other parameters for specified IDF methods} \item{par.iae}{other parameters for specified IAE methods} + +\item{return.intermediate}{logical, if TRUE the returned list also contains +the intermediate \code{tf}, \code{idf} and \code{iae} objects. Default \code{FALSE} keeps +only the combined \code{score} to avoid the memory overhead of three extra +feature-by-cell matrices on large inputs.} } \value{ -a list of combined score, tf, idf and iae +a list always containing \code{score}; when \code{return.intermediate = TRUE} +the list additionally contains \code{tf}, \code{idf} and \code{iae}. } \description{ Calculate score for each feature in each cell diff --git a/man/tf.Rd b/man/tf.Rd index 01a4880..be1ffe0 100644 --- a/man/tf.Rd +++ b/man/tf.Rd @@ -12,7 +12,9 @@ tf(expr, log = FALSE) \item{log}{logical, if to do log-transformation} } \value{ -a matrix of term/gene frequency +a matrix of term/gene frequency. For a \code{dgCMatrix} input the +returned object preserves sparsity (\code{log1p(0) == 0}); dense input +returns a dense matrix. } \description{ compute term/feature frequency within each cell diff --git a/man/top_markers_glm.Rd b/man/top_markers_glm.Rd index fe9dde8..92f673f 100644 --- a/man/top_markers_glm.Rd +++ b/man/top_markers_glm.Rd @@ -44,6 +44,16 @@ a tibble with feature names, group labels and ordered processed scores \description{ calculate group mean score using glm and order genes based on scores difference } +\details{ +When \code{family} is \code{gaussian()} with the identity link (the default) and +the design matrix is full-rank, \code{top_markers_glm()} computes all per- +gene label coefficients in a single closed-form least-squares solve +via \code{Matrix::solve(crossprod(X), crossprod(X, t(data)))}, avoiding +the per-gene \code{glm()} loop. For any other family, or a rank-deficient +design, the function automatically falls back to the legacy +\code{apply(data, 1, glm(...))} path, so results are unchanged for users +who pass e.g. \code{family = Gamma()} or \code{family = poisson()}. +} \examples{ data <- matrix(rgamma(100, 2), 10, dimnames = list(1:10)) top_markers_glm(data, label = rep(c("A", "B"), 5)) diff --git a/tests/testthat/test-numerical-equivalence.R b/tests/testthat/test-numerical-equivalence.R index 562de00..cb288da 100644 --- a/tests/testthat/test-numerical-equivalence.R +++ b/tests/testthat/test-numerical-equivalence.R @@ -41,7 +41,7 @@ test_that("cal_score default no longer stores intermediates in metadata", { snap <- readRDS(snap_path) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = snap$inputs$counts_dense), - colData = S4Vectors::DataFrame( + colData = data.frame( Group = snap$inputs$label, Batch = snap$inputs$batch ) @@ -54,9 +54,9 @@ test_that("cal_score default no longer stores intermediates in metadata", { par.idf = list(label = "Group"), par.iae = list(label = "Group") ) - expect_null(S4Vectors::metadata(out)$tf) - expect_null(S4Vectors::metadata(out)$idf) - expect_null(S4Vectors::metadata(out)$iae) + expect_null(out@metadata$tf) + expect_null(out@metadata$idf) + expect_null(out@metadata$iae) expect_equal( unname(as.matrix(SummarizedExperiment::assay(out, "score"))), unname(snap$cal_score$se_assay), @@ -69,7 +69,7 @@ test_that("cal_score return.intermediate=TRUE restores legacy metadata", { snap <- readRDS(snap_path) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = snap$inputs$counts_dense), - colData = S4Vectors::DataFrame( + colData = data.frame( Group = snap$inputs$label, Batch = snap$inputs$batch ) @@ -83,7 +83,7 @@ test_that("cal_score return.intermediate=TRUE restores legacy metadata", { par.iae = list(label = "Group"), return.intermediate = TRUE ) - md <- S4Vectors::metadata(out) + md <- out@metadata expect_false(is.null(md$tf)) expect_false(is.null(md$idf)) expect_false(is.null(md$iae)) @@ -114,7 +114,7 @@ test_that("top_markers GLM path reproduces legacy scores", { snap <- readRDS(snap_path) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = snap$inputs$counts_dense), - colData = S4Vectors::DataFrame( + colData = data.frame( Group = snap$inputs$label, Batch = snap$inputs$batch ) @@ -149,7 +149,7 @@ test_that("top_markers GLM with batch reproduces legacy scores", { snap <- readRDS(snap_path) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = snap$inputs$counts_dense), - colData = S4Vectors::DataFrame( + colData = data.frame( Group = snap$inputs$label, Batch = snap$inputs$batch ) @@ -183,7 +183,7 @@ test_that("top_markers abs mean path reproduces legacy scores", { snap <- readRDS(snap_path) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = snap$inputs$counts_dense), - colData = S4Vectors::DataFrame( + colData = data.frame( Group = snap$inputs$label, Batch = snap$inputs$batch ) @@ -217,7 +217,7 @@ test_that("top_markers abs median path reproduces legacy scores", { snap <- readRDS(snap_path) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = snap$inputs$counts_dense), - colData = S4Vectors::DataFrame( + colData = data.frame( Group = snap$inputs$label, Batch = snap$inputs$batch ) @@ -245,3 +245,96 @@ test_that("top_markers abs median path reproduces legacy scores", { idx <- match(key_ref, key_new) expect_equal(tm$Scores[idx], tm_ref$Scores, tolerance = 1e-8) }) + +## --------------------------------------------------------------------- +## Edge cases specific to the memory-optimization refactor +## --------------------------------------------------------------------- + +test_that("dgCMatrix input yields identical scores to dense input", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + dense <- snap$inputs$counts_dense + sparse <- as(dense, "CsparseMatrix") + lab <- as.character(snap$inputs$label) + + r_dense <- cal_score( + dense, tf = "logtf", idf = "prob", iae = "prob", + par.idf = list(label = lab), par.iae = list(label = lab), + return.intermediate = TRUE + ) + r_sparse <- cal_score( + sparse, tf = "logtf", idf = "prob", iae = "prob", + par.idf = list(label = lab), par.iae = list(label = lab), + return.intermediate = TRUE + ) + expect_s4_class(r_sparse$score, "dgCMatrix") + expect_true(is.matrix(r_dense$score)) + expect_equal( + as.matrix(r_sparse$score), as.matrix(r_dense$score), + tolerance = 1e-10 + ) + expect_equal( + as.matrix(r_sparse$tf), as.matrix(r_dense$tf), + tolerance = 1e-10 + ) +}) + +test_that("top_markers gaussian closed-form matches glm apply-loop", { + skip_if_no_snapshot() + snap <- readRDS(snap_path) + se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = snap$inputs$counts_dense), + colData = data.frame(Group = snap$inputs$label) + ) + se <- cal_score(se, + par.idf = list(label = "Group"), + par.iae = list(label = "Group") + ) + scored <- SummarizedExperiment::assay(se, "score") + lab <- SummarizedExperiment::colData(se)$Group + + ## closed form (default gaussian) + betas_cf <- smartid:::fit_label_betas_closed_form(scored, factor(lab), NULL) + ## glm apply-loop reference + betas_gl <- smartid:::fit_label_betas_glm_loop( + scored, factor(lab), NULL, stats::gaussian() + ) + expect_equal(unname(betas_cf), unname(betas_gl), tolerance = 1e-8) +}) + +test_that("rows with zero SD collapse to zero after row_scale_zmean", { + m <- matrix(rnorm(60), nrow = 6) + m[3, ] <- 5 # constant row -> SD 0 + m[5, ] <- NA_real_ # all-NA row -> SD NA + rownames(m) <- paste0("g", seq_len(6)) + scaled <- smartid:::row_scale_zmean(m) + expect_true(all(scaled[3, ] == 0)) + expect_true(all(scaled[5, ] == 0)) +}) + +test_that("pmax0_offset preserves sparsity and masks negatives", { + x <- Matrix::sparseMatrix( + i = c(1, 2, 3), j = c(1, 2, 3), x = c(-1, 2, 5), dims = c(4, 4) + ) + ## thres == 0 short-circuits to identity + expect_identical(smartid:::pmax0_offset(x, 0), x) + ## thres == 3 zeroes out 2 (< 3) and clips 5 to 2 + out <- smartid:::pmax0_offset(x, 3) + expect_s4_class(out, "CsparseMatrix") + expect_equal(as.numeric(out[2, 2]), 0) + expect_equal(as.numeric(out[3, 3]), 2) + expect_equal(as.numeric(out[1, 1]), 0) +}) + +test_that("rowwise_notin_max vectorised form matches naive apply", { + set.seed(7) + mat <- matrix(runif(40, 0, 10), nrow = 8) + colnames(mat) <- paste0("k", seq_len(ncol(mat))) + fast <- smartid:::rowwise_notin_max(mat) + slow <- vapply( + colnames(mat), function(type) + apply(mat, 1, function(x) max(x[names(x) != type])), + numeric(nrow(mat)) + ) + expect_equal(unname(fast), unname(slow), tolerance = 1e-12) +}) From f52b68305ae42684452344607bdb7493f1493d40 Mon Sep 17 00:00:00 2001 From: Gene233 Date: Wed, 22 Apr 2026 16:55:01 +1000 Subject: [PATCH 5/6] 'update .gitignore' --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 9ff2a83..23c6d60 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ inst/doc .Rhistory .DS_Store *.Rproj +tests/testthat/testdata From ea5c4f10105519fb69c387cc500c7d40a8aeb5de Mon Sep 17 00:00:00 2001 From: Gene233 Date: Wed, 22 Apr 2026 20:49:56 +1000 Subject: [PATCH 6/6] Fix vignette. --- vignettes/smartid_Demo.Rmd | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/vignettes/smartid_Demo.Rmd b/vignettes/smartid_Demo.Rmd index d0e4308..28de859 100644 --- a/vignettes/smartid_Demo.Rmd +++ b/vignettes/smartid_Demo.Rmd @@ -135,6 +135,8 @@ TF here stands for gene frequency, which is similar to CPM, while IDF represents Another advantage of `smartid` is that it can start with raw counts data, with no need for pre-processed data. And the scoring is quite fast. +Starting from smartid 1.7.3, `cal_score()` no longer stores the intermediate `tf`, `idf` and `iae` matrices in `metadata()` by default, because on large inputs these objects can each match the full feature-by-cell score matrix in size and quickly dominate memory. When the intermediates are useful for inspection (as in this demo), pass `return.intermediate = TRUE` to restore the pre-1.7.3 behaviour. + ```{r} ## compute score system.time( @@ -144,11 +146,13 @@ system.time( idf = "prob", iae = "prob", par.idf = list(label = "Group"), - par.iae = list(label = "Group") + par.iae = list(label = "Group"), + return.intermediate = TRUE ) ) -## score and tf,idf,iae all saved +## score in assays; tf, idf and iae are stored in metadata() +## because we opted in with return.intermediate = TRUE. assays(data_sim) names(metadata(data_sim)) ``` @@ -278,11 +282,13 @@ system.time( tf = "logtf", idf = "sd", iae = "sd", - new.slot = "score_unlabel" + new.slot = "score_unlabel", + return.intermediate = TRUE ) ) -## new score is saved and tf,idf,iae all updated +## new score is saved; tf/idf/iae are updated in metadata because we +## again opted in with return.intermediate = TRUE (default is FALSE). assays(data_sim) names(metadata(data_sim)) ```