Псевдослучайная картина с фиксированной плотностью и зоной исключения
Я хочу создать 2D набор из N точек (обычно 1e2 - 1e4) в квадрате со следующими ограничениями:
между всеми точками должно быть минимальное расстояние (зона исключения жесткого ядра)
количество точек, заполняющих квадрат, дается заранее (или близкая оценка), так как я хочу получить фиксированную плотность (я могу немного отрегулировать размер квадрата впоследствии, если это необходимо).
шаблон должен быть достаточно "случайным"
быстрое решение является предпочтительным
Раньше я использовал rStrauss в пакете spatstat, но никогда не мог понять, как надежно получить заданное количество баллов, и довольно часто эта функция останавливает мою машину на 10 минут, предположительно из-за слишком сложной задачи. Я предполагаю, что для этого может быть более подходящая функция.
## regular grid of 1e2 points in [-10, 10]^2
xy = expand.grid(x=seq(-10, 10, length=10), y=seq(-10, 10, length=10))
N = NROW(xy)
РЕДАКТИРОВАТЬ: как предложено в ответе
xyr = rSSI(r=0.1, N, win = owin(c(-10,10),c(-10,10)), N)
plot(xyr)
2 ответа
rSSI
также в пакете spatstat заботится о ваших проблемах, за исключением, возможно, скорости, в зависимости от ваших стандартов. У него жесткая дистанция торможения, и он достигнет заданного количества очков (или откажется от попыток - вы можете установить порог для отказа), и места размещения будут случайными. Я не думаю, что это особенно быстро, но я смог создать 1e6
точки в единичном квадрате с расстоянием торможения 1e-4
примерно через 30 секунд. Скорость и успех будут сильно зависеть от вашего расстояния запрета относительно количества очков.
Главным образом в качестве предлога, чтобы узнать больше о Rcpp, вот моя попытка сделать небольшую функцию для этого:
require(inline)
require(Rcpp)
randPoints = cxxfunction(signature(r_n='int', r_mindist='float', r_maxiter='int'), body =
'
using namespace std;
int n = as<int> (r_n);
float mindist = as<float> (r_mindist);
int maxiter = as<int> (r_maxiter);
RNGScope scope;
bool tooclose;
int iter;
NumericVector rands (2);
NumericMatrix points (n, 2);
NumericVector dist (2);
for (int i=0; i < n; i++) {
iter = 0;
do {
iter++;
tooclose = false;
rands = runif(2, 0, 1);
for (int idist=0; idist < i; idist++) {
dist = rands - points(idist, _);
dist = dist * dist;
if (sqrt(accumulate(dist.begin(), dist.end(), 0.0)) < mindist) {
tooclose = true;
break;
}
}
} while (tooclose & iter < maxiter);
if (iter == maxiter) {
Rprintf("%d iterations reached\\nOnly %d points found\\n", maxiter, i+1);
break;
}
NumericMatrix::Row target(points, i);
target = rands;
}
return(wrap(points));
'
, plugin='Rcpp')
Тогда вы можете использовать его как:
> x = randPoints(1000, 0.05, 10000)
10000 iterations reached
Only 288 points found
А вот и сюжет:
x = x[as.logical(rowMeans(x != 0)), ]
dev.new(width=4, height=4)
plot(x)