Выборка подграфов разных размеров с использованием igraph

Ответ 1

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

• valid: индексы всех узлов, которые находятся в достаточно большом компоненте
• el, deg, firstPos: представление списка ребер графа (node 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 секунды (ускорение 1200x по сравнению с ранее предложенной функцией josilber и ускорение 4000x по сравнению с исходной функцией random_network)
• Размер = 100: 0,08 секунды (ускорение 8700x по сравнению с ранее предложенной функцией josilber и ускорение 162000x по сравнению с исходной функцией random_network)
• Размер = 1000: 0,13 секунды (ускорение 32000x по сравнению с ранее предложенной функцией josilber и 20,4 миллиона раз быстрее по сравнению с исходной функцией random_network)
• Размер = 5000: 0,32 секунды (ускорение 68000x по сравнению с ранее предложенной функцией josilber и 290 миллионов раз быстрее по сравнению с исходной функцией random_network)

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

Ответ 2

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

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

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

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

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
}

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

• Для размера = 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-кратное ускорение от улучшения кода и дополнительное ускорение 3x от распараллеливания для полного ускорения 60x.

Ответ 3

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

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

  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)
  

• Рассмотрите возможность параллельного выполнения этого параллельного использования нескольких ядер. Например, используя пакет 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)
                                       })
  

Ответ 4

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

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)

Ответ 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