векторизованное сопоставление значений с числовым допуском

Question

Я пытаюсь сопоставить значения в векторе numeri c с дискретным набором значений, но из-за ошибок округления они не обязательно являются точным совпадением. Рассмотрим следующий пример, преобразование сферических координат в декартовы и попытку сопоставить полученное расстояние с исходными значениями r,

phithetar <- expand.grid(theta=seq(0,pi,length=10),
                        phi = seq(0, 2*pi, length=10),
                        r = seq(1,10))

xyz <- transform(phithetar, x = r*sin(theta)*cos(phi),
                                 y = r*sin(theta)*sin(phi),
                                 z = r*cos(theta))

xyz <- transform(xyz,  newr = sqrt(x^2+y^2+z^2))

with(xyz, sum(newr != r))

Это ненулевое значение, так как операции приводят к потере точности.

Вот моя медленная реализация для сопоставления значений в пределах заданного допуска; итерация по x делает его медленным. Альтернативой было бы расширение всех комбинаций (x, y), но это потребовало бы большого объема памяти для больших векторов.

match_closest <- function(x, y, tol=1e-5){

  match_one <- function(.x) {
    difs <- abs(.x - y)
    best <- which.min(difs)
    if(difs[best] > tol) message(.x, " ain't so close to ", y[best])
    y[best]
  }
  vapply(x, match_one, 0)

}

rvalues <- unique(phithetar$r)
xyz <- transform(xyz,  matchedr = match_closest(newr, rvalues))

with(xyz, sum(matchedr != r))

Я пробовал другую идею с cut(),

match_closenough <- function(x, y, tol=1e-5){

  levs <- c(y-0.5*tol, y+0.5*tol)
  f <- cut(x, breaks=levs, include.lowest = TRUE)
  y[f]

}

, но это не совсем работает. Точно так же этот трюк с округлением

match_closenough <- function(x, y, tol=1e-3){

  dig <- -log10(tol)
  y[match(round(x, dig), round(y, dig), nomatch = NA)]

}

вроде работает и быстрее, но ощущается слишком сильно fr agile.

Answer 1

Я думаю, что версию cut() можно заставить работать с некоторыми поправками (хотя факторы всегда кажутся вне моего контроля),

match_madethecut <- function(x, y, tol=1e-5){

  breaks <- sort(c(y-0.5*tol, y+0.5*tol))
  f <- cut(x, breaks=breaks, include.lowest = TRUE)
  rep(sort(y),each=2)[as.numeric(f)]

}

Обратите внимание, что y должен содержать уникальные значения, на практике. Это кажется быстрым, как и ожидалось, поскольку cut() выполняет сравнения на уровне C.

Answer 2

Если ваша цель - просто сопоставить значения в пределах некоторого (очень) малого допуска для учета потери точности во время последовательных операций с плавающей запятой (как в примере), то я думаю, что ваше последнее match_closenough() решение элегантен и быстр, и мне не кажется, что он особенно "fr agile".

Единственным возможным ограничением этого подхода в целом кажется то, что вы ограничены округлением до десятичной степени (а не произвольно выбранными tol значениями), и что вы не можете легко использовать tol значения больше 1. Оба этих случая, конечно, не применимы к вашему примеру, но могут быть ситуации, когда они могут стать актуальными.

Итак: что, если мы настаиваем на выборе произвольных tol значений и хотим иметь векторизованную функцию? Одним из решений было бы записать его в R Cpp (и это, вероятно, не было бы длиннее, чем решение только для R ...).

Но это может быть сделано «с нуля» в R: это решение использует несколько шагов, чтобы никогда не требовать невекторизованную функцию:

# define sloppy_match:
sloppy_match=function(x,y,tol=.Machine$double.eps,nomatch=NA) {
    # combine x & y indices:
        v=c(x,y)
        i=c(rep(1,length(x)),rep(0,length(y)))[order(v)]
        xi=cumsum(i==1)
        yi=cumsum(i==0)
    # find y values above and below each x value
        xrank=rank(x,ties="first")
        yi_below=yi[match(xrank,xi,NA)]
        yi_above=yi_below+1
        yi_below[yi_below==0]=NA
        yi_above[yi_above>length(y)]=NA
        ysort=sort(y)   
        y_below=ysort[yi_below]
        y_above=ysort[yi_above]
    # now choose closest y value to each x:
        closest_y=y_below
        to_change=is.na(y_below)|(x-closest_y>y_above-x)&!is.na(y_above)
        closest_y[to_change]=y_above[to_change]
    # and reject any that aren't close enough
        closest_y[abs(x-closest_y)>tol]=nomatch
    return(closest_y)
}

Попробуйте:

vec1=sample(1:10,5000,replace=T)+rnorm(50)*1e-6
vec2=sample(1:10,5000,replace=T)

test=sloppy_match(vec1,vec2,tol=1e-3)
all(abs(test-vec1)<1e-3)
# TRUE

Как это сравнить по скорости?

library(microbenchmark)
microbenchmark(match(vec1,vec1),
    match_closest(vec1,vec2,tol=1e-5),
    match_closenough(vec1,vec2,tol=1e-5),
    sloppy_match(vec1,vec2,tol=1e-5),
    times=10)
# Unit: microseconds
#                                       expr        min         lq        mean      median         uq        max neval
#                          match(vec1, vec1)    190.672    191.683    217.7785    199.9135    219.783    357.662    10
#     match_closest(vec1, vec2, tol = 1e-05) 186051.018 205712.022 208032.1262 209191.3980 215335.356 221682.351    10
#  match_closenough(vec1, vec2, tol = 1e-05)    511.426    513.772    531.2815    530.7850    544.416    561.935    10
#      sloppy_match(vec1, vec2, tol = 1e-05)   4110.814   4279.376   5859.7425   4301.0610   4608.116  18587.816    10

Таким образом, встроенный match(), очевидно, является самым быстрым, когда значения точно равны, и для достаточно близкого совпадения, если округление до степени меньше единицы допустимо ( и я уверен, что так будет в 99% случаев), match_closenough() самый быстрый.

Наконец, подход sloppy_match(), позволяющий сопоставить произвольно выбранные допуски, несколько медленнее, чем эти, но преимущество векторизации очевидно: он уже почти в 40 раз быстрее, чем не векторизованный подход match_closest(), всего с 5000 элементами, и эта разница становится еще более заметной, когда количество элементов увеличивается.

Итак, в конце концов, вы, вероятно, захотите выбрать один из двух подходов в зависимости от конкретного требования к допуску; но оба они будут в сотни раз быстрее, чем невекторизованный подход для больших векторов.