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) +})