vectorize multcomp_pmt

Author: qddyy

R/Studentized.R

@@ -60,30 +60,32 @@ Studentized <- R6Class(
 
         .define = function() {
             private$.statistic_func <- function(data, group) {
+                i <- private$.group_ij$i
+                j <- private$.group_ij$j
+
                 inv_lengths <- 1 / tabulate(group)
 
-                weights <- 1 / sqrt(outer(inv_lengths, inv_lengths, `+`))
+                weights <- 1 / sqrt(inv_lengths[i] + inv_lengths[j])
 
                 if (private$.scoring == "none") {
-                    N <- length(data)
-                    k <- group[N]
+                    weights <- weights * sqrt(
+                        length(data) - length(inv_lengths)
+                    )
 
                     function(data, group) {
                         means <- rowsum.default(data, group) * inv_lengths
 
-                        weights_ <- weights / sqrt(
-                            sum((data - means[group])^2) / (N - k)
+                        (means[i] - means[j]) * weights / sqrt(
+                            sum((data - means[group])^2)
                         )
-
-                        function(i, j) (means[i] - means[j]) * weights_[i, j]
                     }
                 } else {
                     weights <- weights / sd(data)
 
                     function(data, group) {
                         means <- rowsum.default(data, group) * inv_lengths
 
-                        function(i, j) (means[i] - means[j]) * weights[i, j]
+                        (means[i] - means[j]) * weights
                     }
                 }
             }

inst/include/pmt/impl_ksample_pmt.hpp

@@ -35,5 +35,51 @@ RObject impl_ksample_pmt(
         }
     }
 
+    return static_cast<RObject>(statistic_container);
+}
+
+template <bool progress, typename T>
+RObject impl_multcomp_pmt(
+    const NumericVector data,
+    IntegerVector group,
+    T&& statistic_func,
+    const double n_permu)
+{
+    Stat<progress> statistic_container;
+
+    int c = C(*std::prev(group.end()), 2);
+
+    auto statistic_closure = statistic_func(data, group);
+    auto multcomp_update = [&statistic_container, &statistic_closure, data, group, c]() {
+        double* statistic_ptr = statistic_closure(data, group);
+
+        for (int i = 1; i < c; i++) {
+            statistic_container << *statistic_ptr++;
+        }
+
+        return statistic_container << *statistic_ptr;
+    };
+
+    statistic_container.allocate(c, n_permu != 0 ? n_permu : n_permutation(group));
+
+#ifdef SETJMP
+    SETJMP(statistic_func)
+#endif
+
+    multcomp_update();
+
+    if (!std::isnan(n_permu)) {
+        statistic_container.switch_ptr();
+        if (n_permu == 0) {
+            do {
+                multcomp_update();
+            } while (next_permutation(group));
+        } else {
+            do {
+                random_shuffle(group);
+            } while (multcomp_update());
+        }
+    }
+
     return static_cast<RObject>(statistic_container);
 }

inst/include/pmt/impl_multcomp_pmt.hpp

@@ -20,12 +20,12 @@ class CachedFunc : public Function {
     template <typename... Args>
     auto operator()(Args&&... args) const
     {
-        return fast_invoker(Function::operator()(std::forward<Args>(args)...), std::forward<Args>(args)...);
+        return _fast_invoker(Function::operator()(std::forward<Args>(args)...), std::forward<Args>(args)...);
     }
 
-protected:
+private:
     template <typename... Args>
-    auto fast_invoker(SEXP closure, Args&&... args) const
+    auto _fast_invoker(SEXP closure, Args&&... args) const
     {
 #if defined(R_VERSION) && R_VERSION >= R_Version(4, 5, 0)
         Shield<SEXP> closure_formals(R_ClosureFormals(closure)), closure_body(R_ClosureBody(closure)), closure_envir(R_ClosureEnv(closure));
@@ -96,22 +96,15 @@ SEXP ksample_pmt(
         impl_ksample_pmt<false, CachedFunc<double>>(data, clone(group), statistic_func, n_permu);
 }
 
-#include "pmt/impl_multcomp_pmt.hpp"
-
-template <typename T>
 class MultcompFunc : public CachedFunc<SEXP> {
 public:
     using CachedFunc<SEXP>::CachedFunc;
 
     template <typename... Args>
     auto operator()(Args&&... args) const
     {
-        return [statistic_closure = CachedFunc<SEXP>::operator()(std::forward<Args>(args)...), i = RObject(Rf_allocVector(INTSXP, 1)), j = RObject(Rf_allocVector(INTSXP, 1)), this](auto&&...) {
-            return [pairwise_closure = this->fast_invoker(statistic_closure(), i, j), i_ptr = INTEGER(i), j_ptr = INTEGER(j)](int i, int j) {
-                *i_ptr = i;
-                *j_ptr = j;
-                return as<T>(pairwise_closure());
-            };
+        return [statistic_closure = CachedFunc<SEXP>::operator()(std::forward<Args>(args)...)](auto&&... args) {
+            return REAL(statistic_closure(std::forward<decltype(args)>(args)...));
         };
     }
 };
@@ -125,8 +118,8 @@ SEXP multcomp_pmt(
     const bool progress)
 {
     return progress ?
-        impl_multcomp_pmt<true, MultcompFunc<double>>(data, clone(group), statistic_func, n_permu) :
-        impl_multcomp_pmt<false, MultcompFunc<double>>(data, clone(group), statistic_func, n_permu);
+        impl_multcomp_pmt<true, MultcompFunc>(data, clone(group), statistic_func, n_permu) :
+        impl_multcomp_pmt<false, MultcompFunc>(data, clone(group), statistic_func, n_permu);
 }
 
 #include "pmt/impl_paired_pmt.hpp"