add parallel AF assignment function to speed up get.af.spectra() in AutoSpectral

Author: DrCytometer

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(assign.af.fluor.fast)
 export(fast_kde2d_cpp)
 export(find_local_maxima)
 export(poisson_irls_rcpp_parallel)

R/RcppExports.R

@@ -1,6 +1,10 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+assign_af_fluor <- function(raw_data, spectra, af_spectra, n_threads = 4L) {
+    .Call(`_AutoSpectralRcpp_assign_af_fluor`, raw_data, spectra, af_spectra, n_threads)
+}
+
 #' Fast Kernel Density Estimation in 2D
 #'
 #' @param x Numeric vector of x-values (typically FSC-A)

R/assign_af_fluor_fast.R

@@ -0,0 +1,43 @@
+# assign_af_fluor_fast.r
+
+#' @title Assign AF Spectrum By Fluorophore Projection Fast
+#'
+#' @description
+#' Projects the autofluorescence spectral variants into fluorophore unmixed
+#' space to determine which best fits each cell (event). A fast approximation for
+#' brute force sequential unmixing method in early versions of AutoSpectral.
+#' Provides essentially identical results to minimization of fluorophore signal
+#' (dist0 method). Substantially faster. Use L1 (absolute value) minimization,
+#' which works better than the standard L2 (squared error). Wrapper for the C++
+#' version, which uses OpenMP for parallel processing.
+#'
+#' @param raw.data Expression data from raw fcs files. Cells in rows and
+#' detectors in columns. Columns should be fluorescent data only and must
+#' match the columns in spectra.
+#' @param spectra Spectral signatures of fluorophores, normalized between 0
+#' and 1, with fluorophores in rows and detectors in columns.
+#' @param af.spectra Spectral signatures of autofluorescences, normalized
+#' between 0 and 1, with fluorophores in rows and detectors in columns. Prepare
+#' using `get.af.spectra`.
+#' @param threads Numeric, default is `1`.
+#' 
+#' @return The indices of the best-fitting AF spectrum per cell
+#' 
+#' @export
+
+assign.af.fluor.fast <- function(raw.data, spectra, af.spectra, threads = 1) {
+  # remove AF if present
+  if ("AF" %in% rownames(spectra)) 
+    spectra <- spectra[rownames(spectra) != "AF", , drop = FALSE]
+  
+  # Call C++ function
+  # Ensure inputs are matrices
+  af.idx <- assign_af_fluor(
+    as.matrix(raw.data), 
+    as.matrix(spectra), 
+    as.matrix(af.spectra), 
+    n_threads = threads
+  )
+  
+  return(af.idx)
+}

man/assign.af.fluor.fast.Rd

@@ -11,6 +11,20 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
+// assign_af_fluor
+Rcpp::IntegerVector assign_af_fluor(const arma::mat& raw_data, const arma::mat& spectra, const arma::mat& af_spectra, int n_threads);
+RcppExport SEXP _AutoSpectralRcpp_assign_af_fluor(SEXP raw_dataSEXP, SEXP spectraSEXP, SEXP af_spectraSEXP, SEXP n_threadsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type raw_data(raw_dataSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type spectra(spectraSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type af_spectra(af_spectraSEXP);
+    Rcpp::traits::input_parameter< int >::type n_threads(n_threadsSEXP);
+    rcpp_result_gen = Rcpp::wrap(assign_af_fluor(raw_data, spectra, af_spectra, n_threads));
+    return rcpp_result_gen;
+END_RCPP
+}
 // fast_kde2d_cpp
 List fast_kde2d_cpp(NumericVector x, NumericVector y, int n, NumericVector h, NumericVector x_limits, NumericVector y_limits);
 RcppExport SEXP _AutoSpectralRcpp_fast_kde2d_cpp(SEXP xSEXP, SEXP ySEXP, SEXP nSEXP, SEXP hSEXP, SEXP x_limitsSEXP, SEXP y_limitsSEXP) {
@@ -132,6 +146,7 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_AutoSpectralRcpp_assign_af_fluor", (DL_FUNC) &_AutoSpectralRcpp_assign_af_fluor, 4},
     {"_AutoSpectralRcpp_fast_kde2d_cpp", (DL_FUNC) &_AutoSpectralRcpp_fast_kde2d_cpp, 6},
     {"_AutoSpectralRcpp_find_local_maxima", (DL_FUNC) &_AutoSpectralRcpp_find_local_maxima, 2},
     {"_AutoSpectralRcpp_optimize_unmix", (DL_FUNC) &_AutoSpectralRcpp_optimize_unmix, 12},

src/RcppExports.cpp

@@ -0,0 +1,99 @@
+#include <RcppArmadillo.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::plugins(openmp)]]
+
+using namespace arma;
+
+// [[Rcpp::export]]
+Rcpp::IntegerVector assign_af_fluor(const arma::mat& raw_data, 
+                                    const arma::mat& spectra, 
+                                    const arma::mat& af_spectra, 
+                                    int n_threads = 4) {
+  
+  // Setup dimensions
+  int n_cells = raw_data.n_rows;
+  int n_channels = raw_data.n_cols;
+  int n_fluors = spectra.n_rows;
+  int n_af = af_spectra.n_rows;
+  
+  // Pseudoinverse: (S*S')^-1 * S
+  mat P = solve(spectra * spectra.t(), spectra);
+  mat S_t = spectra.t();
+  mat AF_t = af_spectra.t(); // channels x af
+  
+  // v.library: (n_fluors x n_af) - how much each AF looks like fluors
+  mat v_library = P * AF_t;
+  
+  // r.library: (n_channels x n_af) - the residual AF
+  mat r_library = AF_t - (S_t * v_library);
+  
+  // Denominator for k calculation: colSums(r.library^2)
+  vec r_dots(n_af);
+  for(int j = 0; j < n_af; ++j) {
+    r_dots[j] = dot(r_library.col(j), r_library.col(j));
+  }
+  
+  // Transpose raw_data to make cell access contiguous (column-major memory)
+  mat Y_t = raw_data.t();
+  
+  // Output: Best AF index per cell
+  Rcpp::IntegerVector best_indices(n_cells);
+  
+#ifdef _OPENMP
+  omp_set_num_threads(n_threads);
+#endif
+  
+  // Parallel Cell Loop
+#pragma omp parallel
+{
+  // Thread-local buffers to avoid repeated allocations
+  vec init_fluor(n_fluors);
+  
+#pragma omp for schedule(static)
+  for (int i = 0; i < n_cells; ++i) {
+    const double* y_ptr = Y_t.colptr(i);
+    
+    // Step A: Initial unmix for this cell (P * y)
+    for(int f = 0; f < n_fluors; ++f) {
+      double sum = 0;
+      for(int c = 0; c < n_channels; ++c) {
+        sum += P(f, c) * y_ptr[c];
+      }
+      init_fluor[f] = sum;
+    }
+    
+    double min_l1_error = datum::inf;
+    int best_af_val = 0;
+    
+    // Step B: Iterate through AF variants
+    for (int j = 0; j < n_af; ++j) {
+      // Calculate k_i: (y . r_j) / dot(r_j, r_j)
+      double numerator = 0;
+      const double* r_ptr = r_library.colptr(j);
+      for(int c = 0; c < n_channels; ++c) {
+        numerator += y_ptr[c] * r_ptr[c];
+      }
+      double ki = numerator / r_dots[j];
+      
+      // Calculate L1 error: sum(abs(init_fluor - ki * v_j))
+      double current_l1 = 0;
+      const double* v_ptr = v_library.colptr(j);
+      for(int f = 0; f < n_fluors; ++f) {
+        current_l1 += std::abs(init_fluor[f] - (ki * v_ptr[f]));
+      }
+      
+      if (current_l1 < min_l1_error) {
+        min_l1_error = current_l1;
+        best_af_val = j + 1; // R-style 1-indexing
+      }
+    }
+    best_indices[i] = best_af_val;
+  }
+}
+
+return best_indices;
+}