выборка подграфов разных размеров с помощью igraph

Один из способов ускорить ваш код дальше, чем это возможно в базе R, - это использовать пакет Rcpp. Рассмотрим следующую реализацию полной операции Rcpp. В качестве входных данных он принимает следующее:

• valid: индексы всех узлов, которые находятся в достаточно большом компоненте.
• el, deg, firstPos: представление списка ребер графа (соседями узла i являются el[firstPos[i]], el[firstPos[i]+1], ..., el[firstPos[i]+deg[i]-1]).
• size: размер подграфа для выборки
• nrep: количество повторений
• weights: веса ребер хранятся в V(G)$weight
• RWRNodeweight: веса ребер хранятся в V(G)$RWRNodeweight

library(Rcpp)
cppFunction(
"NumericVector scores(IntegerVector valid, IntegerVector el, IntegerVector deg,
                      IntegerVector firstPos, const int size, const int nrep,
                      NumericVector weights, NumericVector RWRNodeweight) {
  const int n = deg.size();
  std::vector<bool> used(n, false);
  std::vector<bool> neigh(n, false);
  std::vector<int> neighList;
  std::vector<double> scores(nrep);
  for (int outerIter=0; outerIter < nrep; ++outerIter) {
    // Initialize variables
    std::fill(used.begin(), used.end(), false);
    std::fill(neigh.begin(), neigh.end(), false);
    neighList.clear();

    // Random first node
    int recent = valid[rand() % valid.size()];
    used[recent] = true;
    double wrSum = weights[recent] * RWRNodeweight[recent];
    double rrSum = RWRNodeweight[recent] * RWRNodeweight[recent];

    // Each additional node
    for (int idx=1; idx < size; ++idx) {
      // Add neighbors of recent
      for (int p=firstPos[recent]; p < firstPos[recent] + deg[recent]; ++p) {
        if (!neigh[el[p]] && !used[el[p]]) {
          neigh[el[p]] = true;
          neighList.push_back(el[p]);
        }
      }

      // Compute new node to add from all neighbors
      int newPos = rand() % neighList.size();
      recent = neighList[newPos];
      used[recent] = true;
      wrSum += weights[recent] * RWRNodeweight[recent];
      rrSum += RWRNodeweight[recent] * RWRNodeweight[recent];

      // Remove from neighList
      neighList[newPos] = neighList[neighList.size() - 1];
      neighList.pop_back();
    }

    // Compute score from wrSum and rrSum
    scores[outerIter] = wrSum / sqrt(rrSum);
  }
  return NumericVector(scores.begin(), scores.end());
}
")

Теперь все, что нам нужно сделать в базе R, - это сгенерировать аргументы для scores, что можно сделать довольно легко:

josilber.rcpp <- function(size, num.rep, G) {
  n <- length(V(G)$name)

  # Determine which nodes fall in sufficiently large connected components
  comp <- components(G)
  valid <- which(comp$csize[comp$membership] >= size)

  # Construct an edge list representation for use in the Rcpp code
  el <- get.edgelist(G, names=FALSE) - 1
  el <- rbind(el, el[,2:1])
  el <- el[order(el[,1]),]
  deg <- degree(G)
  first.pos <- c(0, cumsum(head(deg, -1)))

  # Run the proper number of replications
  scores(valid-1, el[,2], deg, first.pos, size, num.rep,
         as.numeric(V(G)$weight), as.numeric(V(G)$RWRNodeweight))
}

Время на выполнение 1000 репликаций невероятно быстро по сравнению с исходным кодом и всеми igraph решениями, которые мы видели до сих пор (обратите внимание, что для большей части этого тестирования я тестировал исходные функции josilber и random_network для 1 репликации вместо 1000, потому что тестирование для 1000 займет слишком много времени):

• Размер = 10: 0,06 секунды (ускорение в 1200 раз по сравнению с моей ранее предложенной функцией josilber и ускорение в 4000 раз по сравнению с исходной функцией random_network)
• Размер = 100: 0,08 секунды (ускорение на 8700 раз по сравнению с моей ранее предложенной функцией josilber и ускорение в 162000 раз по сравнению с исходной функцией random_network)
• Размер = 1000: 0,13 секунды (ускорение в 32000 раз по сравнению с моей ранее предложенной функцией josilber и ускорение в 20,4 миллиона раз по сравнению с исходной функцией random_network)
• Размер = 5000: 0,32 секунды (увеличение на 68000 раз по сравнению с моей ранее предложенной функцией josilber и ускорение в 290 миллионов раз по сравнению с исходной функцией random_network)

Короче говоря, Rcpp, вероятно, делает возможным вычисление 1000 реплик для каждого размера от 5 до 500 (эта операция, вероятно, займет в общей сложности около 1 минуты), в то время как вычисление 1000 реплик для каждого из этих размеров с использованием чистый код R, который предлагался до сих пор.

По сути, ваш алгоритм выборки графа можно описать как инициализацию набора узлов как случайно выбранный узел с последующим итеративным добавлением соседа вашего текущего набора до тех пор, пока не будет больше соседей или у вас не будет желаемого размера подмножества.

Дорогостоящая повторяющаяся операция здесь - это определение соседей текущего набора, что вы делаете со следующим:

tmp.neigh <- V(G)[unlist(neighborhood(G,1,seed))]$name
tmp.neigh <- setdiff(tmp.neigh, seed)

Короче говоря, вы смотрите на соседей каждого узла в выбранном подмножестве на каждой итерации. Более эффективный подход - хранить вектор соседей и обновлять его каждый раз, когда вы добавляете новый узел; это будет более эффективно, потому что вам нужно учитывать только соседей нового узла.

josilber <- function(size, num.rep, G) {
  score_fun <- function(vert) sum(vert$weight*vert$RWRNodeweight)/sqrt(sum(vert$RWRNodeweight^2))
  n <- length(V(G)$name)

  # Determine which nodes fall in sufficiently large connected components
  comp <- components(G)
  valid <- which(comp$csize[comp$membership] >= size)

  perm <- replicate(num.rep, {
    first.node <- sample(valid, 1)
    used <- (1:n) == first.node  # Is this node selected?
    neigh <- (1:n) %in% neighbors(G, first.node)  # Does each node neighbor our selections?
    for (iter in 2:size) {
      new.node <- sample(which(neigh & !used), 1)
      used[new.node] <- TRUE
      neigh[neighbors(G, new.node)] <- TRUE
    }
    score_fun(V(G)[used])
  })
  perm
}

Для одной реплики это дает значительное ускорение по сравнению с одной репликой кода в вопросе:

• Для size = 50 одна репликация занимает 0,3 секунды для этого кода и 3,8 секунды для опубликованного кода.
• Для size = 100 одна репликация занимает 0,6 секунды для этого кода и 15,2 секунды для опубликованного кода.
• Для size = 200 одна реплика занимает 1,5 секунды для этого кода и 69,4 секунды для опубликованного кода.
• Для size = 500 одна реплика для этого кода занимает 2,7 секунды (таким образом, 1000 реплик должны занять около 45 минут); Я не тестировал ни одной копии опубликованного кода.

Как упоминалось в других ответах, распараллеливание может еще больше повысить производительность выполнения 1000 реплик для заданного размера графа. Далее используется пакет doParallel для распараллеливания повторяющегося шага (настройка практически идентична той, которую сделал @Chris в его ответе):

library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)
josilber2 <- function(size, num.rep, G) {
  score_fun <- function(vert) sum(vert$weight*vert$RWRNodeweight)/sqrt(sum(vert$RWRNodeweight^2))
  n <- length(V(G)$name)

  # Determine which nodes fall in sufficiently large connected components
  comp <- components(G)
  valid <- which(comp$csize[comp$membership] >= size)

  perm <- foreach (i=1:num.rep, .combine='c') %dopar% {
    library(igraph)
    first.node <- sample(valid, 1)
    used <- (1:n) == first.node  # Is this node selected?
    neigh <- (1:n) %in% neighbors(G, first.node)  # Does each node neighbor our selections?
    for (iter in 2:size) {
      new.node <- sample(which(neigh & !used), 1)
      used[new.node] <- TRUE
      neigh[neighbors(G, new.node)] <- TRUE
    }
    score_fun(V(G)[used])
  }
  perm
}

На моем Macbook Air для запуска josilber(100, 1000, my_graph) требуется 670 секунд (это непараллельная версия), а для josilber2(100, 1000, my_graph) требуется 239 секунд (это параллельная версия, сконфигурированная с 4 рабочими процессами). Таким образом, в случае size=100 мы получили 20-кратное ускорение за счет улучшения кода и дополнительное 3-кратное ускорение за счет распараллеливания, то есть общее ускорение в 60 раз.

У меня нет полного ответа, но вот несколько вещей, которые следует учесть, чтобы ускорить его (при условии, что нет более быстрого подхода с использованием другого метода).

1. Удалите из вашего графа любые узлы, которые не являются частью компонента такого размера, как вы ищете. Это действительно будет зависеть от структуры вашей сети, но похоже, что ваши сети являются генами, поэтому, вероятно, существует много генов с очень низкой степенью, и вы можете получить некоторое ускорение, удалив их. Что-то вроде этого кода:

   cgraph <- clusters(G)
   tooSmall <- which(cgraph$csize < size)
   toKeep <- setdiff(1:length(V(G)), which(cgraph$membership %in% tooSmall))
   graph <- induced.subgraph(G, vids=toKeep)
   

2. Подумайте о том, чтобы запустить это параллельно, чтобы использовать несколько ядер. Например, используя пакет parallel и mclapply.

   library(parallel)
   genesets.length<- seq(5, 500)
   names(genesets.length) <- genesets.length
   genesets.length.null.dis <- mclapply(genesets.length, mc.cores=7,
                                        function(length) {
                                          random_network(size=length, G=my_graph)
                                        })
   

Я думаю, что было бы намного эффективнее использовать функцию клик в igraph, поскольку клика - это подграф полностью связанных узлов. Просто установите min и max равными размеру подграфа, который вы ищете, и он вернет все клики размера 5. Затем вы можете выбрать любое подмножество из них, которое соответствует вашим потребностям. К сожалению, с примером графа Эрдоша-Реньи, который вы сгенерировали, часто бывает, что наибольшая клика меньше 5, поэтому для этого примера это не сработает. Тем не менее, он должен работать нормально для реальной сети, которая демонстрирует большую кластеризацию, чем график Эрдоша-Реньи, что, скорее всего, делает ваша.

library(igraph)
##Should be 0.003, not 0.0003 (145000/choose(10000,2))
my_graph <- erdos.renyi.game(10000, 0.003)

cliques(my_graph,min=5,max=5)

У вас есть ряд проблем с вашим кодом (вы не выделяете заранее векторы и т. Д.). Пожалуйста, посмотрите код, который я придумал ниже. Однако я тестировал его только до подграфа размера 100. Однако при увеличении размера подграфа экономия скорости несколько увеличивается по сравнению с вашим кодом. Вам также следует установить пакет foreach. Я запустил это на ноутбуке с 4 ядрами, 2,1 ГГц.

random_network_new <- function(gsize, G) {
  score_fun <- function(g) {
    subsum <- sum(V(g)$weight * V(g)$RWRNodeweight) / sqrt(sum(V(g)$RWRNodeweight^2))
  }

  genes.idx <- V(G)$name

  perm <- foreach (i=seq_len(1e3), .combine='c') %dopar% {
    seed <- rep(0, length=gsize)
    seed[1] <- sample(genes.idx, 1)

    for (j in 2:gsize) {
      tmp.neigh <- neighbors(G, as.numeric(seed[j-1]))
      tmp.neigh <- setdiff(tmp.neigh, seed)
      if (length(tmp.neigh) > 0) {
        seed[j] <- sample(tmp.neigh, 1)
      } else {
        break
      }
    }
    score_fun(induced.subgraph(G, seed))
  }
  perm
}

Обратите внимание, что я переименовал функцию в random_network_new, а аргумент - в gsize.

system.time(genesets <- random_network_new(gsize=100, G=my_graph))                                            
   user   system  elapsed 
1011.157    2.974  360.925 
system.time(genesets <- random_network_new(gsize=50, G=my_graph))
   user  system elapsed 
822.087   3.119 180.358 
system.time(genesets <- random_network_new(gsize=25, G=my_graph))
   user  system elapsed 
379.423   1.130  74.596 
system.time(genesets <- random_network_new(gsize=10, G=my_graph))
   user  system elapsed 
144.458   0.677  26.508 

Один пример с использованием вашего кода (мой более чем в 10 раз быстрее для подграфа 10; это было бы намного быстрее с более крупными подграфами):

system.time(genesets_slow <- random_network(10, my_graph))
   user  system elapsed 
350.112   0.038 350.492