2016-06-20 2 views
2

Я вычисляю матрицу расстояний, которая будет использоваться для некоторой кластеризации. Размер матрицы расстояния приведет к 4.6 gb, поэтому для его вычисления мне нужен очень быстрый код!Нужно ускорить мою дистанционную функцию в R

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

Функция расстояния принимает в качестве входных данных 2 географических координат и 2 строки и вычисляет расстояние следующим образом:

require(Imap) 

mydist <- function (lat1,lon1,lingua1,lat2,lon2,lingua2,DT){ 
    delta=0.1 
    gamma=3 
    d=sqrt(delta*gdist(lon1,lat1,lon2,lat2)^2 + gamma*(DT[language1 %in% lingua1 & language2 %in%lingua2]$distance)^2) 
} 

Он считывает расстояние двух моих строк из data.table DT, где я хранимую все возможные расстояния моих строк

функция, которая выделяет матрицу является:

require(bigmemory) 

distmatrix <- function(twit2,DT){ 
    N=dim(twit2)[1] 
    distmat = big.matrix(N,N) 
    for(i in 1:N){ 
    for(j in 1:N){ 
     distmat[i,j]=mydist(twit2[i,]$longitude,twit2[i,]$latitude,twit2[i,]$language,twit2[j,]$longitude,twit2[j,]$latitude,twit2[j,]$language,DT) 
    } 
    } 
    return(distmat) 
} 

EDIT: Я работаю над Альтернативный способ, который заключается в использовании библиотеки (ископаемое), в которой реализован векторную версию геодезического вычисления расстояния

Кроме того, я переехал DT в СД2, который теперь является квадратная матрица

library(fossil) 

lingdist <- function(lang1,lang2, DT2){ 
    list=colnames(DT2) 
    i=which(list==lang1) 
    j=which(list==lang2) 
    return(DT2[i,j]) 
} 

distmatrix <- function(twit2,DT2){ 
    N=dim(twit2)[1] 
    long<-as.vector(twit2$longitude) 
    lat <-as.vector(twit2$latitude) 
    lang<-as.vector(twit2$language) 
    distmat = t(as.matrix(earth.dist(as.matrix(cbind(long,lat))))) 
    for(i in 1:N) { 
    for (j in i:N) { 
     distmat[i,j]=sqrt(distmat[i,j]*distmat[i,j] + lingdist(lang[i],lang[j],DT2)) 
    } 
    } 
    return(distmat) 
} 

Этот добился существенного ускорения, 20x, с «малым» вводом (до 5 тыс. строк), но не смог выделить distmat с моим полным фреймворком (строки 24 тыс.)

У вас есть идея, как его решить?

EDIT2: вот небольшая версия моих баз данных

dput(DT2[1:5,1:5]) 
structure(c(0, 0.808204378308, 0.873223132147, 0.885209298235, 
0.849854297278, 0.808204378308, 0, 0.881177373275, 0.854352223232, 
0.854317529225, 0.873223132147, 0.881177373275, 0, 0.834454614055, 
0.861541199715, 0.885209298235, 0.854352223232, 0.834454614055, 
0, 0.76583938666, 0.849854297278, 0.854317529225, 0.861541199715, 
0.76583938666, 0), .Dim = c(5L, 5L), .Dimnames = list(c("1", 
"2", "3", "4", "5"), c("AA.SEMITIC.ARABIC_GULF_SPOKEN", "Alt.TURKIC.TURKISH", 
"An.MESO-PHILIPPINE.TAGALOG", "IE.BALTIC.LITHUANIAN", "IE.CELTIC.WELSH" 
))) 

dput(twit4[1:40,]) 
structure(list(day = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), nil = c(28L, 
28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 71L, 71L, 
20L, 5L, 24L, 49L, 50L, 28L, 28L, 22L, 22L, 21L, 21L, 24L, 20L, 
20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 24L 
), longitude = c(9.2235078, 9.22355903, 9.22362504, 9.22318987, 
9.22355654, 9.22361992, 9.22348964, 9.22366317, 9.22383346, 9.2238841, 
9.22374533, 9.22351081, 9.1361611, 9.1361805, 9.2144687, 9.1871549, 
9.2504309, 9.14652258, 9.16928, 9.22321188, 9.22387642, 9.2237509, 
9.22372656, 9.22278207, 9.2225214, 9.2470243, 9.22405217, 9.22404052, 
9.22405638, 9.22396956, 9.22402622, 9.2239671, 9.2239646, 9.22400299, 
9.22400299, 9.22403204, 9.22396816, 9.22404027, 9.22407831, 9.246786 
), latitude = c(45.45206021, 45.45202558, 45.4523043, 45.45211746, 
45.45204048, 45.45232425, 45.45207132, 45.45205533, 45.45218499, 
45.45216514, 45.45220716, 45.45214255, 45.5053803, 45.5053559, 
45.4871762, 45.4539539, 45.4660934, 45.45278042, 45.455855, 45.45882439, 
45.46055371, 45.47414199, 45.47947343, 45.48080458, 45.48119442, 
45.4658805, 45.49167007, 45.49168084, 45.49160813, 45.49164877, 
45.49165014, 45.49163468, 45.49165405, 45.49169004, 45.49169004, 
45.49160814, 45.49164155, 45.49161845, 45.49160889, 45.4660437 
), language = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 4L, 4L, 1L, 8L, 4L, 4L, 8L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("AA.SEMITIC.ARABIC_GULF_SPOKEN", 
"AA.SEMITIC.HEBREW", "Alt.TURKIC.TURKISH", "An.MESO-PHILIPPINE.TAGALOG", 
"AuA.VIET-MUONG.VIETNAMESE", "IE.ARMENIAN.EASTERN_ARMENIAN", 
"IE.BALTIC.LATVIAN", "IE.BALTIC.LITHUANIAN", "IE.CELTIC.WELSH", 
"IE.GERMANIC.DANISH", "IE.GERMANIC.DUTCH", "IE.GERMANIC.ICELANDIC", 
"IE.GERMANIC.NORWEGIAN_BOKMAAL", "IE.GERMANIC.STANDARD_GERMAN", 
"IE.GERMANIC.SWEDISH", "IE.GREEK.GREEK", "IE.INDIC.HINDI", "IE.IRANIAN.PERSIAN", 
"IE.ROMANCE.FRENCH", "IE.ROMANCE.PORTUGUESE", "IE.ROMANCE.ROMANIAN", 
"IE.ROMANCE.SPANISH", "IE.SLAVIC.BOSNIAN", "IE.SLAVIC.BULGARIAN", 
"IE.SLAVIC.CROATIAN", "IE.SLAVIC.POLISH", "IE.SLAVIC.RUSSIAN", 
"IE.SLAVIC.SERBOCROATIAN", "IE.SLAVIC.SLOVAK", "IE.SLAVIC.SLOVENIAN", 
"IE.SLAVIC.UKRAINIAN", "Jap.JAPANESE.JAPANESE", "Kor.KOREAN.KOREAN", 
"Krt.KARTVELIAN.GEORGIAN", "Oth.CREOLES_AND_PIDGINS.HAITIAN_CREOLE", 
"ST.CHINESE.CANTONESE", "TK.KAM-TAI.THAI", "Ura.FINNIC.ESTONIAN", 
"Ura.FINNIC.FINNISH", "Ura.UGRIC.HUNGARIAN"), class = "factor")), .Names = c("day", 
"nil", "longitude", "latitude", "language"), row.names = c("2", 
"6", "7", "8", "13", "15", "16", "20", "25", "29", "30", "32", 
"84", "86", "195", "266", "322", "467", "495", "521", "524", 
"534", "542", "546", "550", "580", "624", "640", "668", "676", 
"679", "699", "742", "751", "754", "768", "779", "800", "803", 
"857"), class = "data.frame") 
+2

Один взломать будет просто использовать верхнюю треугольную матрицу. Вам действительно не нужно вычислять расстояние между точками (A и B), а затем между точками (B и A). Для этого просто измените значение 'for (j in 1: N)' to 'for (j in i: N)'. – abhiieor

+2

Вы должны сделать подмножество data.frame вне цикла, например, 'long <- twit2 $ longitude', а затем использовать' long [i] 'внутри цикла. Но, скорее всего, достаточная скорость может быть достигнута только с помощью скомпилированного кода, например, с использованием Rcpp. – Roland

+0

Знаете ли вы, что [сообщение] (http://stackoverflow.com/questions/2908822/speed-up-the-loop-operation-in-r)? – Christoph

ответ

0

Это, безусловно, самый быстрый альтернативы я нашел, на основе data.table. Для этого требуется, чтобы данные были организованы в длинном формате, так что каждая строка в вашем наборе данных имеет комбинацию изначального (long long) и пункта назначения (lat long).

4 простых шагов

# load library 
library(geosphere) 
library(data.table) 

### STEP 1. reshape your matrix with languange distances to long format 
setDT(DT)[, language := names(DT)] 
DT_long <- melt.data.table(DT, id.vars="language", variable.name = "language2", value.name = "lingdist") 


### STEP 2. Get all possible combinations of origin and destination in long format 
df <- expand.grid.df(twit4, twit4) 
names(df)[c(3,4,5,8,9,10)] <- c("long_orig", "lat_orig", "language", "long_dest", "lat_dest","language2") 


### STEP 3. Efficiently merge the two datasets 
setkey(DT_long, "language", "language2") 
setkey(df, "language", "language2") 
df <- df[DT_long, nomatch=0] 


### STEP 4. Calculate distances 
df[ , dist := lingdist + distGeo(matrix(c(long_orig, lat_orig), ncol = 2), 
           matrix(c(long_dest, lat_dest), ncol = 2))/1000] 

Это решение должно быть относительно быстро. Узким местом здесь является операция expand.grid.df, которая, вероятно, является частью кода, который занимает больше времени, особенно при работе с большими наборами данных, как в вашем случае. Я уверен, что должна быть более быстрая альтернатива expand.grid.df. Я буду обновлять этот ответ, когда найду его.

+0

это довольно интересное решение, но есть две проблемы: 1) мое определение расстояния различно, так как ti также является частью «лингвистического расстояния» // 2) с большим фреймом данных я не могу process expand.grid.df – user5609462

+0

Я обновил свой ответ, чтобы включить лингвистическое расстояние в операции. Вы хотели бы рассчитать два типа расстояний? Что касается вашей точки 2, 'expand.grid.df' действительно является узким местом здесь, но я уверен, что должна быть более быстрая альтернатива этому –

Смежные вопросы