Better integrated the new multithreaded calculations into the other parts.

Author: pcarbo

DESCRIPTION

@@ -1,7 +1,7 @@
 Encoding: UTF-8
 Type: Package
 Package: fastTopics
-Version: 0.7-41
+Version: 0.7-42
 Date: 2026-04-17
 Title: Fast Algorithms for Fitting Topic Models and Non-Negative
     Matrix Factorizations to Count Data

R/fit_poisson_nmf.R

@@ -455,7 +455,7 @@ fit_poisson_nmf <- function (X, k, fit0, numiter = 100,
     cat(sprintf("Running at most %d %s updates, %s extrapolation ",
                 numiter,method.text,
                 ifelse(control$extrapolate,"with","without")))
-    cat("(fastTopics 0.7-41).\n")
+    cat("(fastTopics 0.7-42).\n")
   }
 
   # INITIALIZE ESTIMATES
@@ -554,7 +554,9 @@ fit_poisson_nmf_main_loop <- function (X, fit, numiter, update.factors,
     progress[i,"loglik.multinom"] <-
       loglik.const - fit$loss - sum(loglik_size_factors(X,fit$F,fit$L))
     progress[i,"dev"]         <- dev.const + 2*fit$loss
-    res <- with(poisson_nmf_kkt(X,fit$F,fit$L),
+    res <- with(poisson_nmf_kkt(X,fit$F,fit$L,
+                                version = ifelse(control$nc == 1,
+                                                 "Rcpp","RcppParallel")),
                 max(abs(rbind(F[update.factors,],
                               L[update.loadings,]))))
     progress[i,"res"]         <- res
@@ -647,7 +649,9 @@ update_poisson_nmf <- function (X, fit, update.factors, update.loadings,
 
   # Compute the value of the objective ("loss") function at the updated
   # estimates.
-  fit$loss      <- sum(cost(X,fit$L,t(fit$F),control$eps))
+  fit$loss <- sum(cost(X,fit$L,t(fit$F),control$eps,
+                       version = ifelse(control$nc == 1,
+                                        "Rcpp","RcppParallel")))
   fit$loss.fnly <- fit$loss
 
   # Output the updated "fit".
@@ -697,7 +701,9 @@ update_poisson_nmf_extrapolated <- function (X, fit, update.factors,
   # Compute the value of the objective (loss) function at the
   # extrapolated solution for the loadings (Ly) and the
   # non-extrapolated solution for the factors (Fn).
-  fit$loss.fnly <- sum(cost(X,fit$Ly,t(Fn),control$eps))
+  fit$loss.fnly <- sum(cost(X,fit$Ly,t(Fn),control$eps,
+                            version = ifelse(control$nc == 1,
+                                             "Rcpp","RcppParallel")))
 
   # Update the extrapolation parameters following Algorithm 3 of
   # Ang & Gillis (2019).

R/init_poisson_nmf.R

@@ -142,8 +142,9 @@ init_poisson_nmf <-
   L <- pmax(L,control$minval)
 
   # Compute the value of the objective ("loss") function at the
-  # initial estimates of the factors and loading.
-  loss <- sum(cost(X,L,t(F),control$eps))
+  # initial estimates of the factors and loadings.
+  loss <- sum(cost(X,L,t(F),control$eps,
+                   version = ifelse(control$nc == 1,"Rcpp","RcppParallel")))
 
   # Restore the BLAS settings.
   blas_set_num_threads(ncb)

R/likelihood.R

@@ -95,13 +95,13 @@ cost <- function (X, A, B, e = 1e-8, family = c("poisson","multinom"),
     B <- matrix(B,1,ncol(X))
 
   # Check and process "model" and "version" input arguments.
-  family  <- match.arg(family)
+  family <- match.arg(family)
   poisson <- family == "poisson"
   if (missing(version)) {
     if (is.matrix(X))
       version <- "R"
     else
-      version <- "Rcpp_parallel"
+      version <- "Rcpp"
   }
 
   # Compute the terms in the log-likelihoods that depend on A or B.
@@ -184,6 +184,7 @@ deviance_poisson_const <- function (X) {
 # Compute the residuals of the first-order Karush-Kuhn-Tucker (KKT)
 # conditions for Poisson non-negative matrix factorization at solution
 # estimate (F,L).
+#
 #' @importFrom Matrix sparseMatrix
 poisson_nmf_kkt <- function (X, F, L, e = 1e-8,
                              version = c("Rcpp_parallel","Rcpp")) {
@@ -198,10 +199,10 @@ poisson_nmf_kkt <- function (X, F, L, e = 1e-8,
     A <- sparseMatrix(i = d$i,j = d$j,x = y,dims = dim(X))
     return(list(F = F*(repmat(colSums(L),ncol(X)) - as.matrix(t(A) %*% L)),
                 L = L*(repmat(colSums(F),nrow(X)) - as.matrix(A %*% F))))
-  } else {
-    result <- poisson_nmf_kkt_sparse_parallel_rcpp(X,L,F,e)
-    return(list(F = F*(repmat(colSums(L),ncol(X)) - result$tAL),
-                L = L*(repmat(colSums(F),nrow(X)) - result$AF)))
+  } else if (version == "Rcpp_parallel") {
+    res <- poisson_nmf_kkt_sparse_parallel_rcpp(X,L,F,e)
+    return(list(F = F*(repmat(colSums(L),ncol(X)) - res$tAL),
+                L = L*(repmat(colSums(F),nrow(X)) - res$AF)))
   }
 }
 

tests/testthat/test_likelihood.R

@@ -9,10 +9,10 @@ test_that(paste("R, Rcpp and Rcpp_parallel versions of cost function return",
 
     # Generate a data set.
     out <- simulate_count_data(10,8,k)
-    X   <- out$X
-    F   <- out$F
-    L   <- out$L
-    Y   <- as(X,"CsparseMatrix")
+    X <- out$X
+    F <- out$F
+    L <- out$L
+    Y <- as(X,"CsparseMatrix")
 
     # Compute the loss function.
     f1 <- cost(X,L,t(F),version = "R")
@@ -23,11 +23,11 @@ test_that(paste("R, Rcpp and Rcpp_parallel versions of cost function return",
     f6 <- cost(Y,L,t(F),version = "Rcpp_parallel")
 
     # The cost function calculations should all give the same result.
-    expect_equal(f1,f2)
-    expect_equal(f1,f3)
-    expect_equal(f1,f4)
-    expect_equal(f1,f5)
-    expect_equal(f1,f6)
+    expect_equal(f1,f2,scale = 1,tolerance = 1e-10)
+    expect_equal(f1,f3,scale = 1,tolerance = 1e-10)
+    expect_equal(f1,f4,scale = 1,tolerance = 1e-10)
+    expect_equal(f1,f5,scale = 1,tolerance = 1e-10)
+    expect_equal(f1,f6,scale = 1,tolerance = 1e-10)
   }
 })
 
@@ -40,19 +40,20 @@ test_that(paste("loglik_poisson_nmf gives correct result for sparse and",
 
     # Generate a data set.
     out <- simulate_count_data(10,8,k)
-    X   <- out$X
+    X <- out$X
+    Y <- as(X,"CsparseMatrix")
     fit <- out[c("F","L")]
     class(fit) <- c("poisson_nmf_fit","list")
-    
+
     # Compute the log-likelikhood.
     f1 <- loglik_poisson_nmf_with_dpois(X,fit)
     f2 <- loglik_poisson_nmf(X,fit,e = 0)
-    f3 <- loglik_poisson_nmf(as(X,"CsparseMatrix"),fit,e = 0)
+    f3 <- loglik_poisson_nmf(Y,fit,e = 0)
     names(f1) <- rownames(X)
 
     # The likelihood calculations should all be the same.
-    expect_equal(f1,f2)
-    expect_equal(f1,f3)
+    expect_equal(f1,f2,scale = 1,tolerance = 1e-10)
+    expect_equal(f1,f3,scale = 1,tolerance = 1e-10)
   }
 })
 

tests/testthat/test_predict.R

@@ -17,7 +17,7 @@ test_that("project_poisson_nmf leaves the F matrix unchanged",{
 
   # Predict L in unseen (test) data points.
   capture.output(out <- project_poisson_nmf(test,F,numiter = 20))
-  expect_equal(rep(1,3),unname(diag(cor(out$F,F))),tolerance = 1e-8)
+  expect_equal(rep(1,3),unname(diag(cor(out$F,F))),scale=1,tolerance=1e-8)
 })
 
 test_that(paste("Running project_poisson_nmf on the training data",
@@ -32,5 +32,6 @@ test_that(paste("Running project_poisson_nmf on the training data",
   capture.output(fit <- init_poisson_nmf(X,F = dat$F,init.method = "random"))
   capture.output(fit <- fit_poisson_nmf(X,fit0 = fit))
   capture.output(out <- project_poisson_nmf(X,fit$F,numiter = 20))
-  expect_equal(rep(1,3),unname(diag(cor(out$L,fit$L))),tolerance = 1e-3)
+  expect_equal(rep(1,3),unname(diag(cor(out$L,fit$L))),scale=1,
+               tolerance=1e-3)
 })