From 9d1811c5dace752f5e703cc9d9de0ebb502506cc Mon Sep 17 00:00:00 2001 From: Josh Buckner Date: Tue, 20 Oct 2026 11:38:44 -0400 Subject: [PATCH] gpuMi, gpuRowTtest, gpuCorPerm, Spearman ranks: the method fixes from widestats gpuMi gave NaN (undefined behaviour, in fact: floor(NaN) cast to int) for a constant variable, built its B-spline weights in the requested precision from narrowed data (a value on a knot, as discrete data put them, could land in the neighbouring bin in single), accepted splineOrder > bins, subtracted entropies, and looped over samples per pair. Now: weights in double from the double data, every pair's joint distribution from one GEMM per tile, the information in log form from each block's own marginals, exactly 0 for a constant variable, never below 0, 1 <= splineOrder <= bins <= 64 checked, defaults Daub et al.'s bins = 10, splineOrder = 3. Help page rewritten: the data need not lie in [0, 1] (the kernel always rescaled), order is not degree. gpuRowTtest with equal groups compared a labelling's swap, computed by different roundings, against the observed t from the first block only; for strongly separated rows the swap could fail the allowance and the p-value came out too small. Labellings are canonical (sample 1 in group 1) and the observed t is recomputed with every block; a within-group sum of squares within 16 eps of its total is zero; a complete separation is +/-Inf with p 0 (Welch's df NA; Student's n - 2) instead of NaN; rows with no test stay out of the maxT maxima. gpuCorPerm's allowance is relative, as gpuRowTtest's. Spearman ranks are taken in double before any narrowing (gpuCor, gpuCorPairs, gpuCorPerm). Reviewed against widestats (weights to 1e-14, Daub's worked example, tiling at several tile sizes, the block loop for every block size); the full suite passes on an A5500. Tests pin the equal-group count and the complete-separation count exactly from the draws, constant variables, knot-valued data in single vs double, and the argument checks. Co-Authored-By: Claude Fable 5.1 --- NEWS.md | 22 +++ R/gpuCorPerm.R | 4 +- R/gpuMi.R | 67 ++++---- R/gpuRowTtest.R | 21 ++- inst/cuda/mi.cu | 245 ++++++++++++----------------- inst/cuda/perm.cu | 5 +- inst/cuda/rowtest.cu | 45 ++++-- man/gpuCorPerm.Rd | 4 +- man/gpuMi.Rd | 114 ++++++++------ man/gpuRowTtest.Rd | 13 +- src/cudaUtils.cpp | 2 +- src/mi.cpp | 260 ++++++++++++------------------- src/mi.h | 10 +- src/rinterface.cpp | 16 +- src/rowtest.cpp | 50 +++--- tests/testthat/test-granger-mi.R | 32 ++++ tests/testthat/test-rowtest.R | 59 +++++++ 17 files changed, 532 insertions(+), 437 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7d73779..b9ec003 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,27 @@ # gputools 2.1.1 +- `gpuMi()` gave `NaN` (or worse) for a constant variable, built its + B-spline weights in single precision when `precision = "single"` (a + value on a knot, as discrete data put them, could land in the + neighbouring bin), accepted `splineOrder > bins`, and computed every + pair's joint distribution in a per-pair loop. It now builds the weights + in double from the double data, takes every pair's joint distribution + from one matrix product per tile, computes the information in log form + (no large entropies subtracted; never below 0), returns exactly 0 for a + constant variable, checks `1 <= splineOrder <= bins`, and defaults to + Daub et al.'s `bins = 10, splineOrder = 3` instead of plain binning into + two bins. Its help page no longer claims the data must lie in [0, 1]. +- `gpuRowTtest()` with equal groups could count a permutation equivalent + to the observed labelling as more extreme, or fail to count its swap, + giving p-values too small for strongly separated rows: labellings are + now canonical and the observed t is recomputed with every block of + permutations. A row constant within each group with the means apart is + `Inf`/`-Inf` with p 0 rather than `NaN`; rows with no test no longer + enter the maxT maxima. +- `gpuCorPerm()`'s rounding allowance is relative, as `gpuRowTtest()`'s. +- Spearman ranks (`gpuCor`, `gpuCorPairs`, `gpuCorPerm`) are taken in + double before any narrowing to single, so distinct values no longer tie. + - `gpuGranger()` reported F = 0, p = 1, "no causality", for a pair it could not resolve (collinear lags, a constant series, a cause identical to the effect), computed in single precision by default on uncentred diff --git a/R/gpuCorPerm.R b/R/gpuCorPerm.R index 95fe441..1adecaa 100644 --- a/R/gpuCorPerm.R +++ b/R/gpuCorPerm.R @@ -41,7 +41,7 @@ gpuCorPerm <- function(x, y = NULL, B = 1000L, method = "pearson", perms <- vapply(seq_len(B), function(b) sample.int(n), integer(n)) - 1L dim(perms) <- c(n, B) # a permuted correlation equal to the observed one counts as at least - # as extreme; this allows for rounding in the comparison + # as extreme; a relative allowance for rounding in the comparison eps <- if(precision == "double") 1e-10 else 1e-5 ans <- .Call(C_gpuCorPerm, x, y, perms, eps, precision) @@ -52,7 +52,7 @@ gpuCorPerm <- function(x, y = NULL, B = 1000L, method = "pearson", # Westfall and Young's single-step maxT: how often the permutation's # largest |r| over all the tests reaches this one maxima <- sort(ans$maxima) - reach <- B - findInterval(abs(r) - eps, maxima, left.open = TRUE) + reach <- B - findInterval(abs(r) * (1 - eps), maxima, left.open = TRUE) adjusted <- (reach + 1) / (B + 1) p[is.na(r)] <- NA adjusted[is.na(r)] <- NA diff --git a/R/gpuMi.R b/R/gpuMi.R index 2e331ed..bf88b83 100644 --- a/R/gpuMi.R +++ b/R/gpuMi.R @@ -1,37 +1,36 @@ -gpuMi <- function(x, y = NULL, bins = 2, splineOrder = 1, - precision = getOption("gputools.precision", "single")) +gpuMi <- function(x, y = NULL, bins = 10, splineOrder = 3, + precision = getOption("gputools.precision", "single")) { - precision <- .precision(precision) - x <- as.matrix(x) - if(!is.null(y)) { - y <- as.matrix(y) - } + precision <- .precision(precision) + x <- as.matrix(x) + storage.mode(x) <- "double" + if (!is.null(y)) { + y <- as.matrix(y) + storage.mode(y) <- "double" + if (nrow(y) != nrow(x)) stop("x and y have different numbers of rows") + } + whole <- function(v) length(v) == 1 && is.numeric(v) && !is.na(v) && + v == round(v) + if (!whole(bins) || !whole(splineOrder) || splineOrder < 1 || + splineOrder > bins) { + stop("bins and splineOrder must be whole numbers with ", + "1 <= splineOrder <= bins") + } + if (bins > 64) stop("bins must be at most 64") + if (nrow(x) < 2) stop("need at least 2 rows") + if (anyNA(x) || (!is.null(y) && anyNA(y))) { + stop("x and y must have no missing values") + } + bins <- as.integer(bins) + splineOrder <- as.integer(splineOrder) + nsamples <- as.integer(nrow(x)) + na <- as.integer(ncol(x)) + nb <- if (is.null(y)) na else as.integer(ncol(y)) - bins <- as.integer(bins) - splineOrder <- as.integer(splineOrder) - - nsamples <- as.integer(nrow(x)) - na <- as.integer(ncol(x)) - a <- as.double(x) - - b <- a - nb <- na - row_labels <- colnames(x) - - if(is.null(y)) { - nb <- as.integer(ncol(x)) - b <- as.double(x) - row_labels <- colnames(x) - } else { - nb <- as.integer(ncol(y)) - b <- as.double(y) - row_labels <- colnames(y) - } - - mi <- .Call(C_gpuMutualInfo, bins, splineOrder, nsamples, a, na, b, nb, - precision) - mi <- matrix(mi, nb, na) - rownames(mi) <- row_labels - colnames(mi) <- colnames(x) - return(mi) + mi <- .Call(C_gpuMutualInfo, bins, splineOrder, nsamples, x, na, y, nb, + precision) + mi <- matrix(mi, nb, na) + rownames(mi) <- if (is.null(y)) colnames(x) else colnames(y) + colnames(mi) <- colnames(x) + mi } diff --git a/R/gpuRowTtest.R b/R/gpuRowTtest.R index d856c72..36ab194 100644 --- a/R/gpuRowTtest.R +++ b/R/gpuRowTtest.R @@ -40,10 +40,17 @@ gpuRowTtest <- function(x, group, var.equal = FALSE, B = 0L, se <- sqrt(a + b) df <- (a + b)^2 / (a^2 / (n1 - 1) + b^2 / (n2 - 1)) } + # a zero standard error with the means apart is a complete separation + # of the groups: t infinite, p 0 (as SciPy reports it); with the means + # equal as well, or a constant row, there is no test t <- (m1 - m2) / se - t[se == 0] <- NaN + t[se == 0 & m1 == m2] <- NaN + p.value <- 2 * pt(-abs(t), df) + # (Welch's df is 0/0 there, so the p-value is set directly) + p.value[is.infinite(t)] <- 0 + if(!var.equal) df[is.infinite(t)] <- NA res <- data.frame(mean1 = m1, mean2 = m2, statistic = t, df = df, - p.value = 2 * pt(-abs(t), df), row.names = rownames(x)) + p.value = p.value, row.names = rownames(x)) if(B > 0) { # the permutations, drawn here so set.seed() reproduces them, on @@ -51,6 +58,16 @@ gpuRowTtest <- function(x, group, var.equal = FALSE, B = 0L, labels <- matrix(0L, n, B + 1L) labels[, 1] <- in1 for(b in seq_len(B)) labels[, b + 1L] <- in1[sample.int(n)] + if(n1 == n2) { + # With equal groups a labelling and its swap give the same |t| + # in exact arithmetic but are computed by different roundings, + # which for strongly separated rows can differ by more than the + # allowance: a swap of the observed labelling could then fail + # to count, making p-values too small. Putting sample 1 in + # group 1 always makes the two the same column, bit for bit. + flip <- labels[1, ] == 0L + labels[, flip] <- 1L - labels[, flip] + } eps <- if(precision == "double") 1e-10 else 1e-5 ans <- .Call(C_gpuRowTtest, x, labels, n1, isTRUE(var.equal), eps, precision) diff --git a/inst/cuda/mi.cu b/inst/cuda/mi.cu index cb1a5dd..08c597a 100644 --- a/inst/cuda/mi.cu +++ b/inst/cuda/mi.cu @@ -1,157 +1,118 @@ -#define NTHREADS 16 - -__global__ void scale(real knot_max, int nx, int nsamples, - real * x, int pitch_x) +// B-spline mutual information (Daub et al. 2004) of every pair of +// variables, as src/mi.cpp drives it: +// +// - each variable's values get B-spline weights over `bins` bins +// (miWeights, always in double from the double data: rounding the data +// to single first moves a value on a knot, as discrete data put them, +// into the neighbouring bin), laid out as W (n x nvars*bins, +// column-major: sample k of variable v's bin i at W[k + (v*bins + i)*n]) +// so that one matrix product W_x' W_y / n (cuBLAS, from the host) holds +// the joint distribution of every pair as a bins x bins block; +// - miPairs reads each block and computes the mutual information in log +// form, sum P log2(P / (P_a P_b)) with the marginals P_a, P_b the +// block's own row and column sums: the same as H_a + H_b - H_ab, but +// without subtracting large entropies, and consistent with any rounding +// in the product. 0 log 0 counts as 0. +// +// bins + order - 1 is at most MI_MAXSLOTS (the R code limits bins to 64). + +#define MI_MAXSLOTS 128 + +// the knot vector 0 (order times), 1, ..., bins - order, span (order +// times), span = bins - order + 1; knot i of bins + order +__device__ double miKnot(int i, int bins, int order) { - int - col_idx = blockDim.x * blockIdx.x + threadIdx.x; - - if(col_idx >= nx) return; - - real - min, max, - * col = x + col_idx * pitch_x; - - // find the min and the max - min = max = col[0]; - for(int i = 1; i < nsamples; i++) { - if(col[i] < min) min = col[i]; - if(col[i] > max) max = col[i]; - } - - real delta = max - min; - for(int i = 0; i < nsamples; i++) - col[i] = (knot_max * (col[i] - min)) / delta; -} - -__device__ real do_fraction(real numer, real denom) { - real result = real(0); - - if((numer == denom) && (numer != real(0))) - result = real(1); - else if(denom != real(0)) - result = numer / denom; - - return result; + if (i < order) return 0.0; + if (i <= bins - 1) return (double) (i - order + 1); + return (double) (bins - order + 1); } -// bins must be initialized to zero before calling get_bin_scores -__global__ void get_bin_scores(int nbins, int order, - int nknots, real * knots, int nsamples, - int nx, real * x, int pitch_x, - real * bins, int pitch_bins) +// Each variable's weights: one thread per variable, over its n samples. +// A constant variable is put wholly at the domain's start (every value in +// the first bin). +__global__ void miWeights(int n, int nvars, int bins, int order, + const double * x, double * w) { - int - col_x = blockDim.x * blockIdx.x + threadIdx.x; - - if(col_x >= nx) - return; - - real - ld, rd, z, - term1, term2, - * in_col = x + col_x * pitch_x, - * bin_col = bins + col_x * pitch_bins; - int i0; - - for(int k = 0; k < nsamples; k++, bin_col += nbins) { - z = in_col[k]; - i0 = (int)floor(z) + order - 1; - if(i0 >= nbins) - i0 = nbins - 1; - - bin_col[i0] = real(1); - for(int i = 2; i <= order; i++) { - for(int j = i0 - i + 1; j <= i0; j++) { - rd = do_fraction(knots[j + i] - z, knots[j + i] - knots[j + 1]); - - if((j < 0) || (j >= nbins) || (j >= nknots) || (j + i - 1 < 0) || (j > nknots)) - term1 = real(0); - else { - ld = do_fraction(z - knots[j], - knots[j + i - 1] - knots[j]); - term1 = ld * bin_col[j]; - } - - if((j + 1 < 0) || (j + 1 >= nbins) || (j + 1 >= nknots) || (j + i < 0) || (j + i >= nknots)) - term2 = real(0); - else { - rd = do_fraction(knots[j + i] - z, - knots[j + i] - knots[j + 1]); - term2 = rd * bin_col[j + 1]; - } - bin_col[j] = term1 + term2; + int v = blockIdx.x * blockDim.x + threadIdx.x; + if (v >= nvars) return; + const double * col = x + (size_t) v * n; + double lo = col[0], hi = col[0]; + for (int k = 1; k < n; k++) { + if (col[k] < lo) lo = col[k]; + if (col[k] > hi) hi = col[k]; + } + double width = hi > lo ? hi - lo : 1.0; + int span = bins - order + 1, slots = bins + order - 1; + double wk[MI_MAXSLOTS]; + for (int k = 0; k < n; k++) { + double z = (col[k] - lo) / width * span; + // order 1: the knot interval holding the value, among the non-empty + // ones order - 1 .. bins - 1; the far end of the domain belongs to the + // last + int cell = (int) floor(z); + if (cell < 0) cell = 0; + if (cell > span - 1) cell = span - 1; + cell += order - 1; + for (int i = 0; i < slots; i++) wk[i] = 0.0; + wk[cell] = 1.0; + // Cox-de Boor: B_{i,k} = (z - t_i)/(t_{i+k-1} - t_i) B_{i,k-1} + // + (t_{i+k} - z)/(t_{i+k} - t_{i+1}) B_{i+1,k-1}, 0/0 as 0 + for (int kk = 2; kk <= order; kk++) { + int count = bins + order - kk; + for (int i = 0; i < count; i++) { + double ti = miKnot(i, bins, order), tik = miKnot(i + kk, bins, order); + double lden = miKnot(i + kk - 1, bins, order) - ti; + double rden = tik - miKnot(i + 1, bins, order); + double left = lden > 0 ? (z - ti) / lden * wk[i] : 0.0; + double right = rden > 0 ? (tik - z) / rden * wk[i + 1] : 0.0; + wk[i] = left + right; } } + for (int i = 0; i < bins; i++) + w[k + ((size_t) v * bins + i) * n] = wk[i]; } } -__global__ void get_entropy(int nbins, int nsamples, int nx, - real * bin_scores, int pitch_bin_scores, real * entropies) +// Single-precision copy of the weights, for the product in single. +__global__ void miNarrow(const double * in, float * out, size_t count) { - int - col_x = blockDim.x * blockIdx.x + threadIdx.x; - - if(col_x >= nx) - return; - - real - * in_col = bin_scores + col_x * pitch_bin_scores, - entropy = real(0), prob, logp; - - for(int i = 0; i < nbins; i++) { - prob = real(0); - for(int j = 0; j < nsamples; j++) - prob += in_col[j * nbins + i]; - prob /= (double) nsamples; - - if(prob <= real(0)) - logp = real(0); - else - logp = log2(prob); - - entropy += prob * logp; - } - entropies[col_x] = -entropy; + size_t e = blockIdx.x * (size_t) blockDim.x + threadIdx.x; + if (e < count) out[e] = (float) in[e]; } -__global__ void get_mi(int nbins, int nsamples, - int nx, real * x_bin_scores, int pitch_x_bin_scores, - real * entropies_x, - int ny, real * y_bin_scores, int pitch_y_bin_scores, - real * entropies_y, - real * mis, int pitch_mis) +// The mutual information of each pair of a tile from the product j +// (tx*bins x ty*bins, leading dimension ldj) of their weights, scaled by +// invn = 1/n to the joint distribution. Writes mi[(j0 + b) + (i0 + a)*ldm] +// for x variable a and y variable b of the tile, in bits. +__global__ void miPairs(int bins, int tx, int ty, const real * j, int ldj, + real invn, int i0, int j0, real * mi, int ldm) { - int - col_x = blockDim.x * blockIdx.x + threadIdx.x, - col_y = blockDim.y * blockIdx.y + threadIdx.y; - - if((col_x >= nx) || (col_y >= ny)) - return; - - real - prob, logp, mi = real(0), - * x_bins = x_bin_scores + col_x * pitch_x_bin_scores, - * y_bins = y_bin_scores + col_y * pitch_y_bin_scores; - - // calculate joint entropy - for(int i = 0; i < nbins; i++) { - for(int j = 0; j < nbins; j++) { - prob = real(0); - for(int k = 0; k < nsamples; k++) - prob += x_bins[k * nbins + i] * y_bins[k * nbins + j]; - prob /= (real)nsamples; - - if(prob <= real(0)) - logp = real(0); - else - logp = log2(prob); - - mi += prob * logp; + size_t e = blockIdx.x * (size_t) blockDim.x + threadIdx.x; + if (e >= (size_t) tx * ty) return; + int a = e % tx, b = e / tx; + const real * block = j + (size_t) a * bins + (size_t) b * bins * ldj; + real rows[MI_MAXSLOTS], cols[MI_MAXSLOTS]; + for (int i = 0; i < bins; i++) { rows[i] = real(0); cols[i] = real(0); } + for (int c = 0; c < bins; c++) + for (int r = 0; r < bins; r++) { + real p = block[r + (size_t) c * ldj] * invn; + rows[r] += p; + cols[c] += p; } - } - - // calculate mi from entropies - mi += entropies_x[col_x] + entropies_y[col_y]; - (mis + col_y * pitch_mis)[col_x] = mi; + // in logs, each clamped below at the smallest normal number: where an + // entry is positive its row and column sums are at least as large, so + // every log is finite + const real tiny = sizeof(real) == 4 ? real(1.17549435e-38) + : real(2.2250738585072014e-308); + real sum = real(0); + for (int c = 0; c < bins; c++) + for (int r = 0; r < bins; r++) { + real p = block[r + (size_t) c * ldj] * invn; + if (p > real(0)) { + real lr = rows[r] > tiny ? rows[r] : tiny; + real lc = cols[c] > tiny ? cols[c] : tiny; + sum += p * (log2(p) - log2(lr) - log2(lc)); + } + } + mi[(size_t) (j0 + b) + (size_t) (i0 + a) * ldm] = sum; } diff --git a/inst/cuda/perm.cu b/inst/cuda/perm.cu index c61cdf9..0dabddb 100644 --- a/inst/cuda/perm.cu +++ b/inst/cuda/perm.cu @@ -13,7 +13,8 @@ __global__ void gpuPermuteRows(const real * z, int n, int q, const int * perm, } // One permutation's correlations rb (p x q) against the observed r0: -// counts[e] += 1 where |rb| >= |r0| - eps. Entries that are no test (NaN, +// counts[e] += 1 where |rb| >= |r0| (1 - eps), a relative allowance for +// rounding in the comparison. Entries that are no test (NaN, // or the diagonal of x with itself) are set to 0 in rb, so that its // largest |value| is the permutation's maximum statistic. __global__ void gpuCountExceed(real * rb, const real * r0, int p, int q, @@ -26,5 +27,5 @@ __global__ void gpuCountExceed(real * rb, const real * r0, int p, int q, rb[e] = real(0); return; } - if(fabs(v) >= fabs(o) - eps) counts[e]++; + if(fabs(v) >= fabs(o) * (real(1) - eps)) counts[e]++; } diff --git a/inst/cuda/rowtest.cu b/inst/cuda/rowtest.cu index 8fbca76..67ccbe4 100644 --- a/inst/cuda/rowtest.cu +++ b/inst/cuda/rowtest.cu @@ -6,7 +6,12 @@ // sums of squares do not cancel. // t (genes x nb): Welch's t, or Student's with the pooled variance if -// pooled is set; NaN where a variance is zero +// pooled is set. A within-group sum of squares within rounding of zero +// (16 eps of the sum of squares it is taken from) is zero: a row constant +// within the group would otherwise get a large finite t where the true t +// is infinite. A zero standard error with the means apart gives +/-Inf +// (a labelling that separates the groups completely); with the means +// equal, or a constant row (NaN after scaling), there is no test: NaN. __global__ void gpuRowT(const real * s1, const real * q1, const real * s, const real * q, int genes, int nb, int n1, int n2, int pooled, real * t) @@ -16,7 +21,11 @@ __global__ void gpuRowT(const real * s1, const real * q1, const real * s, int g = e % genes; real a1 = s1[e], b1 = q1[e], a2 = s[g] - a1, b2 = q[g] - b1; real m1 = a1 / n1, m2 = a2 / n2; - real ss1 = fmax(b1 - a1 * m1, real(0)), ss2 = fmax(b2 - a2 * m2, real(0)); + const real tiny = real(16) * (sizeof(real) == 4 ? real(1.1920929e-7) + : real(2.220446049250313e-16)); + real ss1 = b1 - a1 * m1, ss2 = b2 - a2 * m2; + if(ss1 <= tiny * b1) ss1 = real(0); + if(ss2 <= tiny * b2) ss2 = real(0); real se2; if(pooled) { real v = (ss1 + ss2) / (n1 + n2 - 2); @@ -24,26 +33,32 @@ __global__ void gpuRowT(const real * s1, const real * q1, const real * s, } else { se2 = ss1 / (n1 - 1) / n1 + ss2 / (n2 - 1) / n2; } - t[e] = se2 > real(0) ? (m1 - m2) / sqrt(se2) : REAL_NAN; + real diff = m1 - m2; + if(se2 > real(0)) t[e] = diff / sqrt(se2); + else if(diff > real(0)) t[e] = REAL_INF; + else if(diff < real(0)) t[e] = -REAL_INF; + else t[e] = REAL_NAN; } -// For a block of c labellings, global columns first .. first + c - 1 of -// which column 0 is the observed one (skipped): counts[g] += 1 where -// |t_b| >= |t0| (1 - eps), and maxima[b - 1] = the largest |t_b| over the -// genes (NaN counting as 0). One block of 256 threads per column. -__global__ void gpuRowTCount(const real * t, const real * t0, int genes, - int first, real eps, int * counts, double * maxima) +// For a block of labellings whose column 0 is the observed one, recomputed +// with the block (so a labelling equivalent to it is compared with the +// same arithmetic), and columns 1 .. c the permutations first .. first + +// c - 1: counts[g] += 1 where |t_b| >= |t0| (1 - eps), and maxima[first + +// b - 1] = the largest |t_b| over the rows that have a test (a row with +// no observed test, NaN, counts as 0 everywhere). One block of 256 +// threads per permutation column. +__global__ void gpuRowTCount(const real * t, int genes, int first, real eps, + int * counts, double * maxima) { __shared__ double best[256]; - int b = first + blockIdx.x; - if(b == 0) return; - const real * tb = t + (size_t) blockIdx.x * genes; + const real * t0 = t; + const real * tb = t + (size_t) (blockIdx.x + 1) * genes; double m = 0; for(int g = threadIdx.x; g < genes; g += blockDim.x) { real v = tb[g], o = t0[g]; - if(isnan(v)) continue; + if(isnan(o) || isnan(v)) continue; m = fmax(m, (double) fabs(v)); - if(!isnan(o) && fabs(v) >= fabs(o) * (real(1) - eps)) + if(fabs(v) >= fabs(o) * (real(1) - eps)) atomicAdd(counts + g, 1); } best[threadIdx.x] = m; @@ -54,5 +69,5 @@ __global__ void gpuRowTCount(const real * t, const real * t0, int genes, best[threadIdx.x + stride]); __syncthreads(); } - if(threadIdx.x == 0) maxima[b - 1] = best[0]; + if(threadIdx.x == 0) maxima[first + blockIdx.x] = best[0]; } diff --git a/man/gpuCorPerm.Rd b/man/gpuCorPerm.Rd index 962b2ee..62d05b3 100644 --- a/man/gpuCorPerm.Rd +++ b/man/gpuCorPerm.Rd @@ -45,7 +45,9 @@ gpuCorPerm(x, y = NULL, B = 1000L, method = "pearson", never zero. The adjusted p-value counts instead the permutations whose largest \eqn{|r_b|} over all the tests reaches \eqn{|r|}: single-step maxT (Westfall and Young 1993), which controls the family-wise error - rate. + rate. A permuted \eqn{|r_b|} within a relative allowance of \eqn{|r|} + (1e-10 in double, 1e-5 in single) counts as reaching it, so that a + permutation equivalent to the observed data always counts. The whole \code{ncol(x)} by \code{ncol(y)} matrix is held on the GPU, twice; for wider data, split the columns of \code{y} (the maxT diff --git a/man/gpuMi.Rd b/man/gpuMi.Rd index 1a4ed7e..dd8c560 100644 --- a/man/gpuMi.Rd +++ b/man/gpuMi.Rd @@ -1,68 +1,88 @@ \name{gpuMi} \alias{gpuMi} -\title{B spline based mutual information} +\title{B-Spline Mutual Information of Every Pair of Variables on a GPU} \description{ - This function estimates the mutual information for permutations - of pairs of columns of a matrix using a B spline approach on a - GPU device. Please note, the data must be values from the interval - [0.0, 1.0]. + Estimates the mutual information between every column of \code{x} and + every column of \code{y} (of \code{x} itself when \code{y} is missing) + by the B-spline estimator of Daub et al. (2004), on the GPU. Unlike a + correlation it detects any kind of dependence, not only linear or + monotone association. } \usage{ -gpuMi(x, y = NULL, bins = 2, splineOrder = 1, - precision = getOption("gputools.precision", "single")) +gpuMi(x, y = NULL, bins = 10, splineOrder = 3, + precision = getOption("gputools.precision", "single")) } \arguments{ - \item{x}{a matrix of floating point numbers from the interval [0.0, 1.0]. - Each column represents a list of samples of a random variable. The - mutual information between each column of x and each column of y will - be computed. If y is NULL then each pair of columns of x will be - compared.} - \item{y}{a matrix of floating point numbers from the interval [0.0, 1.0]. - Each column represents a list of samples of a random variable. The - mutual information between each column of x and each column of y will - be computed. If y is NULL then each pair of columns of x will be - compared.} - \item{bins}{a single integer value representing the number of equal - intervals that [0.0, 1.0] will be divided into in order to determine - the bins in which to place each value of the columns of x and y. In - the case of splineOrder = 1, this determines the histogram for - traditional mutual information. For splineOrder > 1, a single value - may be placed in multiple adjoining bins with varying weights on - membership.} - \item{splineOrder}{a single integer value giving the degree of the spline - polynomials used to define both the number of bins a single value will - be placed in and the weight of membership given to the value.} - \item{precision}{\code{"single"} or \code{"double"}: the floating point precision - the GPU computes in. Single precision is the default, set with - \code{options(gputools.precision = "double")}; it runs many times - faster on most consumer GPUs. Results are returned as doubles either way.} + \item{x}{a numeric matrix, one variable per column, a sample per row, + with no missing values. The values may lie in any range: each variable + is scaled onto the spline domain by its own minimum and maximum.} + \item{y}{an optional numeric matrix of further variables, the same + samples in the same rows.} + \item{bins}{the number of bins each variable's range is divided into + (\eqn{M} in Daub et al., who find 10 to work well), at most 64.} + \item{splineOrder}{the B-spline order (\eqn{k} in Daub et al.), from 1 to + \code{bins}: 1 puts each value wholly in one bin (plain binning), 2 + spreads it linearly over two, 3 (which Daub et al. recommend) + quadratically over three, and so on. The polynomial degree is + \code{splineOrder - 1}.} + \item{precision}{\code{"single"} or \code{"double"}: the floating point + precision of the products on the GPU. The B-spline weights themselves + are always computed in double from the data as given, since rounding + a value on a knot (as discrete data put them) to single would move it + into the neighbouring bin. Results are returned as doubles either way.} } -\value{ - a matrix of single precision floating point values of order ncol(y) by - ncol(x). Entry $(i, j)$ of this matrix represents the mutual information - calculation for $(y_i, x_j)$. +\details{ + Each variable is scaled linearly onto the spline domain \eqn{[0, M - k + + 1]}, its minimum at 0 and its maximum at the far end, with the knot + vector \eqn{0, \ldots, 0, 1, 2, \ldots, M - k, M - k + 1, \ldots, M - k + + 1} (each end repeated \eqn{k} times), as in Daub et al.; each value + then has B-spline weights over the \eqn{M} bins, non-negative and + summing to 1. With \eqn{w_i(s, a)} the weight of sample \eqn{s} of + variable \eqn{i} in bin \eqn{a}, the bin probabilities are \eqn{P_i(a) = + \mathrm{mean}_s\, w_i(s, a)}{P_i(a) = mean_s w_i(s, a)} and the joint + ones \eqn{P_{ij}(a, b) = \mathrm{mean}_s\, w_i(s, a) w_j(s, b)}{P_ij(a, + b) = mean_s w_i(s, a) w_j(s, b)}, so that \eqn{MI = H_i + H_j - H_{ij}} + with \eqn{H = -\sum P \log_2 P}{H = -sum P log2 P}. + + On the GPU, the weights of many variables side by side make one matrix + product hold the joint distribution of every pair, a \eqn{M \times M}{M + x M} block per pair. Since each sample's weights sum to 1, a block's row + and column sums are the two marginals, and the information is computed + from the block alone as \eqn{\sum P_{ij} \log_2 (P_{ij} / (P_i + P_j))}{sum P_ij log2(P_ij / (P_i P_j))}: the same quantity, without + subtracting large entropies to get a small difference, and consistent + with any rounding in the product. \eqn{0 \log 0}{0 log 0} counts as 0, + and an estimate that rounds just below 0 is 0. +} + +\value{a numeric matrix with a row per column of \code{y} (of \code{x} + when \code{y} is missing) and a column per column of \code{x}: entry + \eqn{(j, i)} is the mutual information of \code{y[, j]} with + \code{x[, i]}, in bits. A constant variable shares no information: 0. + With \code{y} missing, the diagonal is each variable's information with + itself: its entropy for \code{splineOrder = 1}, and less for higher + orders, whose spread weights blur a variable's joint distribution with + itself. } \references{ - Carten O. Daub, Ralf Steuer, Joachim Selbig, and Sebastian Kloska. - 2004. Estimating mutual information using B-spline functions -- an - improved similarity measure for analysing gene expression data. - \emph{BMC Bioinformatics}. 5:118. Available from - \doi{10.1186/1471-2105-5-118} + Daub, C. O., Steuer, R., Selbig, J. and Kloska, S. (2004). Estimating + mutual information using B-spline functions -- an improved similarity + measure for analysing gene expression data. \emph{BMC Bioinformatics} + 5, 118. \doi{10.1186/1471-2105-5-118} } \examples{ \dontshow{if (gpuAvailable()) withAutoprint(\{ # examplesIf} - # get 3 random variables each with 20 samples - x <- matrix(runif(60), 20, 3) - y <- matrix(runif(60), 20, 3) - # do something interesting - y[,2] <- 3.0 * (x[,1] + x[,3]) - z <- gpuMi(x, y, bins = 10, splineOrder = 3) - print(z) +set.seed(1) +u <- runif(2000, -1, 1) +x <- cbind(u = u, usq = u^2, noise = runif(2000, -1, 1)) +mi <- gpuMi(x) +mi["usq", "u"] > 10 * mi["noise", "u"] # u and u^2: dependent, uncorrelated +abs(cor(u, u^2)) < 0.05 \dontshow{\}) # examplesIf} } diff --git a/man/gpuRowTtest.Rd b/man/gpuRowTtest.Rd index 2376f89..4d00ecf 100644 --- a/man/gpuRowTtest.Rd +++ b/man/gpuRowTtest.Rd @@ -50,8 +50,17 @@ gpuRowTtest(x, group, var.equal = FALSE, B = 0L, keeps the dependence between genes, which the adjustment exploits, and makes no assumption of normality. - A constant row, or one constant within both groups, has no test: - its statistic is \code{NaN} and its p-values \code{NA}. + A permuted \eqn{|t_b|} within a relative allowance of \eqn{|t|} + (1e-10 in double, 1e-5 in single) counts as reaching it, so that a + relabelling equivalent to the observed one always counts; with equal + groups a labelling and its swap are made the same labelling, as they + are the same test. A row constant within each group but with the + groups' means apart separates them completely: its statistic is + \code{Inf} or \code{-Inf} and its p-values 0 (as SciPy reports it), + and in the permutations every relabelling that separates the groups + as completely counts. A constant row, or one constant within both + groups with equal means, has no test: its statistic is \code{NaN} and + its p-values \code{NA}. } \value{ diff --git a/src/cudaUtils.cpp b/src/cudaUtils.cpp index ffe5d62..a0df9cf 100644 --- a/src/cudaUtils.cpp +++ b/src/cudaUtils.cpp @@ -57,7 +57,7 @@ const std::vector & kernelFiles() { "rowtest", { "gpuRowT", "gpuRowTCount" } }, { "perm", { "gpuPermuteRows", "gpuCountExceed" } }, { "kmeans", { "gpuAssignNearest", "gpuOneHot" } }, - { "mi", { "scale", "get_bin_scores", "get_entropy", "get_mi" } }, + { "mi", { "miWeights", "miNarrow", "miPairs" } }, { "qrdecomp", { "getColNorms" } }, }; diff --git a/src/mi.cpp b/src/mi.cpp index 4482fc2..36ed5b5 100644 --- a/src/mi.cpp +++ b/src/mi.cpp @@ -1,178 +1,114 @@ -#include -#include -#include -#include -#include -#include -#include +// B-spline mutual information of every pair of variables: Daub et al. +// (2004), with the weights computed in double from the double data and +// the joint distributions of every pair from one matrix product per tile +// (see inst/cuda/mi.cu for the kernels and the formulas). T is the +// precision of the products; the method and its guards follow the Python +// package widestats (mi.py). -#include "cudaUtils.h" +#include +#include +#include -#define NTHREADS 16 +#include "cuda_runtime.h" -template -int initKnots(int nbins, int order, T ** knots) { - int - om1 = order - 1, - degree = nbins - 1, - dpo = degree + order, - nknots = dpo + 1; +#include "R.h" - *knots = R_Calloc(nknots, T); - for(int i = 0; i < nknots; i++) { - if(i <= om1) - (*knots)[i] = T(0); - else if(i <= degree) - (*knots)[i] = (*knots)[i-1] + T(1); - else - (*knots)[i] = (*knots)[degree] + T(1); - } - return nknots; -} +#include "cudaUtils.h" +#include "gpublas.h" +#include "mi.h" +namespace { +// the weights of nvars variables, n x nvars*bins in T on the device template -void bSplineMutualInfo(int nbins, int order, int nsamples, - int nx, const T * x, - int ny, const T * y, - T * out_mi) +T * weights(int n, int nvars, int bins, int order, const double * x) { - size_t - pitch[2], pitch_bins[2], - col_bytes = (size_t)nsamples * sizeof(T); - int - nknots, nblocks[2], size[2] = { nx, ny }; - T - * knots, * dknots, - * stage[2], * dx[2], * dentropy[2], * dbins[2]; - const T - * data[2] = { x, y }; - - nknots = initKnots(nbins, order, &knots); - T knot_max = knots[nknots - 1]; - cudaMalloc((void **)&dknots, nknots * sizeof(T)); - cudaMemcpy(dknots, knots, nknots * sizeof(T), cudaMemcpyHostToDevice); - R_Free(knots); - - checkCudaError("bSplineMutualInfoSingle: 1"); - - for(int i = 0; i < 2; i++) { - cudaMallocPitch((void **)&(dx[i]), pitch + i, col_bytes, size[i]); - cudaMallocHost((void **)&(stage[i]), size[i] * col_bytes); - cudaMalloc((void **)&(dentropy[i]), size[i] * sizeof(T)); - cudaMallocPitch((void **)&(dbins[i]), pitch_bins + i, - nbins * col_bytes, size[i]); - cudaMemset2D(dbins[i], pitch_bins[i], 0, nbins * col_bytes, size[i]); - - nblocks[i] = size[i] / NTHREADS; - if(nblocks[i] * NTHREADS < size[i]) - nblocks[i]++; + double * dX = gpuAlloc((size_t) n * nvars); + CUDART_CALL(cudaMemcpy(dX, x, (size_t) n * nvars * sizeof(double), + cudaMemcpyHostToDevice)); + size_t count = (size_t) n * nvars * bins; + double * dW = gpuAlloc(count); + int threads = 128; + { + int nn = n, nv = nvars, nb = bins, no = order; + void * args[] = { &nn, &nv, &nb, &no, &dX, &dW }; + cudaLaunch("miWeights", args, dim3((nvars + threads - 1) / threads), + dim3(threads), Precision::Double); } - - checkCudaError("bSplineMutualInfoSingle: 2"); - - cudaStream_t stream[2]; - for(int i = 0; i < 2; i++) - cudaStreamCreate(stream + i); - for(int i = 0; i < 2; i++) { - cudaMemcpyAsync(stage[i], data[i], size[i] * col_bytes, - cudaMemcpyHostToHost, stream[i]); - cudaMemcpy2DAsync(dx[i], pitch[i], stage[i], col_bytes, col_bytes, - size[i], cudaMemcpyHostToDevice, stream[i]); + cudaFree(dX); + if (sizeof(T) == sizeof(double)) return (T *) dW; + // the product in single: narrow the weights once + T * dWs = gpuAlloc(count); + { + size_t cnt = count; + void * args[] = { &dW, &dWs, &cnt }; + cudaLaunch("miNarrow", args, + dim3((unsigned) ((count + 255) / 256)), dim3(256), + Precision::Double); } - for(int i = 0; i < 2; i++) { - dim3 - grid(nblocks[i]), - block(NTHREADS); - - // pitches go to the kernels in elements, as the int they take - int xpitch = pitch[i] / sizeof(T); - - void * scaleArgs[] = { - &knot_max, - size + i, - &nsamples, - dx + i, - &xpitch - }; - cudaLaunch("scale", scaleArgs, - grid, block, precisionOf(), stream[i]); + cudaFree(dW); + return dWs; +} - int pitch_bins_i = pitch_bins[i] / sizeof(T); - void * gbsArgs[] = { - &nbins, &order, - &nknots, &dknots, - &nsamples, size + i, - dx + i, &xpitch, - dbins + i, - &pitch_bins_i - }; - cudaLaunch("get_bin_scores", gbsArgs, - grid, block, precisionOf(), stream[i]); +} // namespace - void * entropyArgs[] = { - &nbins, - &nsamples, - size + i, - dbins + i, - &pitch_bins_i, - dentropy + i - }; - cudaLaunch("get_entropy", entropyArgs, - grid, block, precisionOf(), stream[i]); +template +void bSplineMutualInfo(int bins, int order, int n, int nx, const double * x, + int ny, const double * y, double * out) +{ + const bool self = y == NULL || y == x; + if (self) { y = x; ny = nx; } + T * dWx = weights(n, nx, bins, order, x); + T * dWy = self ? dWx : weights(n, ny, bins, order, y); + T * dMi = gpuAlloc((size_t) nx * ny); + // a tile's product holds (tile*bins)^2 entries: about 2^24 of them + int tile = std::max(1, (int) (4096 / bins)); + T * dJ = gpuAlloc((size_t) tile * bins * tile * bins); + T invn = (T) (1.0 / n); + for (int i0 = 0; i0 < nx; i0 += tile) { + int tx = std::min(tile, nx - i0); + for (int j0 = 0; j0 < ny; j0 += tile) { + int ty = std::min(tile, ny - j0); + int ldj = tx * bins; + // every pair's joint distribution (times n) in this tile: a bins x + // bins block of W_x' W_y per pair + blas::gemm('T', 'N', tx * bins, ty * bins, n, T(1), + dWx + (size_t) i0 * bins * n, n, + dWy + (size_t) j0 * bins * n, n, T(0), dJ, ldj); + int nb = bins, ttx = tx, tty = ty, ii0 = i0, jj0 = j0, ldm = ny; + size_t pairs = (size_t) tx * ty; + void * args[] = { &nb, &ttx, &tty, &dJ, &ldj, &invn, &ii0, &jj0, &dMi, + &ldm }; + cudaLaunch("miPairs", args, dim3((unsigned) ((pairs + 255) / 256)), + dim3(256), precisionOf()); + } } - checkCudaError("bSplineMutualInfoSingle: 3"); - - cudaFree(dknots); - for(int i = 0; i < 2; i++) { - cudaFreeHost(stage[i]); - cudaFree(dx[i]); + std::vector h((size_t) nx * ny); + CUDART_CALL(cudaMemcpy(h.data(), dMi, h.size() * sizeof(T), + cudaMemcpyDeviceToHost)); + cudaFree(dJ); cudaFree(dMi); + if (!self) cudaFree(dWy); + cudaFree(dWx); + // out[j + i*ny] is the information of y_j with x_i. A constant variable + // shares no information: exactly 0, not rounding; an estimate of 0 can + // round just below it + std::vector cx(nx, 0), cy(ny, 0); + for (int i = 0; i < nx; i++) { + const double * c = x + (size_t) i * n; + cx[i] = std::equal(c + 1, c + n, c); } - - size_t pitch_mi; - T * dmi; - cudaMallocPitch((void **)&dmi, &pitch_mi, ny * sizeof(T), nx); - - dim3 - gridDim(nblocks[0], nblocks[1]), blockDim(NTHREADS, NTHREADS); - - int - elem_pitch_bins[2] = { (int) (pitch_bins[0] / sizeof(T)), - (int) (pitch_bins[1] / sizeof(T)) }, - dpitch_mi = pitch_mi / sizeof(T); - - void * miArgs[] = { - &nbins, &nsamples, - &nx, dbins, elem_pitch_bins, dentropy, - &ny, dbins + 1, elem_pitch_bins + 1, dentropy + 1, - &dmi, &dpitch_mi - }; - cudaLaunch("get_mi", miArgs, gridDim, blockDim, precisionOf()); - checkCudaError("bSplineMutualInfoSingle: 4"); - - for(int i = 0; i < 2; i++) - cudaFree(dbins[i]); - - T * mi_stage; - cudaMallocHost((void **)&mi_stage, nx * ny * sizeof(T)); - - cudaMemcpy2D(mi_stage, ny * sizeof(T), dmi, pitch_mi, - ny * sizeof(T), nx, cudaMemcpyDeviceToHost); - checkCudaError("bSplineMutualInfoSingle: 5"); - cudaFree(dmi); - - memcpy(out_mi, mi_stage, nx * ny * sizeof(T)); - cudaFreeHost(mi_stage); - checkCudaError("bSplineMutualInfoSingle: 6"); + for (int j = 0; j < ny; j++) { + const double * c = y + (size_t) j * n; + cy[j] = self ? cx[j] : std::equal(c + 1, c + n, c); + } + for (int i = 0; i < nx; i++) + for (int j = 0; j < ny; j++) { + size_t e = (size_t) j + (size_t) i * ny; + out[e] = (cx[i] || cy[j]) ? 0.0 : std::max(0.0, (double) h[e]); + } } -template int initKnots(int nbins, int order, float ** knots); -template int initKnots(int nbins, int order, double ** knots); -template void bSplineMutualInfo(int nbins, int order, int nsamples, - int nx, const float * x, - int ny, const float * y, - float * out_mi); -template void bSplineMutualInfo(int nbins, int order, int nsamples, - int nx, const double * x, - int ny, const double * y, - double * out_mi); +template void bSplineMutualInfo(int, int, int, int, const double *, + int, const double *, double *); +template void bSplineMutualInfo(int, int, int, int, const double *, + int, const double *, double *); diff --git a/src/mi.h b/src/mi.h index 6e3094a..6521041 100644 --- a/src/mi.h +++ b/src/mi.h @@ -1,10 +1,12 @@ #ifndef _MI_H_ #define _MI_H_ +// B-spline mutual information (bits) of every column of x (n x nx) with +// every column of y (n x ny; x itself when y is NULL), B-spline weights +// over `bins` bins of order `order`, the products in precision T. +// out[j + i*ny] is the information of y's column j with x's column i. template -void bSplineMutualInfo(int nbins, int order, int nsamples, - int nx, const T * x, - int ny, const T * y, - T * out_mi); +void bSplineMutualInfo(int bins, int order, int n, int nx, const double * x, + int ny, const double * y, double * out); #endif /* _MI_H_ */ diff --git a/src/rinterface.cpp b/src/rinterface.cpp index c0ce799..f12f90b 100644 --- a/src/rinterface.cpp +++ b/src/rinterface.cpp @@ -539,13 +539,19 @@ SEXP grangerD(SEXP x, int rows, int nx, SEXP y, int ny, int lag) return out; } +// the weights are always computed in double; T is the precision of the +// products template SEXP mutualInfoT(int bins, int order, int nsamples, SEXP x, int na, SEXP y, int nb) { - std::vector a = narrow(x, "x"), b = narrow(y, "y"), mi(na * nb); - bSplineMutualInfo(bins, order, nsamples, na, a.data(), nb, b.data(), - mi.data()); + std::vector a = narrow(x, "x"), b; + const double * py = NULL; + if (!isNull(y)) { b = narrow(y, "y"); py = b.data(); } + else nb = na; + std::vector mi((size_t) na * nb); + bSplineMutualInfo(bins, order, nsamples, na, a.data(), nb, py, + mi.data()); return doubles(mi.data(), mi.size()); } @@ -797,7 +803,9 @@ SEXP gpuRowTtest(SEXP x, SEXP labels, SEXP n1, SEXP pooled, SEXP eps, SEXP gpuRank(SEXP x, SEXP precision) { requireDevice(); - return isDouble(precision) ? rankT(x) : rankT(x); + // ranks are always taken in double: narrowed to single first, distinct + // values that round to the same float would tie + return rankT(x); } SEXP gpuSetDevice(SEXP device) diff --git a/src/rowtest.cpp b/src/rowtest.cpp index 6e066f5..6826b62 100644 --- a/src/rowtest.cpp +++ b/src/rowtest.cpp @@ -25,14 +25,19 @@ void rowTtests(int genes, int n, const T * x, int nb, const int * labels, } for (size_t e = 0; e < g.size(); ++e) g[e] = (T) labels[e]; - // labellings a block at a time, so that genes x block fits in ~128 MB + // labellings a block at a time, so that genes x block fits in ~128 MB; + // each block's product also carries the observed labelling as its + // column 0, recomputed with the block (the same arithmetic, even where + // the library picks kernels by width), so that a permutation equivalent + // to the observed labelling always counts as reaching it int block = (int) std::max(1, std::min(nb, ((size_t) 1 << 25) / genes)); T * dX = gpuAlloc(gn), * dX2 = gpuAlloc(gn), * dS = gpuAlloc(genes), * dQ = gpuAlloc(genes), * dG = gpuAlloc((size_t) n * nb), - * dS1 = gpuAlloc((size_t) genes * block), - * dQ1 = gpuAlloc((size_t) genes * block), - * dT = gpuAlloc((size_t) genes * block), * dT0 = gpuAlloc(genes); + * dGb = gpuAlloc((size_t) n * (block + 1)), + * dS1 = gpuAlloc((size_t) genes * (block + 1)), + * dQ1 = gpuAlloc((size_t) genes * (block + 1)), + * dT = gpuAlloc((size_t) genes * (block + 1)); int * dCounts = gpuAlloc(genes); double * dMax = gpuAlloc(nb > 1 ? nb - 1 : 1); CUDART_CALL(cudaMemcpy(dX, x, gn * sizeof(T), cudaMemcpyHostToDevice)); @@ -48,31 +53,38 @@ void rowTtests(int genes, int n, const T * x, int nb, const int * labels, int n2 = n - n1, pool = pooled; T teps = (T) eps; - for (int b0 = 0; b0 < nb; b0 += block) { + std::vector ht(genes); + // the observed labelling (column 0 of labels) heads every block + CUDART_CALL(cudaMemcpy(dGb, dG, n * sizeof(T), cudaMemcpyDeviceToDevice)); + for (int b0 = 1; b0 < std::max(nb, 2); b0 += block) { + // permutations b0 .. b0 + c - 1 (none when nb == 1: the observed only) int c = std::min(block, nb - b0); - T * dGc = dG + (size_t) b0 * n; - blas::gemm('N', 'N', genes, c, n, T(1), dX, genes, dGc, n, T(0), dS1, + if (c > 0) + CUDART_CALL(cudaMemcpy(dGb + n, dG + (size_t) b0 * n, + (size_t) c * n * sizeof(T), + cudaMemcpyDeviceToDevice)); + int cols = 1 + std::max(c, 0); + blas::gemm('N', 'N', genes, cols, n, T(1), dX, genes, dGb, n, T(0), dS1, genes); - blas::gemm('N', 'N', genes, c, n, T(1), dX2, genes, dGc, n, T(0), dQ1, + blas::gemm('N', 'N', genes, cols, n, T(1), dX2, genes, dGb, n, T(0), dQ1, genes); - size_t entries = (size_t) genes * c; + size_t entries = (size_t) genes * cols; dim3 grid((unsigned) ((entries + 255) / 256)), threads(256); - void * tArgs[] = { &dS1, &dQ1, &dS, &dQ, &genes, &c, &n1, &n2, &pool, + void * tArgs[] = { &dS1, &dQ1, &dS, &dQ, &genes, &cols, &n1, &n2, &pool, &dT }; cudaLaunch("gpuRowT", tArgs, grid, threads, precisionOf()); - if (b0 == 0) - CUDART_CALL(cudaMemcpy(dT0, dT, genes * sizeof(T), - cudaMemcpyDeviceToDevice)); - if (nb > 1) { + if (b0 == 1) + CUDART_CALL(cudaMemcpy(ht.data(), dT, genes * sizeof(T), + cudaMemcpyDeviceToHost)); + if (c > 0) { + int first = b0 - 1; dim3 cgrid(c); - void * cArgs[] = { &dT, &dT0, &genes, &b0, &teps, &dCounts, &dMax }; + void * cArgs[] = { &dT, &genes, &first, &teps, &dCounts, &dMax }; cudaLaunch("gpuRowTCount", cArgs, cgrid, threads, precisionOf()); } + if (nb <= 1) break; } - std::vector ht(genes); - CUDART_CALL(cudaMemcpy(ht.data(), dT0, genes * sizeof(T), - cudaMemcpyDeviceToHost)); for (int r = 0; r < genes; ++r) t0[r] = ht[r]; CUDART_CALL(cudaMemcpy(counts, dCounts, genes * sizeof(int), cudaMemcpyDeviceToHost)); @@ -81,7 +93,7 @@ void rowTtests(int genes, int n, const T * x, int nb, const int * labels, cudaMemcpyDeviceToHost)); cudaFree(dX); cudaFree(dX2); cudaFree(dS); cudaFree(dQ); cudaFree(dG); - cudaFree(dS1); cudaFree(dQ1); cudaFree(dT); cudaFree(dT0); + cudaFree(dGb); cudaFree(dS1); cudaFree(dQ1); cudaFree(dT); cudaFree(dCounts); cudaFree(dMax); } diff --git a/tests/testthat/test-granger-mi.R b/tests/testthat/test-granger-mi.R index 099f213..e1c7d4a 100644 --- a/tests/testthat/test-granger-mi.R +++ b/tests/testthat/test-granger-mi.R @@ -213,3 +213,35 @@ test_that("the reference reproduces the paper's own numbers", { expect_equal(entropyBits(colMeans(w)), log2(10) - 0.6 * log2(3) - 0.4 * log2(4)) }) + +test_that("gpuMi gives a constant variable exactly 0, and checks its arguments", { + skip_if_no_gpu() + set.seed(9) + X <- matrix(rnorm(100 * 3), 100, 3) + X[, 2] <- 3 + m <- gpuMi(X) + expect_false(anyNA(m)) + expect_true(all(m[2, ] == 0) && all(m[, 2] == 0)) + expect_gt(m[1, 1], 0) + expect_true(all(m >= 0)) + expect_equal(gpuMi(X, X), m, ignore_attr = TRUE) + expect_error(gpuMi(X, bins = 3, splineOrder = 4), "splineOrder <= bins") + expect_error(gpuMi(X, bins = 2.5), "whole numbers") + expect_error(gpuMi(X, bins = 65), "at most 64") + expect_error(gpuMi(X, X[1:50, ]), "different numbers of rows") + expect_error(gpuMi(X[1, , drop = FALSE]), "at least 2 rows") +}) + +test_that("gpuMi in single precision keeps knot-valued data in their bins", { + skip_if_no_gpu() + # discrete data sit exactly on the knots; weights rounded to single + # once moved values into the neighbouring bin, so single and double + # disagreed by far more than single's rounding + set.seed(10) + X <- matrix(sample(0:9, 300 * 4, replace = TRUE), 300, 4) + for (b in c(5, 10, 20)) { + expect_equal(gpuMi(X, bins = b, splineOrder = 3, precision = "single"), + gpuMi(X, bins = b, splineOrder = 3, precision = "double"), + tolerance = 1e-5, label = paste("bins", b)) + } +}) diff --git a/tests/testthat/test-rowtest.R b/tests/testthat/test-rowtest.R index a253e15..52a8b25 100644 --- a/tests/testthat/test-rowtest.R +++ b/tests/testthat/test-rowtest.R @@ -89,3 +89,62 @@ test_that("gpuRowTtest names, constant rows and refusals", { x[1, 1] <- NA expect_error(gpuRowTtest(x, group), "missing") }) + +test_that("gpuRowTtest counts a labelling equivalent to the observed one exactly", { + skip_if_no_gpu() + # equal groups: a labelling and its swap are the same test. For a row + # that separates the groups completely, the permutations equal to the + # observed labelling or its swap are exactly the ones that reach it; + # the count of those is known from the draws themselves + set.seed(84) + n <- 10 + x <- matrix(rnorm(50 * n), 50, n) + group <- rep(1:2, each = 5) + x[1, group == 1] <- x[1, group == 1] + 1000 + for (pr in precisions) { + B <- 3000 + set.seed(85) + res <- gpuRowTtest(x, group, B = B, precision = pr) + set.seed(85) + in1 <- group == 1 + same <- 0 + for (b in seq_len(B)) { + lab <- in1[sample.int(n)] + if (all(lab == in1) || all(lab == !in1)) same <- same + 1 + } + expect_equal(res$p.perm[1], (same + 1) / (B + 1), info = pr) + expect_equal(res$p.adjusted[1], (same + 1) / (B + 1), info = pr) + } +}) + +test_that("gpuRowTtest reports a complete separation as infinite", { + skip_if_no_gpu() + set.seed(86) + x <- matrix(rnorm(4 * 8), 4, 8) + group <- rep(c("a", "b"), each = 4) + x[2, ] <- rep(c(1, 2), each = 4) # constant within each group, apart + x[3, ] <- 7 # constant: no test + B <- 500 + set.seed(87) + res <- gpuRowTtest(x, group, B = B) + expect_equal(res$statistic[2], -Inf) + expect_equal(res$p.value[2], 0) + expect_true(is.na(res$df[2])) + # only a relabelling that separates the groups as completely (the + # observed one or its swap) reaches an infinite t: count those draws + set.seed(87) + in1 <- group == "a" + same <- 0 + for (b in seq_len(B)) { + lab <- in1[sample.int(8)] + if (all(lab == in1) || all(lab == !in1)) same <- same + 1 + } + expect_equal(res$p.perm[2], (same + 1) / (B + 1)) + expect_equal(res$p.adjusted[2], (same + 1) / (B + 1)) + expect_true(is.nan(res$statistic[3]) && is.na(res$p.perm[3])) + expect_true(all(is.finite(res$statistic[c(1, 4)]))) + # Student's df is defined there + pooled <- gpuRowTtest(x, group, var.equal = TRUE) + expect_equal(pooled$statistic[2], -Inf) + expect_equal(pooled$df[2], 6) +})