Oferuję Rrozwiązanie, które jest nieco zakodowane w Rcelu zilustrowania, jak można podejść do niego na innych platformach.
Problemem R(podobnie jak niektórych innych platform, szczególnie tych, które preferują funkcjonalny styl programowania) jest to, że ciągła aktualizacja dużej tablicy może być bardzo droga. Zamiast tego algorytm zachowuje własną strukturę prywatnych danych, w której (a) wymienione są wszystkie komórki, które zostały dotychczas wypełnione, i (b) wszystkie komórki, które można wybrać (na obwodzie wypełnionych komórek) są wymienione. Mimo że manipulowanie tą strukturą danych jest mniej wydajne niż bezpośrednie indeksowanie do tablicy, utrzymując zmodyfikowane dane w małym rozmiarze, prawdopodobnie zajmie to znacznie mniej czasu obliczeniowego. (Nie podjęto również żadnych wysiłków, aby go zoptymalizować R. Wstępna alokacja wektorów stanu powinna zaoszczędzić trochę czasu wykonania, jeśli wolisz nadal pracować wewnątrz R).
Kod jest komentowany i powinien być łatwy do odczytania. Aby algorytm był jak najbardziej kompletny, nie wykorzystuje żadnych dodatków poza końcem do wykreślenia wyniku. Jedyną trudną częścią jest to, że dla wydajności i prostoty woli indeksować do siatek 2D za pomocą indeksów 1D. Konwersja zachodzi w neighborsfunkcji, która potrzebuje indeksowania 2D, aby dowiedzieć się, czym mogą być dostępni sąsiedzi komórki, a następnie przekształca je w indeks 1D. Ta konwersja jest standardowa, więc nie będę komentować jej dalej, poza zaznaczeniem, że na innych platformach GIS możesz chcieć odwrócić role indeksów kolumn i wierszy. (W Rindeksy wierszy zmieniają się przed indeksami kolumn).
Aby to zilustrować, ten kod pobiera siatkę xprzedstawiającą ląd i rzekę cechę niedostępnych punktów, rozpoczyna się w określonej lokalizacji (5, 21) na tej siatce (w pobliżu dolnego zakola rzeki) i rozwija ją losowo, aby objąć 250 punktów . Całkowity czas wynosi 0,03 sekundy. (Gdy rozmiar tablicy zostanie zwiększony o współczynnik od 10 000 do 3000 wierszy o 5000 kolumn, czas wzrośnie tylko do 0,09 sekundy - współczynnik tylko około 3 - co pokazuje skalowalność tego algorytmu.) Zamiast po prostu wyprowadzając siatkę zer, 1 i 2, wyświetla sekwencję, z którą przydzielono nowe komórki. Na rycinie najwcześniejsze komórki są zielone, przechodząc przez złoto do kolorów łososiowych.
Powinno być oczywiste, że stosuje się ośmiopunktowe sąsiedztwo każdej komórki. W przypadku innych dzielnic wystarczy zmodyfikować nbrhoodwartość na początku expand: jest to lista przesunięć indeksu względem dowolnej komórki. Na przykład sąsiedztwo „D4” można określić jako matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).
Oczywiste jest również, że ta metoda rozprzestrzeniania ma swoje problemy: pozostawia dziury w tyle. Jeśli nie było to zamierzone, istnieją różne sposoby rozwiązania tego problemu. Na przykład trzymaj dostępne komórki w kolejce, aby najwcześniej znalezione komórki były również najwcześniej wypełnione. Nadal można zastosować pewną randomizację, ale dostępne komórki nie będą już wybierane z jednakowymi (równymi) prawdopodobieństwami. Innym, bardziej skomplikowanym sposobem byłby wybór dostępnych komórek z prawdopodobieństwem, które zależą od liczby ich wypełnionych sąsiadów. Gdy komórka zostanie otoczona, jej szansa na wybranie jest tak wysoka, że niewiele otworów pozostanie niewypełnionych.
Zakończę komentując, że nie jest to całkiem automat komórkowy (CA), który nie przechodziłby komórka po komórce, ale zamiast tego aktualizowałby całe warstwy komórek w każdym pokoleniu. Różnica jest subtelna: w przypadku CA prawdopodobieństwo wyboru komórek nie byłoby jednolite.
#
# Expand a patch randomly within indicator array `x` (1=unoccupied) by
# `n.size` cells beginning at index `start`.
#
expand <- function(x, n.size, start) {
if (x[start] != 1) stop("Attempting to begin on an unoccupied cell")
n.rows <- dim(x)[1]
n.cols <- dim(x)[2]
nbrhood <- matrix(c(-1,-1, -1,0, -1,1, 0,-1, 0,1, 1,-1, 1,0, 1,1), nrow=2)
#
# Adjoin one more random cell and update `state`, which records
# (1) the immediately available cells and (2) already occupied cells.
#
grow <- function(state) {
#
# Find all available neighbors that lie within the extent of `x` and
# are unoccupied.
#
neighbors <- function(i) {
n <- c((i-1)%%n.rows+1, floor((i-1)/n.rows+1)) + nbrhood
n <- n[, n[1,] >= 1 & n[2,] >= 1 & n[1,] <= n.rows & n[2,] <= n.cols,
drop=FALSE] # Remain inside the extent of `x`.
n <- n[1,] + (n[2,]-1)*n.rows # Convert to *vector* indexes into `x`.
n <- n[x[n]==1] # Stick to valid cells in `x`.
n <- setdiff(n, state$occupied)# Remove any occupied cells.
return (n)
}
#
# Select one available cell uniformly at random.
# Return an updated state.
#
j <- ceiling(runif(1) * length(state$available))
i <- state$available[j]
return(list(index=i,
available = union(state$available[-j], neighbors(i)),
occupied = c(state$occupied, i)))
}
#
# Initialize the state.
# (If `start` is missing, choose a value at random.)
#
if(missing(start)) {
indexes <- 1:(n.rows * n.cols)
indexes <- indexes[x[indexes]==1]
start <- sample(indexes, 1)
}
if(length(start)==2) start <- start[1] + (start[2]-1)*n.rows
state <- list(available=start, occupied=c())
#
# Grow for as long as possible and as long as needed.
#
i <- 1
indices <- c(NA, n.size)
while(length(state$available) > 0 && i <= n.size) {
state <- grow(state)
indices[i] <- state$index
i <- i+1
}
#
# Return a grid of generation numbers from 1, 2, ... through n.size.
#
indices <- indices[!is.na(indices)]
y <- matrix(NA, n.rows, n.cols)
y[indices] <- 1:length(indices)
return(y)
}
#
# Create an interesting grid `x`.
#
n.rows <- 3000
n.cols <- 5000
x <- matrix(1, n.rows, n.cols)
ij <- sapply(1:n.cols, function(i)
c(ceiling(n.rows * 0.5 * (1 + exp(-0.5*i/n.cols) * sin(8*i/n.cols))), i))
x[t(ij)] <- 0; x[t(ij - c(1,0))] <- 0; x[t(ij + c(1,0))] <- 0
#
# Expand around a specified location in a random but reproducible way.
#
set.seed(17)
system.time(y <- expand(x, 250, matrix(c(5, 21), 1)))
#
# Plot `y` over `x`.
#
library(raster)
plot(raster(x[n.rows:1,], xmx=n.cols, ymx=n.rows), col=c("#2020a0", "#f0f0f0"))
plot(raster(y[n.rows:1,] , xmx=n.cols, ymx=n.rows),
col=terrain.colors(255), alpha=.8, add=TRUE)
Przy niewielkich modyfikacjach możemy zapętlić się, expandaby utworzyć wiele klastrów. Wskazane jest rozróżnienie klastrów według identyfikatora, który tutaj będzie działał 2, 3, ... itd.
Najpierw zmień, expandaby zwracać (a) NAw pierwszym wierszu, jeśli wystąpił błąd, i (b) wartości indiceszamiast w macierzy y. (Nie marnuj czasu na tworzenie nowej matrycy yprzy każdym wywołaniu.) Po wprowadzeniu tej zmiany zapętlenie jest łatwe: wybierz losowy początek, spróbuj rozszerzyć się wokół niego, akumuluj indeksy klastrów, indicesjeśli się powiedzie, i powtarzaj, aż się zakończy. Kluczową częścią pętli jest ograniczenie liczby iteracji w przypadku, gdy nie można znaleźć wielu sąsiadujących klastrów: odbywa się to za pomocą count.max.
Oto przykład, w którym 60 centrów skupień jest wybieranych losowo równomiernie.
size.clusters <- 250
n.clusters <- 60
count.max <- 200
set.seed(17)
system.time({
n <- n.rows * n.cols
cells.left <- 1:n
cells.left[x!=1] <- -1 # Indicates occupancy of cells
i <- 0
indices <- c()
ids <- c()
while(i < n.clusters && length(cells.left) >= size.clusters && count.max > 0) {
count.max <- count.max-1
xy <- sample(cells.left[cells.left > 0], 1)
cluster <- expand(x, size.clusters, xy)
if (!is.na(cluster[1]) && length(cluster)==size.clusters) {
i <- i+1
ids <- c(ids, rep(i, size.clusters))
indices <- c(indices, cluster)
cells.left[indices] <- -1
}
}
y <- matrix(NA, n.rows, n.cols)
y[indices] <- ids
})
cat(paste(i, "cluster(s) created.", sep=" "))
Oto wynik po nałożeniu na siatkę 310 na 500 (wykonany wystarczająco mały i gruby, aby klastry były widoczne). Wykonanie zajmuje dwie sekundy; na siatce 3100 na 5000 (100 razy większa) zajmuje więcej czasu (24 sekundy), ale czas jest odpowiednio skalowany. (Na innych platformach, takich jak C ++, taktowanie nie powinno zależeć od rozmiaru siatki.)
