R - 查找按位二进制邻居（一次翻转一位）

Posted 2021-03-29

在使用大型矩阵时，是否有更有效的方法来匹配矩阵行？我有一个向量，其值对应于2 ^ N行的矩阵。 N通常很大，例如> 20。每行是N = {0,1}值的唯一组合，表示决策空间上的“位置”。即，对于N = 3，行将是0 0 0,0 0 1,0 1 0,1 0 0，...，1 1 1

我需要确定位置是否是局部最大值，即N个相邻位置是否具有较低值。例如，相应地，对于位置0 0 0，相邻位置是1 0 0,0 1 0和0 0 1。我编写了以下解决方案来完成这项工作，但对于大N来说非常慢。

library(prodlim) #for row.match command

set.seed(1234)
N=10

space = as.matrix(expand.grid(rep(list(0:1), N))) #creates all combinations of 0-1 along N-dimensions

performance = replicate(2^N, runif(1, min=0, max=1)) #corresponding values for each space-row (position)

#determine whether a space position is a local maxima, that is, the N neighboring positions are smaller in performance value


system.time({
local_peaks_pos = matrix(NA,nrow=2^N, ncol=1)
for(v in 1:2^N)
{

  for(q in 1:N)
  {
    temp_local_pos = space[v,1:N]
    temp_local_pos[q] = abs(temp_local_pos[q]-1)

    if(performance[row.match(temp_local_pos[1:N], space[,1:N])] > performance[v])
    {
      local_peaks_pos[v,1] = 0
      break
    }

  }

}
local_peaks_pos[is.na(local_peaks_pos)] = 1
})

  user  system elapsed 
   9.94    0.05   10.06 

答案

正如Gabe在他的评论中提到的，你可以利用你的决策空间可以被解释为单个整数的事实：

set.seed(1234L)
N <- 10L
performance <- runif(2^N)
powers_of_two <- as.integer(rev(2L ^ (0L:(N - 1L))))

is_local_max <- sapply(0L:(2^N - 1), function(i) {
  multipliers <- as.integer(rev(intToBits(i)[1L:N])) * -1L
  multipliers[multipliers == 0L] <- 1L
  neighbors <- i + powers_of_two * multipliers
  # compensate that R vectors are 1-indexed
  !any(performance[neighbors + 1L] > performance[i + 1L])
})

# compensate again
local_peaks_int <- which(is_local_max) - 1L
local_peaks_binary <- t(sapply(local_peaks_int, function(int) {
  as.integer(rev(intToBits(int)[1L:N]))
}))

> head(local_peaks_binary)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    0    0    1    0     0
[2,]    0    0    0    0    1    0    0    1    1     0
[3,]    0    0    0    0    1    1    1    1    0     0
[4,]    0    0    0    1    0    0    0    1    1     1
[5,]    0    0    0    1    0    1    0    1    0     1
[6,]    0    0    0    1    1    0    1    1    1     0

在十进制中，multipliers包含powers_of_two的符号，因此，当添加到当前整数时，它表示二进制位翻转。例如，如果原始二进制文件是0 0并且我们翻转一位来获得1 0，就好像我们在十进制中添加了2 ^ 1，但如果它原来是1 0并且我们翻转一位来获得0 0，那么我们减去2 ^ 1十进制。

local_peaks_binary中的每一行都是来自决策空间的二进制，其中最低有效位在右侧。因此，例如，第一个局部峰值是小数4。

有关整数到二进制的映射，请参阅this question。

编辑：如果你想并行执行：

library(doParallel)
set.seed(1234L)
N <- 20L
performance <- runif(2^N)
powers_of_two <- as.integer(rev(2 ^ (0:(N - 1))))

num_cores <- detectCores()
workers <- makeCluster(num_cores)
registerDoParallel(workers)

chunks <- splitIndices(length(performance), num_cores)
chunks <- lapply(chunks, "-", 1L)
local_peaks_int <- foreach(chunk=chunks, .combine=c, .multicombine=TRUE) %dopar% {
  is_local_max <- sapply(chunk, function(i) {
    multipliers <- as.integer(rev(intToBits(i)[1L:N])) * -1L
    multipliers[multipliers == 0L] <- 1L
    neighbors <- i + powers_of_two * multipliers
    # compensate that R vectors are 1-indexed
    !any(performance[neighbors + 1L] > performance[i + 1L])
  })

  # return
  chunk[is_local_max]
}

local_peaks_binary <- t(sapply(local_peaks_int, function(int) {
  as.integer(rev(intToBits(int)[1L:N]))
}))

stopCluster(workers); registerDoSEQ()

在我的系统中，上述操作在~2.5秒内完成，该系统有4个CPU内核。

---

这是一个使用多线程的C ++版本，但至少在我的4线程系统中， 它看起来并不比Gabe的Fortran版本快 。但是，当我尝试在新会话中运行Gabe的Fortran代码时，我得到N <- 29L：cannot allocate vector of size 4.0 Gb的以下错误。

编辑：显然我改变了一些重要的东西，因为经过再次测试，C ++版本实际上看起来更快。

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel)]]
#include <cstddef> // size_t
#include <vector>
#include <Rcpp.h>
#include <RcppParallel.h>

using namespace std;
using namespace Rcpp;
using namespace RcppParallel;

class PeakFinder : public Worker
{
public:
  PeakFinder(const NumericVector& performance, vector<int>& peaks, const int N)
    : performance_(performance)
    , peaks_(peaks)
    , N_(N)
  { }

  void operator()(size_t begin, size_t end) {
    vector<int> peaks;
    for (size_t i = begin; i < end; i++) {
      bool is_local_peak = true;
      unsigned int mask = 1;
      for (int exponent = 0; exponent < N_; exponent++) {
        unsigned int neighbor = static_cast<unsigned int>(i) ^ mask; // bitwise XOR
        if (performance_[i] < performance_[neighbor]) {
          is_local_peak = false;
          break;
        }

        mask <<= 1;
      }

      if (is_local_peak)
        peaks.push_back(static_cast<int>(i));
    }

    mutex_.lock();
    peaks_.insert(peaks_.end(), peaks.begin(), peaks.end());
    mutex_.unlock();
  }

private:
  const RVector<double> performance_;
  vector<int>& peaks_;
  const int N_;
  tthread::mutex mutex_;
};

// [[Rcpp::export]]
IntegerVector local_peaks(const NumericVector& performance, const int N) {
    vector<int> peaks;
    PeakFinder peak_finder(performance, peaks, N);
    // each thread call will check at least 1024 values
    parallelFor(0, performance.length(), peak_finder, 1024);

    IntegerVector result(peaks.size());
    int i = 0;
    for (int peak : peaks) {
        result[i++] = peak;
    }
    return result;
}

在local-peaks.cpp中保存C ++代码后，此代码：

library(Rcpp)
library(RcppParallel)

sourceCpp("local-peaks.cpp")

set.seed(1234L)
N <- 29L
performance <- runif(2^N)
system.time({
    local_peaks_int <- local_peaks(performance, N)
})

在~2秒内完成（不考虑分配performance所需的时间）。

如果你确实需要二进制表示，你可以像这样更改local_peaks（参见this question）：

// [[Rcpp::export]]
IntegerMatrix local_peaks(const NumericVector& performance, const int N) {
  vector<int> peaks;
  PeakFinder peak_finder(performance, peaks, N);
  // each thread call will check at least 1024 values
  parallelFor(0, performance.length(), peak_finder, 1024);

  // in case you want the same order every time, #include <algorithm> and uncomment next line
  // sort(peaks.begin(), peaks.end());

  IntegerMatrix result(peaks.size(), N);
  int i = 0;
  for (int peak : peaks) {
    for (int j = 0; j < N; j++) {
      result(i, N - j - 1) = peak & 1;
      peak >>= 1;
    }

    i++;
  }

  return result;
}

另一答案

这是一个遵循与示例代码相同的一般结构的解决方案。 intToBits和packBits映射到每个整数的二进制表示形式（减去一个从零开始）。内循环翻转每个N位以获得邻居。在我的笔记本电脑上，N=10只需几分之一秒，N=20只需一分钟。注释代码存储来自已经测试的邻居的一些信息，以便不重做计算。取消注释这些线使它在N=20大约35秒内运行。

loc_max <- rep(1, 2^N)
for (v in 1:2^N){
  ## if (loc_max[v] == 0) next
  vbits <- intToBits(v-1)
  for (q in 1:N){
    tmp <- vbits
    tmp[q] <- !vbits[q]
    pos <- packBits(tmp, type = "integer") + 1
    if (performance[pos] > performance[v]){
      loc_max[v] <- 0
      break
    ## } else {
    ##   loc_max[pos] <- 0
    }
  }
}

identical(loc_max, local_peaks_pos[, 1])
## [1] TRUE

编辑：听起来你需要尽可能快的速度，所以这里有另一个建议，依赖于编译的代码运行速度明显快于我的第一个例子。 N=20只有几分之一秒，N=29只有不到20秒（我可以放入笔记本电脑内存中的最大例子）。

这是使用单核;您既可以将其并行化，也可以在单个核心中运行它，并将您的蒙特卡罗模拟并行化。

library(inline)

loopcode <-
"  integer v, q, pos
   do v = 0, (2**N)-1
      do q = 0, N-1
         if ( btest(v,q) ) then
            pos = ibclr(v, q)
         else
            pos = ibset(v, q)
         end if
         if (performance(pos) > performance(v)) then
            loc_max(v) = 0
            exit
         end if
      end do
   end do
"

loopfun <- cfunction(sig = signature(performance="numeric", loc_max="integer", n="integer"),
                     dim=c("(0:(2**n-1))", "(0:(2**n-1))", ""),
                     loopcode,
                     language="F95")

N <- 20
performance = runif(2^N, min=0, max=1)
system.time({
  floop <- loopfun(performance, rep(1, 2^N), N)
})
##  user  system elapsed
## 0.049   0.003   0.052

N <- 29
performance = runif(2^N, min=0, max=1)
system.time({
  floop <- loopfun(performance, rep(1, 2^N), N)
})
##   user  system elapsed
## 17.892   1.848  19.741

我不认为预先计算邻居会对此有所帮助，因为我猜测访问如此大型阵列的不同部分的比较是最耗时的部分。