(быстрое) попарное сравнение столбцов матрицы, элементы которых имеют формат a/b

Это то, что я мог предоставить до сих пор:

fun1 <- function (S) {
  n <- ncol(S)
  ref2 <- combn(colnames(S), 2)
  ref1 <- paste(ref2[1, ], ref2[2, ], sep = "&")
  z <- matrix(0, choose(n, 2), 3L, dimnames = list(ref1, 0:2))
  k <- 1L
  for (j in 1:(n - 1)) {
    x <- scan(text = S[, j], what = integer(), sep = "/", quiet = TRUE)
    for (i in (j + 1):n) {
      y <- scan(text = S[, i], what = integer(), sep = "/", quiet = TRUE)
      count <- tabulate(.colSums(x == y, 2L, length(x) / 2L) + 1L)
      z[k, ] <- count / sum(count)
      k <- k + 1L
      }
    }
  z
  }

Выглядит плохо, так как имеет гнездо двойного цикла, написанное на R, но самое внутреннее ядро ​​чрезвычайно эффективно благодаря использованию scan, .colSums и tabulate. Общее количество итераций равно choose(ncol(S), 2), что не слишком много для вашей матрицы из 150 столбцов. Я могу заменить fun1 версией Rcpp, если хотите.

## your data
S <- structure(c("0/0", "1/1", "1/2", "0/0", "0/0", "0/1", "1/1", 
"0/3", "0/0", "0/0", "0/0", "0/1", "1/1", "2/2", "0/0", "1/1", 
"0/1", "2/2", "0/0", "0/0"), .Dim = c(5L, 4L), .Dimnames = list(
NULL, c("A", "B", "C", "D")))

fun1(S)
#      0   1   2
#A&B 0.2 0.2 0.6
#A&C 0.2 0.4 0.4
#A&D 0.2 0.4 0.4
#B&C 0.4 0.4 0.2
#B&D 0.2 0.4 0.4
#C&D 0.6 0.0 0.4

Производительность

Ха, когда я на самом деле проверил свою функцию на матрице 15000 x 150, я обнаружил, что:

1. Я мог бы переместить scan из гнезда цикла для ускорения, то есть я мог бы сканировать символьную матрицу в целочисленную матрицу за один раз;
2. scan(text = blabla) занимает вечность, а scan(file = blabla) работает быстро, поэтому, возможно, стоит прочитать данные из текстового файла;
3. работа с текстовым файлом чувствительна к формату файла, поэтому написать надежный код сложно.

Я создал версию fun2 с доступом к файлам и версию fun3 с использованием Rcpp для гнезда цикла. Оказывается, что:

1. чтение из файла действительно лучше (но мы должны предоставить файл в формате, разделенном "/");
2. Реализация Rcpp цикла полезна.

Я вернулся и разместил их здесь (см. редакция 2), и я увидел user20650, начиная с strsplit. Я исключил strsplit из своего варианта, когда начал, потому что я думаю, что работа со строкой может быть медленной. Да, он медленный, но все же быстрее, чем scan. Поэтому я написал fun4, используя strsplit, и соответствующий fun5 с помощью Rcpp (см. редакция 3). Профилирование говорит, что 60% времени выполнения тратится на strsplit, так что это действительно убийца производительности. Затем я заменил strsplit, unlist, as.integer и matrix одной более простой реализацией на C++. Это дает 10-кратный импульс!! Что ж, это разумно, если хорошенько подумать. Используя atoi (или strtol) из библиотеки C <stdlib.h>, мы можем напрямую преобразовывать строки в целые числа, поэтому все операции над строками исключаются!

Короче говоря, я предоставляю только окончательную, самую быструю версию.

library(Rcpp)

cppFunction("IntegerMatrix getInt (CharacterMatrix Char) {
  int m = Char.nrow(), n = Char.ncol();
  IntegerMatrix Int(2 * m, n);
  char *s1, *s2;
  int i, *iptr = &Int(0, 0);
  for (i = 0; i < m * n; i++) {
    s1 = (char *)Char[i]; s2 = s1;
    while(*s2 != '/') s2++; *iptr++ = atoi(s1);
    s2++; *iptr++ = atoi(s2);
    }
  return Int;
  }")

cppFunction('NumericMatrix pairwise(NumericMatrix z, IntegerMatrix Int) {
  int m = Int.nrow() / 2, n = Int.ncol();
  int i, j, k, *x, *y, count[3], *end; bool b1 = 0, b2 = 0;
  double M = 1 / (double)m;
  for (k = 0, j = 0; j < (n - 1); j++) {
    end = &Int(2 * m, j);
    for (i = j + 1; i < n; i++, k++) {
      x = &Int(0, j); y = &Int(0, i);
      count[0] = 0; count[1] = 0; count[2] = 0;
      for (; x < end; x += 2, y += 2) {
        b1 = (x[0] == y[0]);
        b2 = (x[1] == y[1]);
        count[(int)b1 + (int)b2]++;
        }
      z(k, 0) = (double)count[0] * M;
      z(k, 1) = (double)count[1] * M;
      z(k, 2) = (double)count[2] * M;
      }
    }
  return z;
  }')

fun7 <- function (S) {
  ## separate rows using Rcpp; `Int` is an integer matrix
  n <- ncol(S)
  Int <- getInt(S)
  m <- nrow(Int) / 2
  ## initialize the resulting matrix `z`
  ref2 <- combn(colnames(S), 2)
  ref1 <- paste(ref2[1, ], ref2[2, ], sep = "&")
  z <- matrix(0, choose(n, 2), 3L, dimnames = list(ref1, 0:2))
  ## use Rcpp for pairwise summary
  pairwise(z, Int)
  }

Давайте сгенерируем случайную матрицу 15000 x 150 и попробуем.

sim <- function (m, n) {
  matrix(sample(c("0/0", "0/1", "1/0", "1/1"), m * n, TRUE), m, n,
         dimnames = list(NULL, 1:n))
  }

S <- sim(15000, 150)
system.time(oo <- fun7(S))
#   user  system elapsed 
#  1.324   0.000   1.325

О, это молниеносно!

Можно ли исключить значения «0/0» между двумя парами, чтобы получить пропорции? т.е. при сравнении А и В исключить при А=В= 0/0 и получить пропорции для остальных?

Такая адаптация проста на уровне C/C++. Просто дополнительный if тест.

## a new C++ function `pairwise_exclude00`
cppFunction('NumericMatrix pairwise_exclude00(NumericMatrix z, IntegerMatrix Int) {
  int m = Int.nrow() / 2, n = Int.ncol();
  int i, j, k, *x, *y, count[3], size, *end;
  bool b1 = 0, b2 = 0, exclude = 0;
  double M; 
  for (k = 0, j = 0; j < (n - 1); j++) {
    end = &Int(2 * m, j);
    for (i = j + 1; i < n; i++, k++) {
      x = &Int(0, j); y = &Int(0, i);
      count[0] = 0; count[1] = 0; count[2] = 0; size = 0;
      for (; x < end; x += 2, y += 2) {
        b1 = (x[0] == y[0]);
        b2 = (x[1] == y[1]);
        exclude = (x[0] == 0) & (x[1] == 0) & b1 & b2;
        if (!exclude) {
          count[(int)b1 + (int)b2]++;
          size++;
          }
        }
      M = 1 / (double)size;
      z(k, 0) = (double)count[0] * M;
      z(k, 1) = (double)count[1] * M;
      z(k, 2) = (double)count[2] * M;
      }
    }
  return z;
  }')

## re-define `fun7` with a new logical argument `exclude00`
fun7 <- function (S, exclude00) {
  ## separate rows using Rcpp; `Int` is an integer matrix
  n <- ncol(S)
  Int <- getInt(S)
  m <- nrow(Int) / 2
  ## initialize the resulting matrix `z`
  ref2 <- combn(colnames(S), 2)
  ref1 <- paste(ref2[1, ], ref2[2, ], sep = "&")
  z <- matrix(0, choose(n, 2), 3L, dimnames = list(ref1, 0:2))
  ## use Rcpp for pairwise summary
  if (exclude00) pairwise_exclude00(z, Int)
  else pairwise(z, Int)
  }

Используя пример S в вашем вопросе:

fun7(S, TRUE)
#            0         1         2
#A&B 0.3333333 0.3333333 0.3333333
#A&C 0.3333333 0.6666667 0.0000000
#A&D 0.3333333 0.6666667 0.0000000
#B&C 0.5000000 0.5000000 0.0000000
#B&D 0.3333333 0.6666667 0.0000000
#C&D 0.7500000 0.0000000 0.2500000

Здесь используются идеи из ответа 李哲源; особенно tabulate -- дает небольшое ускорение. Для данных 15000x160 требуется ~ 14 секунд на моем старом ноутбуке.

# split strings and form matrix for each column
ap =  matrix(unlist(strsplit(m, "/")), nc=2, byrow=TRUE)
ap = split.data.frame(ap, rep(colnames(m), each=nrow(m))) # maybe a way to use array?

# get 2-way combination of column names
co = combn(colnames(m), 2)

# test equality of each matrix
ap = apply(co, 2, function(x) tabulate(rowSums(ap[[x[1]]]==ap[[x[2]]])+1, 3))

# output
data.frame(t(co), t(ap)/nrow(m))

данные

m = as.matrix(read.table(header=T, text='       A     B       C     D
 "0/0" "0/1"   "0/0" "1/1"
 "1/1" "1/1"   "0/1" "0/1"
 "1/2" "0/3"   "1/1" "2/2"
 "0/0" "0/0"   "2/2" "0/0"
 "0/0" "0/0"   "0/0" "0/0"'))

m = do.call(cbind, replicate(40 , m, simplify = FALSE))
m = do.call(rbind, replicate(3000, m, simplify = FALSE))
colnames(m) =  paste0("A", 1:160)

Вы можете создать 3 функции для указания условий 0,1,2, а затем перебирать имена столбцов, чтобы иметь разные пары, и применять функции для создания результирующего data.frame:

library(tidyr)
matrix <- read.csv("matrix.csv", stringsAsFactors = F)
n <-nrow(matrix)
c <- ncol(matrix)
zero <- function(A, B){ res <- sum(!grepl("0", A) & !grepl("0", B))/n }
one <- function(A, B) {
  A <- unlist(str_split(A, "/"))
  B <- unlist(str_split(B, "/"))
  comp <-data.frame(cbind(A==B, c(1,2), id= sort(rep(1:n,2))))%>%spread(V2, V1)
  res <- sum(sum(comp[,2]+comp[,3])>0)/n} 
two <- function(A, B){res <- sum(A==B)/n}  

res <-data.frame()
k <-1
for (i in 1:(c-1)){
  for (j in (i+1):c){
    A<-matrix[,i]
    B<-matrix[,j]
    res[k,1] <- colnames(matrix)[i]
    res[k,2] <- colnames(matrix)[j]
    res[k,3] <- zero(A,B)
    res[k,4] <- one(A,B)
    res[k,5] <- two(A,B)
    k <-k+1
  }
}
colnames(res) <-c("G1", "G2", "0", "1", "2")