Merge branch 'dev' into dev-cuda

.gitmodules

@@ -1,6 +1,3 @@
 [submodule "bench/nvbench"]
 	path = bench/nvbench
 	url = https://github.com/NVIDIA/nvbench.git
-[submodule "thirdparty/Catch2"]
-	path = thirdparty/Catch2
-	url = https://github.com/catchorg/Catch2.git

CMakeLists.txt

@@ -3,6 +3,9 @@ cmake_policy(SET CMP0048 NEW)
 
 project(PROPR CXX)
 
+set(PROPR_ENABLE_CUDA "AUTO" CACHE STRING "Enable CUDA backend: AUTO, ON, or OFF")
+set_property(CACHE PROPR_ENABLE_CUDA PROPERTY STRINGS AUTO ON OFF)
+set(PROPR_HAS_CUDA OFF)
 
 if ("${CMAKE_BUILD_TYPE}" STREQUAL "" OR ${CMAKE_BUILD_TYPE} STREQUAL "NOMODE")
     message("WORKING ON NO MODE")
@@ -31,31 +34,48 @@ else ()
     message("R              : NOT FOUND")
 endif ()
 
-check_language(CUDA)
-if (CMAKE_CUDA_COMPILER)
-    set(USE_CUDA ON)
-    message("Build CUDA Support")
-    set(PROPERTIES CUDA_ARCHITECTURES 90)
-    enable_language(CUDA)
-    include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
-    set(USE_CUDA ON)
-    set(CMAKE_CUDA_STANDARD 17)
-    set(CMAKE_CUDA_STANDARD_REQUIRED ON)
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --extended-lambda --expt-relaxed-constexpr -lineinfo -O3 --use_fast_math --extra-device-vectorization")
-    message("CUDA           : FOUND")
-    find_package(CUDAToolkit REQUIRED)
-    execute_process(
-        COMMAND nvidia-smi --query-gpu=compute_cap --format=csv
-        COMMAND tail -n 1
-        COMMAND tr -d .
-        OUTPUT_VARIABLE NV_CC
-        OUTPUT_STRIP_TRAILING_WHITESPACE
-    )
-    message("Detected Compute Capability : " ${NV_CC})
-    set(PROPERTIES CUDA_ARCHITECTURES ${NV_CC})
-    set(CMAKE_CUDA_ARCHITECTURES ${NV_CC})
+if (NOT (PROPR_ENABLE_CUDA STREQUAL "AUTO" OR PROPR_ENABLE_CUDA STREQUAL "ON" OR PROPR_ENABLE_CUDA STREQUAL "OFF"))
+    message(FATAL_ERROR "PROPR_ENABLE_CUDA must be one of AUTO, ON, OFF")
+endif ()
+
+if (PROPR_ENABLE_CUDA STREQUAL "ON" OR PROPR_ENABLE_CUDA STREQUAL "AUTO")
+    check_language(CUDA)
+    if (CMAKE_CUDA_COMPILER)
+        set(USE_CUDA ON)
+        set(PROPR_HAS_CUDA ON)
+        message("Build CUDA Support")
+        enable_language(CUDA)
+        include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
+        set(CMAKE_CUDA_STANDARD 17)
+        set(CMAKE_CUDA_STANDARD_REQUIRED ON)
+        set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --extended-lambda --expt-relaxed-constexpr -lineinfo -O3 --extra-device-vectorization") #--use_fast_math 
+        message("CUDA           : FOUND")
+        find_package(CUDAToolkit REQUIRED)
+
+        execute_process(
+            COMMAND nvidia-smi --query-gpu=compute_cap --format=csv,noheader
+            COMMAND head -n 1
+            COMMAND tr -d .
+            OUTPUT_VARIABLE NV_CC
+            OUTPUT_STRIP_TRAILING_WHITESPACE
+            ERROR_QUIET
+            RESULT_VARIABLE NVSIMI_RC
+        )
+
+        if (NVSIMI_RC EQUAL 0 AND NOT NV_CC STREQUAL "")
+            set(CMAKE_CUDA_ARCHITECTURES ${NV_CC})
+            message("Detected Compute Capability : " ${NV_CC})
+        else ()
+            set(CMAKE_CUDA_ARCHITECTURES 75)
+            message("Detected Compute Capability : unavailable, defaulting to 75")
+        endif ()
+    elseif (PROPR_ENABLE_CUDA STREQUAL "ON")
+        message(FATAL_ERROR "CUDA requested (PROPR_ENABLE_CUDA=ON) but CUDA toolkit/compiler not found.")
+    else ()
+        message("CUDA           : NOT FOUND (building CPU-only backend)")
+    endif ()
 else ()
-    message(FATAL_ERROR "CUDA           : NOT FOUND")
+    message("CUDA           : DISABLED (PROPR_ENABLE_CUDA=OFF)")
 endif ()
 
 if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
@@ -73,9 +93,10 @@ option(PROPR_BUILD_TESTS "Option to enable building tests" OFF)
 add_subdirectory(src)
 message("Install Path   : ${CMAKE_INSTALL_PREFIX}")
 message("PROPR Tests    : ${PROPR_BUILD_TESTS}")
+message("PROPR CUDA     : ${PROPR_HAS_CUDA}")
 
 install(DIRECTORY inst/include DESTINATION ${CMAKE_INSTALL_PREFIX}/PROPR
         FILES_MATCHING
             PATTERN "*.hpp"
             PATTERN "*.h"
-            PATTERN "*.cuh")
+            PATTERN "*.cuh")

NAMESPACE

@@ -20,6 +20,7 @@ export(logratio_with_alpha)
 export(logratio_without_alpha)
 export(pcor.bshrink)
 export(propd)
+export(propdGenewise)
 export(propr)
 export(ratios)
 export(results_to_matrix)

OLDCMakeLists.txt

@@ -73,6 +73,14 @@ set(PROPR_LIB_SOURCES
 
 set(PROPR_BENCHMARK_SOURCES bench/src/main_benchmark.cu )
 
+set(PROPR_TEST_SOURCES
+        tests/cpp_test/test_lrm.cpp
+        tests/cpp_test/test_lrv.cpp
+        tests/cpp_test/test_main.cpp
+        tests/cpp_test/test_omega.cpp
+)
+
+
 add_library(propr SHARED ${PROPR_LIB_SOURCES})
 
 target_compile_features(propr PUBLIC cxx_std_14)
@@ -91,8 +99,9 @@ target_include_directories(propr PUBLIC
         ${CUDAToolkit_INCLUDE_DIRS}
         # Local project headers (e.g., src/include/interface, src/include/kernels etc.)
         inst/include
-        # Prereq headers (if they are part of the main library)
+        # Prereq headers (e.g., Catch2 for tests, if they are part of the main library)
         # The original CMakeLists.txt had `include_directories(prereq)`.
+        # Assuming `prereq` includes `prereq/catch2` and other necessary headers.
         ${CMAKE_CURRENT_SOURCE_DIR}/prereq
 )
 
@@ -120,3 +129,16 @@ set_target_properties(propr_benchmark PROPERTIES
         CUDA_STANDARD 14
         CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}"
 )
+
+add_executable(propr_tests ${PROPR_TEST_SOURCES})
+target_link_libraries(propr_tests PRIVATE propr)
+target_include_directories(propr_tests PRIVATE
+        ${CMAKE_CURRENT_SOURCE_DIR}/prereq/catch2
+        inst/include
+)
+
+target_compile_features(propr_tests PUBLIC cxx_std_14)
+set_target_properties(propr_tests PROPERTIES CUDA_STANDARD 14 CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}")
+
+enable_testing()
+add_test(NAME run_propr_cpp_tests COMMAND propr_tests)

R/2-propd-genewise.R

@@ -0,0 +1,112 @@
+#' Convert pairwise propd results into genewise
+#' 
+#' This function converts pairwise propd results into genewise results. The resulting
+#' genewise results are direct indirectors of genes being differentially expressed.
+#' 
+#' @param propd A \code{\link{propd}} object, with FDR values from updateF. Note:
+#' for the moment, only theta results with F-stats are supported. Later we will look 
+#' on the option to get FDR values based on permutations.
+#' @param pairwise_fdr FDR threashold to consider a pairwise relationship as significant.
+#' Default is 0.05.
+#' @param backend Backend to use for connectivity statistics (`"auto"`, `"cpu"`, or `"cuda"`).
+#' @return A data frame with genewise results that can be used to identify differentially
+#' expressed genes. It contains the following columns:
+#'  - "id": gene identifier
+#'  - "lfc": Log Fold Change of the gene, using the geometric mean of all genes as reference.
+#'  - "lrmD": Log Ratio Mean Difference of the gene. Equivalent to the LFC, but using a subset 
+#'     of genes as reference (only the ones that are significantly connected to the gene).
+#'  - "connectivity": number of significant pairwise relationships the gene has.
+#'  - "wconnectivity": weighted connectivity, which sums the strength of significant pairwise 
+#'     relationships (defined as 1/theta) for the gene.
+#' - "FDR_mean": average FDR value across all pairwise comparisons for the gene, which can be 
+#'     used as a genewise significance statistic. Note that this is a simple average and may 
+#'     not be the most robust method for determining genewise significance, but it provides a 
+#'     starting point for identifying genes of interest based on their pairwise relationships. 
+#'     Future updates may include more sophisticated methods for calculating genewise p-values 
+#'     or FDR values based on the pairwise results.
+#' 
+#' @rdname propdGenewise
+#' @export
+propdGenewise <- function(propd, pairwise_fdr = 0.05, backend='auto') {
+
+  # for the moment it only works with theta results with F-stats,
+  # but later we will look into the option to get FDR values based on permutations.
+  if (!"FDR" %in% colnames(propd@results)) {
+    stop("Please run updateF on the propd object to get FDR values before running propdGenewise.")
+  }
+
+  # not working for alpha != NA too
+  if (!is.na(propd@alpha)) {
+    stop("propdGenewise currently only works for alpha = NA. Future updates may include support for other alpha values.")
+  }
+
+  # not working for more than 2 groups too
+  if (length(unique(propd@group)) > 2) {
+    stop("propdGenewise currently only works for 2 groups. Future updates may include support for more than 2 groups.")
+  }
+
+  # get features and number of features
+  features <- colnames(propd@counts)
+  nfeatures <- length(features)
+
+  pair <- propd@results$Pair
+  partner <- propd@results$Partner
+
+  if (!is.numeric(pair) || !is.numeric(partner)) {
+    pair <- match(pair, features)
+    partner <- match(partner, features)
+    if (any(is.na(pair)) || any(is.na(partner))) {
+      stop("Some features in Pair/Partner are not present in propd@counts.")
+    }
+  }
+
+  ## ---- Connectivity + weighted connectivity + FDR mean (CPU/CUDA dispatch) ----
+  stats <- genewiseConnectivityRcpp(
+      partner = as.integer(partner),
+      pair = as.integer(pair),
+      theta = as.numeric(propd@results$theta),
+      fdr = as.numeric(propd@results$FDR),
+      num_genes = nfeatures,
+      pairwise_fdr = pairwise_fdr,
+      backend = backend
+  )
+
+  connectivity <- as.integer(stats$connectivity)
+  wconnectivity <- as.numeric(stats$wconnectivity)
+  fdr_mean <- as.numeric(stats$FDR_mean)
+
+  ## ---- Build FDR adjacency for lrmD ----
+  fdr_mat <- results_to_matrix(propd@results, what = "FDR", features = features)
+  adj <- (fdr_mat > 0) & (fdr_mat < pairwise_fdr)
+
+  ## ---- Build lrm matrices ----
+  lrm1_all <- results_to_matrix(propd@results, what = "lrm1", features = features)
+  lrm2_all <- results_to_matrix(propd@results, what = "lrm2", features = features)
+  # results_to_matrix returns symmetric matrices, but lrm values are directed:
+  # lrm(Partner, Pair) = mean(log(x_Partner / x_Pair)), with Partner > Pair.
+  # Negate the upper triangle so that mat[g, j] = mean(log(x_g / x_j)) for all g, j.
+  lrm1_all[upper.tri(lrm1_all)] <- -lrm1_all[upper.tri(lrm1_all)]
+  lrm2_all[upper.tri(lrm2_all)] <- -lrm2_all[upper.tri(lrm2_all)]
+
+  ## ---- LFC (CLR-based log fold change) ----
+  # lrm differences represent log fold changes; averaging across all genes as
+  # reference is equivalent to using the geometric mean (CLR transformation)
+  lrm_diff <- lrm1_all - lrm2_all
+  lfc <- rowMeans(lrm_diff, na.rm = TRUE) / log(2)
+
+  ## ---- lrmD (LFC using only significant partners as reference) ----
+  lrm_diff <- ifelse(adj, lrm_diff, NA) # keep only significant pairwise relationships
+  lrmD <- apply(lrm_diff, 1, median, na.rm = TRUE) / log(2)
+
+  ## ---- Compile results into a data frame ----
+  data.frame(
+    id = features,
+    lfc = lfc,
+    lrmD = lrmD,
+    connectivity = connectivity,
+    wconnectivity = wconnectivity,
+    FDR_mean = fdr_mean,
+    stringsAsFactors = FALSE,
+    row.names = NULL
+  )
+}