Случайная выборка символьного вектора без элементов, префиксных друг другу
Рассмотрим вектор символов, pool
, элементы которого (с нулевым заполнением) двоичные числа с цифрами max_len
.
max_len <- 4
pool <- unlist(lapply(seq_len(max_len), function(x)
do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
pool
## [1] "0" "1" "00" "10" "01" "11" "000" "100" "010" "110"
## [11] "001" "101" "011" "111" "0000" "1000" "0100" "1100" "0010" "1010"
## [21] "0110" "1110" "0001" "1001" "0101" "1101" "0011" "1011" "0111" "1111"
Я хотел бы попробовать n
этих элементов с ограничением, что ни один из выбранных элементов не является префиксами любого из другие элементы сэмплирования (т.е. если мы проецируем 1101
, мы запрещаем 1
, 11
и 110
, тогда как если мы опробуем 1
, мы запретим те элементы, начинающиеся с 1
, такие как 10
, 11
, 100
и т.д.).
Ниже приведена моя попытка использовать while
, но, конечно, это медленно, когда n
велико (или когда он приближается к 2^max_len
).
set.seed(1)
n <- 10
chosen <- sample(pool, n)
while(any(rowSums(outer(paste0('^', chosen), chosen, Vectorize(grepl))) > 1)) {
prefixes <- rowSums(outer(paste0('^', chosen), chosen, Vectorize(grepl))) > 1
pool <- pool[rowSums(Vectorize(grepl, 'pattern')(
paste0('^', chosen[!prefixes]), pool)) == 0]
chosen <- c(chosen[!prefixes], sample(pool, sum(prefixes)))
}
chosen
## [1] "0100" "0101" "0001" "0011" "1000" "111" "0000" "0110" "1100" "0111"
Это можно немного улучшить, изначально удалив из pool
те элементы, включение которых означает, что в pool
осталось меньше элементов, чтобы взять общий образец размера n
. Например, когда max_len = 4
и n > 9
мы можем сразу удалить 0
и 1
из pool
, так как, включив либо, максимальный образец будет равен 9 (либо 0
, и восемь 4-символьных элементов начиная с 1
или 1
и восьми 4-символьных элементов, начинающихся с 0
).
На основе этой логики мы могли бы затем опустить элементы из pool
, прежде чем брать исходный образец, с чем-то вроде:
pool <- pool[
nchar(pool) > tail(which(n > (2^max_len - rev(2^(0:max_len))[-1] + 1)), 1)]
Может ли кто-нибудь подумать о лучшем подходе? Я чувствую, что я пропускаю нечто гораздо более простое.
ИЗМЕНИТЬ
Чтобы прояснить мои намерения, я опишу пул как набор ветвей, где узлы и концы являются узлами (элементами pool
). Предположим, что желтый node на следующем рисунке (т.е. 010) был нарисован. Теперь вся красная "ветвь", состоящая из узлов 0, 01 и 010, удаляется из пула. Это то, что я имел в виду, запрещая выборку узлов, которые "префиксные" узлы уже находятся в нашем примере (а также те, которые уже предваряются узлами в нашем примере).
![enter image description here]()
Если выборка node была на полпути вдоль ветки, например, 01 на следующем рисунке, то все красные узлы (0, 01, 010 и 011) запрещены, так как префиксы 0 и 01 префиксов и 010, и 011.
![enter image description here]()
Я не хочу пробовать ни 1 или 0 на каждом соединении (т.е. ходить вдоль ветвей, переворачивающих монеты на вилках) - это хорошо, чтобы иметь как в образце, так и: (1) родители (или грандиозные, родители и т.д.), или дети (внуки и т.д.) из node уже не отбираются; и (2) при выборке node останутся достаточные узлы для достижения желаемого образца размером n
.
На втором рисунке выше, если 010 был первым выбором, все узлы в черных узлах все еще (в настоящее время) действительны, предполагая n <= 4
. Например, если n==4
и мы выбрали node 1 next (и поэтому наши выборки теперь включали 01 и 1), мы впоследствии запретили бы node 00 (из-за правила 2 выше), но все же могли бы выбрать 000 и 001, давая нам наш 4-элементный образец. Если n==5
, с другой стороны, node 1 будет отменено на этом этапе.
Ответы
Ответ 1
Введение
Это числовое изменение строкового алгоритма, реализованного в другом ответе. Это быстрее и не требует создания или сортировки пула.
Схема алгоритма
Мы можем использовать целые числа для представления ваших двоичных строк, что значительно упрощает задачу создания пула и последовательного устранения значений. Например, при max_len==3
мы можем взять число 1--
(где -
представляет отступы), чтобы означать 4
в десятичной форме. Кроме того, мы можем определить, что числа, требующие устранения, если мы выберем это число, - это числа между 4
и 4 + 2 ^ x - 1
. Здесь x
- количество элементов заполнения (в этом случае 2), поэтому числа, подлежащие исключению, находятся между 4
и 4 + 2 ^ 2 - 1
(или между 4
и 7
, выраженным как 100
, 110
и 111
).
Чтобы точно соответствовать вашей проблеме, нам нужна небольшая морщина, так как вы обрабатываете числа, которые потенциально одинаковы в двоичном формате, как отдельные для некоторых частей вашего алгоритма. Например, 100
, 10-
и 1--
- все одинаковое число, но в вашей схеме нужно обрабатывать по-разному. В мире max_len==3
у нас есть 8 возможных чисел, но возможны 14 представлений:
0 - 000: 0--, 00-
1 - 001:
2 - 010: 01-
3 - 011:
4 - 100: 1--, 10-
5 - 101:
6 - 110: 11-
7 - 111:
Итак, 0 и 4 имеют три возможных кодировки, 2 и 6 имеют два, а все остальные - только одно. Нам нужно создать целочисленный пул, который представляет более высокую вероятность выбора для чисел с множественным представлением, а также механизм отслеживания количества пробелов, которые включает число. Мы можем сделать это, добавив несколько бит в конце номера, чтобы указать весы, которые мы хотим. Итак, наши числа становятся (мы используем здесь два бита):
jbaum | int | bin | bin.enc | int.enc
0-- | 0 | 000 | 00000 | 0
00- | 0 | 000 | 00001 | 1
000 | 0 | 000 | 00010 | 2
001 | 1 | 001 | 00100 | 3
01- | 2 | 010 | 01000 | 4
010 | 2 | 010 | 01001 | 5
011 | 3 | 011 | 01101 | 6
1-- | 4 | 100 | 10000 | 7
10- | 4 | 100 | 10001 | 8
100 | 4 | 100 | 10010 | 9
101 | 5 | 101 | 10100 | 10
11- | 6 | 110 | 11000 | 11
110 | 6 | 110 | 11001 | 12
111 | 7 | 111 | 11100 | 13
Несколько полезных свойств:
-
enc.bits
представляет количество бит, которое нам нужно для кодирования (два в этом случае)
-
int.enc %% enc.bits
указывает, сколько ячеек явно указано
-
int.enc %/% enc.bits
возвращает int
-
int * 2 ^ enc.bits + explicitly.specified
возвращает int.enc
Обратите внимание, что explicitly.specified
здесь находится между 0
и max_len - 1
в нашей реализации, поскольку всегда указана не менее одной цифры. Теперь у нас есть кодировка, которая полностью отражает вашу структуру данных, используя только целые числа. Мы можем выбирать из целых чисел и воспроизводить желаемые результаты с правильными весами и т.д. Одним из ограничений этого подхода является то, что мы используем 32-битные целые числа в R, и мы должны зарезервировать некоторые биты для кодирования, поэтому мы ограничиваемся пулами max_len==25
или около того. Вы могли бы пойти больше, если бы использовали целые числа, заданные с плавающей точкой с двойной точностью, но мы этого не делали.
Избегание дубликатов сообщений
Есть два грубых способа гарантировать, что мы не выбираем одно и то же значение дважды
- Отслеживайте, какие значения остаются доступными для выбора, и произвольно выбирайте из них
- Образец случайным образом из всех возможных значений, а затем проверьте, было ли это значение выбрано ранее, и если оно имеет, снова образец
В то время как первый вариант кажется самым чистым, он на самом деле очень дорогостоящий вычислительный. Это требует либо выполнения векторного сканирования всех возможных значений для каждого выбора, чтобы предварительно дисквалифицировать выбранные значения, либо создать сокращающий вектор, содержащий недисквалифицированные значения. Опция усадки является более эффективной, чем векторное сканирование, если вектор сжимается по ссылке через код C, но даже тогда он требует повторных переводов потенциально больших частей вектора, и для этого требуется C.
Здесь мы используем метод # 2. Это позволяет нам случайным образом перетасовывать всевозможные возможные значения один раз, а затем выбирать каждое значение последовательно, проверять, что оно не было дисквалифицировано, а если оно есть, выбрать другое и т.д. Это работает, потому что тривиально проверять, значение было выбрано в результате нашей кодировки значений; мы можем определить местоположение значения в отсортированной таблице на основе только значения. Таким образом, мы записываем статус каждого значения в отсортированную таблицу и можем либо обновлять, либо искать это состояние через прямой доступ к индексу (не требуется проверка).
Примеры
Реализация этого алгоритма в базе R доступна как сущность. Эта конкретная реализация только тянет полные ничьи. Вот пример из 10 рисунков из 8 элементов из пула max_len==4
:
# each column represents a draw from a `max_len==4` pool
set.seed(6); replicate(10, sample0110b(4, 8))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] "1000" "1" "0011" "0010" "100" "0011" "0" "011" "0100" "1011"
[2,] "111" "0000" "1101" "0000" "0110" "0100" "1000" "00" "0101" "1001"
[3,] "0011" "0110" "1001" "0100" "0000" "0101" "1101" "1111" "10" "1100"
[4,] "0100" "0010" "0000" "0101" "1101" "101" "1011" "1101" "0110" "1101"
[5,] "101" "0100" "1100" "1100" "0101" "1001" "1001" "1000" "1111" "1111"
[6,] "110" "0111" "1011" "111" "1011" "110" "1111" "0100" "0011" "000"
[7,] "0101" "0101" "111" "011" "1010" "1000" "1100" "101" "0001" "0101"
[8,] "011" "0001" "01" "1010" "0011" "1110" "1110" "1001" "110" "1000"
Мы также первоначально имели две реализации, которые полагались на метод # 1, чтобы избежать дублирования, один в базе R и один на C, но даже версия C не так быстро, как новая базовая версия R, когда n
большой. Эта функция реализует способность рисовать неполные рисунки, поэтому мы приводим их здесь для справки:
Сравнительные тесты
Вот набор тестов, сравнивающих несколько функций, которые отображаются в этом Q/A. Время в миллисекундах. Версия brodie.b
- это версия, описанная в этом ответе. brodie
является исходной реализацией, brodie.C
является оригиналом с некоторым C. Все это обеспечивает соблюдение требования для полных образцов. brodie.str
- это строковая версия в другом ответе.
size n jbaum josilber frank tensibai brodie.b brodie brodie.C brodie.str
1 4 10 11 1 3 1 1 1 1 0
2 4 50 - - - 1 - - - 1
3 4 100 - - - 1 - - - 0
4 4 256 - - - 1 - - - 1
5 4 1000 - - - 1 - - - 1
6 8 10 1 290 6 3 2 2 1 1
7 8 50 388 - 8 8 3 4 3 4
8 8 100 2,506 - 13 18 6 7 5 5
9 8 256 - - 22 27 13 14 12 6
10 8 1000 - - - 27 - - - 7
11 16 10 - - 615 688 31 61 19 424
12 16 50 - - 2,123 2,497 28 276 19 1,764
13 16 100 - - 4,202 4,807 30 451 23 3,166
14 16 256 - - 11,822 11,942 40 1,077 43 8,717
15 16 1000 - - 38,132 44,591 83 3,345 130 27,768
Это довольно хорошо масштабируется для больших пулов
system.time(sample0110b(18, 100000))
user system elapsed
8.441 0.079 8.527
Примечания:
- frank, и brodie (минус brodie.str) не требуют предварительного создания пулов, что повлияет на сравнение (см. ниже).
- Josilber - это версия LP
- jbaum - пример OP
- tensibai слегка изменен для выхода вместо отказа, если пул пуст.
- не настроен для запуска python, поэтому не может сравниться/учетная запись для буферизации
-
-
представляют собой либо неосуществимые варианты, либо слишком медленные по времени разумно.
Тайминга не включают рисование пулов (0.8
, 2.5
, 401
миллисекунд соответственно для размера 4
, 8
и 16
)), что необходимо для jbaum
, josilber
и brodie.str
запускает или сортирует их (0.1
, 2.7
, 3700
миллисекунды для размера 4
, 8
и 16
), что необходимо для brodie.str
в дополнение к розыгрышу. Если вы хотите включить эти или нет, зависит от того, сколько раз вы запускаете функцию для определенного пула. Кроме того, почти наверняка есть лучшие способы генерации/сортировки пулов.
Это среднее время трех прогонов с microbenchmark
. Код доступен как сущность, но учтите, что вам нужно будет загрузить sample0110b
, sample0110
, sample01101
и sample01
.
Ответ 2
Я нашел проблему интересной, поэтому я попробовал и получил это с очень низкими навыками в R (поэтому это может быть улучшено):
Новая отредактированная версия, благодаря @Franck:
library(microbenchmark)
library(lineprof)
max_len <- 16
pool <- unlist(lapply(seq_len(max_len), function(x)
do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
n<-100
library(stringr)
tree_sample <- function(samples,pool) {
results <- vector("integer",samples)
# Will be used on a regular basis, compute it in advance
PoolLen <- str_length(pool)
# Make a mask vector based on the length of each entry of the pool
masks <- strtoi(str_pad(str_pad("1",PoolLen,"right","1"),max_len,"right","0"),base=2)
# Make an integer vector from "0" right padded orignal: for max_len=4 and pool entry "1" we get "1000" => 8
# This will allow to find this entry as parent of 10 and 11 which become "1000" and "1100", as integer 8 and 12 respectively
# once bitwise "anded" with the repective mask "1000" the first bit is striclty the same, so it a parent.
integerPool <- strtoi(str_pad(pool,max_len,"right","0"),base=2)
# Create a vector to filter the available value to sample
ok <- rep(TRUE,length(pool))
#Precompute the result of the bitwise and betwwen our integer pool and the masks
MaskedPool <- bitwAnd(integerPool,masks)
while(samples) {
samp <- sample(pool[ok],1) # Get a sample
results[samples] <- samp # Store it as result
ok[pool == samp] <- FALSE # Remove it from available entries
vsamp <- strtoi(str_pad(samp,max_len,"right","0"),base=2) # Get the integer value of the "0" right padded sample
mlen <- str_length(samp) # Get sample len
#Creation of unitary mask to remove childs of sample
mask <- strtoi(paste0(rep(1:0,c(mlen,max_len-mlen)),collapse=""),base=2)
# Get the result of bitwise And between the integerPool and the sample mask
FilterVec <- bitwAnd(integerPool,mask)
# Get the bitwise and result of the sample and it mask
Childm <- bitwAnd(vsamp,mask)
ok[FilterVec == Childm] <- FALSE # Remove from available entries the childs of the sample
ok[MaskedPool == bitwAnd(vsamp,masks)] <- FALSE # compare the sample with all the masks to remove parents matching
samples <- samples -1
}
print(results)
}
microbenchmark(tree_sample(n,pool),times=10L)
Основная идея состоит в том, чтобы использовать сравнение битмаски, чтобы узнать, является ли один образец родителем другого (общая часть бит), если да подавить этот элемент из пула.
Теперь требуется 1,4 с, чтобы нарисовать 100 образцов из пула длиной 16 на моей машине.
Ответ 3
Вы можете отсортировать пул, чтобы решить, какие элементы дисквалифицировать. Например, глядя на сортированный пул из трех элементов:
[1] "0" "00" "000" "001" "01" "010" "011" "1" "10" "100" "101" "11"
[13] "110" "111"
Я могу сказать, что я могу дисквалифицировать все, что следует за моим выбранным элементом, который имеет больше символов, чем мой элемент, до первого элемента, который имеет одинаковое число или меньшее количество символов. Например, если я выбираю "01", я сразу вижу, что следующие два элемента ( "010", "011" ) нужно удалить, но не после этого, потому что "1" имеет меньше символов. Снятие "0" после этого легко. Вот реализация:
library(fastmatch) # could use `match`, but we repeatedly search against same hash
# `pool` must be sorted!
sample01 <- function(pool, n) {
picked <- logical(length(pool))
chrs <- nchar(pool)
pick.list <- character(n)
pool.seq <- seq_along(pool)
for(i in seq(n)) {
# Make sure pool not exhausted
left <- which(!picked)
left.len <- length(left)
if(!length(left)) break
# Sample from pool
seq.left <- seq.int(left)
pool.left <- pool[left]
chrs.left <- chrs[left]
pick <- sample(length(pool.left), 1L)
# Find all the elements with more characters that are disqualified
# and store their indices in `valid` (bad name...)
valid.tmp <- chrs.left > chrs.left[[pick]] & seq.left > pick
first.invalid <- which(!valid.tmp & seq.left > pick)
valid <- if(length(first.invalid)) {
pick:(first.invalid[[1L]] - 1L)
} else pick:left.len
# Translate back to original pool indices since we're working on a
# subset in `pool.left`
pool.seq.left <- pool.seq[left]
pool.idx <- pool.seq.left[valid]
val <- pool[[pool.idx[[1L]]]]
# Record the picked value, and all the disqualifications
pick.list[[i]] <- val
picked[pool.idx] <- TRUE
# Disqualify shorter matches
to.rem <- vapply(
seq.int(nchar(val) - 1), substr, character(1L), x=val, start=1L
)
to.rem.idx <- fmatch(to.rem, pool, nomatch=0)
picked[to.rem.idx] <- TRUE
}
pick.list
}
И функция создания отсортированных пулов (точно так же, как ваш код, но возвращает сортировку):
make_pool <- function(size)
sort(
unlist(
lapply(
seq_len(size),
function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))
) ) )
Затем, используя пул max_len
3 (полезный для визуального осмотра вещей, ведет себя так, как ожидалось):
pool3 <- make_pool(3)
set.seed(1)
sample01(pool3, 8)
# [1] "001" "1" "010" "011" "000" "" "" ""
sample01(pool3, 8)
# [1] "110" "111" "011" "10" "00" "" "" ""
sample01(pool3, 8)
# [1] "000" "01" "11" "10" "001" "" "" ""
sample01(pool3, 8)
# [1] "011" "101" "111" "001" "110" "100" "000" "010"
Обратите внимание, что в последнем случае мы получаем все трехзначные двоичные комбинации (2 ^ 3), потому что случайно мы сохранили выборку из трехзначных. Кроме того, с помощью всего лишь 3-х размерного пула есть много выборок, которые предотвращают полную 8 ничьих; вы можете обратиться к этому с вашим предложением об устранении комбинаций, которые предотвращают полную ничью из пула.
Это довольно быстро. Рассматривая пример max_len==9
, который занял 2 секунды с альтернативным решением:
pool9 <- make_pool(9)
microbenchmark(sample01(pool9, 4))
# Unit: microseconds
# expr min lq median uq max neval
# sample01(pool9, 4) 493.107 565.015 571.624 593.791 983.663 100
Около половины миллисекунды. Вы можете разумно попробовать довольно большие пулы:
pool16 <- make_pool(16) # 131K entries
system.time(sample01(pool16, 100))
# user system elapsed
# 3.407 0.146 3.552
Это не очень быстро, но мы говорим о пуле с элементами 130K. Существует также возможность дополнительной оптимизации.
Обратите внимание, что шаг сортировки становится относительно медленным для больших пулов, но я не считаю это, потому что вам нужно когда-либо делать это один раз, и вы, вероятно, могли бы разработать разумный алгоритм для создания предварительно отсортированного пула.
Существует также возможность более быстрого целочисленного подхода к двоичному основанию, который я изучил в теперь удаленном ответе, но для этого требуется более справедливая работа, чтобы привязать ее к тому, что вы ищете.
Ответ 4
Отображение идентификаторов в строках.. Вы можете сопоставить числа с вашими векторами 0/1, как @BrodieG упомянул:
# some key objects
n_pool = sum(2^(1:max_len)) # total number of indices
cuts = cumsum(2^(1:max_len-1)) # new group starts
inds_by_g = mapply(seq,cuts,cuts*2) # indices grouped by length
# the mapping to strings (one among many possibilities)
library(data.table)
get_01str <- function(id,max_len){
cuts = cumsum(2^(1:max_len-1))
g = findInterval(id,cuts)
gid = id-cuts[g]+1
data.table(g,gid)[,s:=
do.call(paste,c(list(sep=""),lapply(
seq(g[1]),
function(x) (gid-1) %/% 2^(x-1) %% 2
)))
,by=g]$s
}
Поиск идентификаторов для удаления. Мы последовательно отбрасываем id
из пула выборки:
# the mapping from one index to indices of nixed strings
get_nixstrs <- function(g,gid,max_len){
cuts = cumsum(2^(1:max_len-1))
gids_child = {
x = gid%%2^sequence(g-1)
ifelse(x,x,2^sequence(g-1))
}
ids_child = gids_child+cuts[sequence(g-1)]-1
ids_parent = if (g==max_len) gid+cuts[g]-1 else {
gids_par = vector(mode="list",max_len)
gids_par[[g]] = gid
for (gg in seq(g,max_len-1))
gids_par[[gg+1]] = c(gids_par[[gg]],gids_par[[gg]]+2^gg)
unlist(mapply(`+`,gids_par,cuts-1))
}
c(ids_child,ids_parent)
}
Индексы сгруппированы на g
, количество символов, nchar(get_01str(id))
. Поскольку индексы сортируются по g
, g=findInterval(id,cuts)
- более быстрый маршрут.
Индекс в группе g
, 1 < g < max_len
имеет один "дочерний" индекс размера g-1
и два родительских индекса размера g+1
. Для каждого дочернего node мы берем его дочерний элемент node, пока не нажмем g==1
; и для каждого родителя node мы берем их пару родительских узлов, пока не нажмем g==max_len
.
Структура дерева простейшая с точки зрения идентификатора внутри группы, gid
. gid
отображает два родителя, gid
и gid+2^g
; и реверсируя это сопоставление, вы обнаружите его.
Sampling
drawem <- function(n,max_len){
cuts = cumsum(2^(1:max_len-1))
inds_by_g = mapply(seq,cuts,cuts*2)
oklens = (1:max_len)[ n <= 2^max_len*(1-2^(-(1:max_len)))+1 ]
okinds = unlist(inds_by_g[oklens])
mysamp = rep(0,n)
for (i in 1:n){
id = if (length(okinds)==1) okinds else sample(okinds,1)
g = findInterval(id,cuts)
gid = id-cuts[g]+1
nixed = get_nixstrs(g,gid,max_len)
# print(id); print(okinds); print(nixed)
mysamp[i] = id
okinds = setdiff(okinds,nixed)
if (!length(okinds)) break
}
res <- rep("",n)
res[seq.int(i)] <- get_01str(mysamp[seq.int(i)],max_len)
res
}
Часть oklens
объединяет идею OP для исключения строк, которые гарантируют, что выборка невозможна. Однако, даже делая это, мы можем следовать пути выборки, который оставляет нас без дополнительных параметров. Принимая пример OP max_len=4
и n=10
, мы знаем, что мы должны отказаться от рассмотрения 0
и 1
, но что произойдет, если наши первые четыре розыгрыша будут 00
, 01
, 11
и 10
? Ну, я думаю, нам не повезло. Вот почему вы должны определить вероятности выборки. (У ОП есть другая идея, для определения на каждом шаге, какие узлы приведут к невозможному состоянию, но это кажется высоким порядком.)
Иллюстрация
# how the indices line up
n_pool = sum(2^(1:max_len))
pdt <- data.table(id=1:n_pool)
pdt[,g:=findInterval(id,cuts)]
pdt[,gid:=1:.N,by=g]
pdt[,s:=get_01str(id,max_len)]
# example run
set.seed(4); drawem(5,5)
# [1] "01100" "1" "0001" "0101" "00101"
set.seed(4); drawem(8,4)
# [1] "1100" "0" "111" "101" "1101" "100" "" ""
Тесты (старше, чем в ответе @BrodieG)
require(rbenchmark)
max_len = 8
n = 8
benchmark(
jos_lp = {
pool <- unlist(lapply(seq_len(max_len),
function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
sample.lp(pool, n)},
bro_string = {pool <- make_pool(max_len);sample01(pool,n)},
fra_num = drawem(n,max_len),
replications=5)[1:5]
# test replications elapsed relative user.self
# 2 bro_string 5 0.05 2.5 0.05
# 3 fra_num 5 0.02 1.0 0.02
# 1 jos_lp 5 1.56 78.0 1.55
n = 12
max_len = 12
benchmark(
bro_string={pool <- make_pool(max_len);sample01(pool,n)},
fra_num=drawem(n,max_len),
replications=5)[1:5]
# test replications elapsed relative user.self
# 1 bro_string 5 0.54 6.75 0.51
# 2 fra_num 5 0.08 1.00 0.08
Другие ответы. Есть еще два ответа:
jos_enum = {pool <- unlist(lapply(seq_len(max_len),
function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
get.template(pool, n)}
bro_num = sample011(max_len,n)
Я отказался от метода перечисления @josilber, потому что он занял слишком много времени; и @BrodieG числовой/индексный метод, потому что он не работал в то время, но теперь делает. См. Обновленный ответ @BrodieG для более точного тестирования.
Скорость против правильности. В то время как ответы @josilber намного медленнее (и для метода перечисления, по-видимому, намного более интенсивного в памяти), с первой попытки гарантируется выборка размера n
. С помощью строкового метода @BrodieG или этого ответа вам придется снова и снова повторять выбор в надежде на рисование полного n
. С большим max_len
, это не должно быть такой проблемой, я думаю.
Этот ответ масштабируется лучше, чем bro_string
, потому что он не требует построения pool
спереди.
Ответ 5
Он находится в питоне, а не в r, но jbaums сказал, что все будет хорошо.
Итак, вот мой вклад, см. комментарии в источнике для объяснения важнейших частей.
Я все еще работаю над аналитическим решением для определения количества возможных комбинаций c
для дерева глубин t
и S
образцов, поэтому я могу улучшить функцию combs
. Может, у кого-то это есть?
Это действительно узкое место сейчас.
Отбор 100 узлов из дерева с глубиной 16 занимает около 8 мс на моем ноутбуке. Не в первый раз, но он ускоряется до определенной точки, чем больше вы выбираете, потому что combBuffer заполняется.
import random
class Tree(object):
"""
:param level: The distance of this node from the root.
:type level: int
:param parent: This trees parent node
:type parent: Tree
:param isleft: Determines if this is a left or a right child node. Can be
omitted if this is the root node.
:type isleft: bool
A binary tree representing possible strings which match r'[01]{1,n}'. Its
purpose is to be able to sample n of its nodes where none of the sampled
nodes' ids is a prefix for another one.
It is possible to change Tree.maxdepth and then reuse the root. All
children are created ON DEMAND, which means everything is lazily evaluated.
If the Tree gets too big anyway, you can call 'prune' on any node to delete
its children.
>>> t = Tree()
>>> t.sample(8, toString=True, depth=3)
['111', '110', '101', '100', '011', '010', '001', '000']
>>> Tree.maxdepth = 2
>>> t.sample(4, toString=True)
['11', '10', '01', '00']
"""
maxdepth = 10
_combBuffer = {}
def __init__(self, level=0, parent=None, isleft=None):
self.parent = parent
self.level = level
self.isleft = isleft
self._left = None
self._right = None
@classmethod
def setMaxdepth(cls, depth):
"""
:param depth: The new depth
:type depth: int
Sets the maxdepth of the Trees. This basically is the depth of the root
node.
"""
if cls.maxdepth == depth:
return
cls.maxdepth = depth
@property
def left(self):
"""This tree left child, 'None' if this is a leave node"""
if self.depth == 0:
return None
if self._left is None:
self._left = Tree(self.level+1, self, True)
return self._left
@property
def right(self):
"""This tree right child, 'None' if this is a leave node"""
if self.depth == 0:
return None
if self._right is None:
self._right = Tree(self.level+1, self, False)
return self._right
@property
def depth(self):
"""
This tree depth. (maxdepth-level)
"""
return self.maxdepth-self.level
@property
def id(self):
"""
This tree id, string of '0 and '1 equal to the path from the root
to this subtree. Where '1' means going left and '0' means going right.
"""
# level 0 is the root node, it has no id
if self.level == 0:
return ''
# This takes at most Tree.maxdepth recursions. Therefore
# it is save to do it this way. We could also save each nodes
# id once it is created to avoid recreating it every time, however
# this won't save much time but use quite some space.
return self.parent.id + ('1' if self.isleft else '0')
@property
def leaves(self):
"""
The amount of leave nodes, this tree has. (2**depth)
"""
return 2**self.depth
def __str__(self):
return self.id
def __len__(self):
return 2*self.leaves-1
def prune(self):
"""
Recursively prune this tree children.
"""
if self._left is not None:
self._left.prune()
self._left.parent = None
self._left = None
if self._right is not None:
self._right.prune()
self._right.parent = None
self._right = None
def combs(self, n):
"""
:param n: The amount of samples to be taken from this tree
:type n: int
:returns: The amount of possible combinations to choose n samples from
this tree
Determines recursively the amount of combinations of n nodes to be
sampled from this tree.
Subsequent calls with same n on trees with same depth will return the
result from the previous computation rather than computing it again.
>>> t = Tree()
>>> Tree.maxdepth = 4
>>> t.combs(16)
1
>>> Tree.maxdepth = 3
>>> t.combs(6)
58
"""
# important for the amount of combinations is only n and the depth of
# this tree
key = (self.depth, n)
# We use the dict to save computation time. Calling the function with
# equal values on equal nodes just returns the alrady computed value if
# possible.
if key not in Tree._combBuffer:
leaves = self.leaves
if n < 0:
N = 0
elif n == 0 or self.depth == 0 or n == leaves:
N = 1
elif n == 1:
return (2*leaves-1)
else:
if n > leaves/2:
# if n > leaves/2, at least n-leaves/2 have to stay on
# either side, otherweise the other one would have to
# sample more nodes than possible.
nMin = n-leaves/2
else:
nMin = 0
# The rest n-2*nMin is the amount of samples that are free to
# fall on either side
free = n-2*nMin
N = 0
# sum up the combinations of all possible splits
for addLeft in range(0, free+1):
nLeft = nMin + addLeft
nRight = n - nLeft
N += self.left.combs(nLeft)*self.right.combs(nRight)
Tree._combBuffer[key] = N
return N
return Tree._combBuffer[key]
def sample(self, n, toString=False, depth=None):
"""
:param n: How may samples to take from this tree
:type n: int
:param toString: If 'True' result will direclty be turned into a list
of strings
:type toString: bool
:param depth: If not None, will overwrite Tree.maxdepth
:type depth: int or None
:returns: List of n nodes sampled from this tree
:throws ValueError: when n is invalid
Takes n random samples from this tree where none of the sample ids is
a prefix for another one's.
For an example see Tree docstring.
"""
if depth is not None:
Tree.setMaxdepth(depth)
if toString:
return [str(e) for e in self.sample(n)]
if n < 0:
raise ValueError('Negative sample size is not possible!')
if n == 0:
return []
leaves = self.leaves
if n > leaves:
raise ValueError(('Cannot sample {} nodes, with only {} ' +
'leaves!').format(n, leaves))
# Only one sample to choose, that is nice! We are free to take any node
# from this tree, including this very node itself.
if n == 1 and self.level > 0:
# This tree has 2*leaves-1 nodes, therefore
# the probability that we keep the root node has to be
# 1/(2*leaves-1) = P_root. Lets create a random number from the
# interval [0, 2*leaves-1).
# It will be 0 with probability 1/(2*leaves-1)
P_root = random.randint(0, len(self)-1)
if P_root == 0:
return [self]
else:
# The probability to land here is 1-P_root
# A child tree size is (leaves-1) and since it obeys the same
# rule as above, the probability for each of its nodes to
# 'survive' is 1/(leaves-1) = P_child.
# However all nodes must have equal probability, therefore to
# make sure that their probability is also P_root we multiply
# them by 1/2*(1-P_root). The latter is already done, the
# former will be achieved by the next condition.
# If we do everything right, this should hold:
# 1/2 * (1-P_root) * P_child = P_root
# Lets see...
# 1/2 * (1-1/(2*leaves-1)) * (1/leaves-1)
# (1-1/(2*leaves-1)) * (1/(2*(leaves-1)))
# (1-1/(2*leaves-1)) * (1/(2*leaves-2))
# (1/(2*leaves-2)) - 1/((2*leaves-2) * (2*leaves-1))
# (2*leaves-1)/((2*leaves-2) * (2*leaves-1)) - 1/((2*leaves-2) * (2*leaves-1))
# (2*leaves-2)/((2*leaves-2) * (2*leaves-1))
# 1/(2*leaves-1)
# There we go!
if random.random() < 0.5:
return self.right.sample(1)
else:
return self.left.sample(1)
# Now comes the tricky part... n > 1 therefore we are NOT going to
# sample this node. Its probability to be chosen is 0!
# It HAS to be 0 since we are definitely sampling from one of its
# children which means that this node will be blocked by those samples.
# The difficult part now is to prove that the sampling the way we do it
# is really random.
if n > leaves/2:
# if n > leaves/2, at least n-leaves/2 have to stay on either
# side, otherweise the other one would have to sample more
# nodes than possible.
nMin = n-leaves/2
else:
nMin = 0
# The rest n-2*nMin is the amount of samples that are free to fall
# on either side
free = n-2*nMin
# Let have a look at an example, suppose we were to distribute 5
# samples among two children which have 4 leaves each.
# Each child has to get at least 1 sample, so the free samples are 3.
# There are 4 different ways to split the samples among the
# children (left, right):
# (1, 4), (2, 3), (3, 2), (4, 1)
# The amount of unique sample combinations per child are
# (7, 1), (11, 6), (6, 11), (1, 7)
# The amount of total unique samples per possible split are
# 7 , 66 , 66 , 7
# In case of the first and last split, all samples have a probability
# of 1/7, this was already proven above.
# Lets suppose we are good to go and the per sample probabilities for
# the other two cases are (1/11, 1/6) and (1/6, 1/11), this way the
# overall per sample probabilities for the splits would be:
# 1/7 , 1/66 , 1/66 , 1/7
# If we used uniform random to determine the split, all splits would be
# equally probable and therefore be multiplied with the same value (1/4)
# But this would mean that NOT every sample is equally probable!
# We need to know in advance how many sample combinations there will be
# for a given split in order to find out the probability to choose it.
# In fact, due to the restrictions, this becomes very nasty to
# determine. So instead of solving it analytically, I do it numerically
# with the method 'combs'. It gives me the amount of possible sample
# combinations for a certain amount of samples and a given tree depth.
# It will return 146 for this node and 7 for the outer and 66 for the
# inner splits.
# What we now do is, we take a number from [0, 146).
# if it is smaller than 7, we sample from the first split,
# if it is smaller than 7+66, we sample from the second split,
# ...
# This way we get the probabilities we need.
r = random.randint(0, self.combs(n)-1)
p = 0
for addLeft in xrange(0, free+1):
nLeft = nMin + addLeft
nRight = n - nLeft
p += (self.left.combs(nLeft) * self.right.combs(nRight))
if r < p:
return self.left.sample(nLeft) + self.right.sample(nRight)
assert False, ('Something really strange happend, p did not sum up ' +
'to combs or r was too big')
def main():
"""
Do a microbenchmark.
"""
import timeit
i = 1
main.t = Tree()
template = ' {:>2} {:>5} {:>4} {:<5}'
print(template.format('i', 'depth', 'n', 'time (ms)'))
N = 100
for depth in [4, 8, 15, 16, 17, 18]:
for n in [10, 50, 100, 150]:
if n > 2**depth:
time = '--'
else:
time = timeit.timeit(
'main.t.sample({}, depth={})'.format(n, depth), setup=
'from __main__ import main', number=N)*1000./N
print(template.format(i, depth, n, time))
i += 1
if __name__ == "__main__":
main()
Результат теста:
i depth n time (ms)
1 4 10 0.182511806488
2 4 50 --
3 4 100 --
4 4 150 --
5 8 10 0.397620201111
6 8 50 1.66054964066
7 8 100 2.90236949921
8 8 150 3.48146915436
9 15 10 0.804011821747
10 15 50 3.7428188324
11 15 100 7.34910964966
12 15 150 10.8230614662
13 16 10 0.804491043091
14 16 50 3.66818904877
15 16 100 7.09567070007
16 16 150 10.404779911
17 17 10 0.865840911865
18 17 50 3.9999294281
19 17 100 7.70257949829
20 17 150 11.3758206367
21 18 10 0.915451049805
22 18 50 4.22935962677
23 18 100 8.22361946106
24 18 150 12.2081303596
10 образцов размера 10 из дерева с глубиной 10:
['1111010111', '1110111010', '1010111010', '011110010', '0111100001', '011101110', '01110010', '01001111', '0001000100', '000001010']
['110', '0110101110', '0110001100', '0011110', '0001111011', '0001100010', '0001100001', '0001100000', '0000011010', '0000001111']
['11010000', '1011111101', '1010001101', '1001110001', '1001100110', '10001110', '011111110', '011001100', '0101110000', '001110101']
['11111101', '110111', '110110111', '1101010101', '1101001011', '1001001100', '100100010', '0100001010', '0100000111', '0010010110']
['111101000', '1110111101', '1101101', '1101000000', '1011110001', '0111111101', '01101011', '011010011', '01100010', '0101100110']
['1111110001', '11000110', '1100010100', '101010000', '1010010001', '100011001', '100000110', '0100001111', '001101100', '0001101101']
['111110010', '1110100', '1101000011', '101101', '101000101', '1000001010', '0111100', '0101010011', '0101000110', '000100111']
['111100111', '1110001110', '1100111111', '1100110010', '11000110', '1011111111', '0111111', '0110000100', '0100011', '0010110111']
['1101011010', '1011111', '1011100100', '1010000010', '10010', '1000010100', '0111011111', '01010101', '001101', '000101100']
['111111110', '111101001', '1110111011', '111011011', '1001011101', '1000010100', '0111010101', '010100110', '0100001101', '0010000000']
Ответ 6
Если вы не хотите генерировать набор всех возможных кортежей, а затем произвольно отбирать (что, как вы заметили, может быть неосуществимо для больших размеров ввода), другой вариант состоит в том, чтобы нарисовать один образец с целым программированием. В принципе, вы можете назначить каждому элементу в pool
случайное значение, а затем выбрать допустимый кортеж с наибольшей суммой значений. Это должно дать каждому кортежу равную вероятность выбора, потому что они имеют одинаковый размер, и их значения были выбраны случайным образом. Ограничения модели гарантируют, что ни одна из запрещенных пар кортежей не выбрана и что выбрано правильное количество элементов.
Здесь решение с пакетом lpSolve
:
library(lpSolve)
sample.lp <- function(pool, max_len) {
pool <- sort(pool)
pml <- max(nchar(pool))
runs <- c(rev(cumsum(2^(seq(pml-1)))), 0)
banned.from <- rep(seq(pool), runs[nchar(pool)])
banned.to <- banned.from + unlist(lapply(runs[nchar(pool)], seq_len))
banned.constr <- matrix(0, nrow=length(banned.from), ncol=length(pool))
banned.constr[cbind(seq(banned.from), banned.from)] <- 1
banned.constr[cbind(seq(banned.to), banned.to)] <- 1
mod <- lp(direction="max",
objective.in=runif(length(pool)),
const.mat=rbind(banned.constr, rep(1, length(pool))),
const.dir=c(rep("<=", length(banned.from)), "=="),
const.rhs=c(rep(1, length(banned.from)), max_len),
all.bin=TRUE)
pool[which(mod$solution == 1)]
}
set.seed(144)
pool <- unlist(lapply(seq_len(4), function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
sample.lp(pool, 4)
# [1] "0011" "010" "1000" "1100"
sample.lp(pool, 8)
# [1] "0000" "0100" "0110" "1001" "1010" "1100" "1101" "1110"
Это похоже на масштабные пулы. Например, для получения образца длиной 20 из пула размером 510 требуется менее 2 секунд:
pool <- unlist(lapply(seq_len(8), function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
length(pool)
# [1] 510
system.time(sample.lp(pool, 20))
# user system elapsed
# 0.232 0.008 0.239
Если вам действительно нужны действительно огромные размеры проблем, вы можете перейти от не-открытых исходников, которые отправляются с lpSolve
в коммерческий решатель, такой как gurobi или cplex (не бесплатный вообще, но бесплатный для академического использования).
Ответ 7
Один из подходов состоит в том, чтобы просто генерировать все возможные кортежи соответствующего размера с помощью итеративного подхода:
- Создайте все кортежи размером 1 (все элементы в
pool
)
- Возьмите крестик с элементами в
pool
- Удалить один кортеж одним и тем же элементом
pool
более одного раза
- Удалите любой точный дубликат другого кортежа
- Удалите любой кортеж с помощью пары, которая не может использоваться вместе
- Промыть и повторить, пока вы не получите соответствующий размер кортежа
Это выполняется для заданных размеров (pool
длины 30, max_len
4):
get.template <- function(pool, max_len) {
banned <- which(outer(paste0('^', pool), pool, Vectorize(grepl)), arr.ind=T)
banned <- banned[banned[,1] != banned[,2],]
banned <- paste(banned[,1], banned[,2])
vals <- matrix(seq(length(pool)))
for (k in 2:max_len) {
vals <- cbind(vals[rep(1:nrow(vals), each=length(pool)),],
rep(1:length(pool), nrow(vals)))
# Can't sample same value more than once
vals <- vals[apply(vals, 1, function(x) length(unique(x)) == length(x)),]
# Sort rows to ensure unique only
vals <- t(apply(vals, 1, sort))
vals <- unique(vals)
# Can't have banned pair
combos <- combn(ncol(vals), 2)
for (k in seq(ncol(combos))) {
c1 <- combos[1,k]
c2 <- combos[2,k]
vals <- vals[!paste(vals[,c1], vals[,c2]) %in% banned,]
}
}
return(matrix(pool[vals], nrow=nrow(vals)))
}
max_len <- 4
pool <- unlist(lapply(seq_len(max_len), function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x)))))
system.time(template <- get.template(pool, 4))
# user system elapsed
# 4.549 0.050 4.614
Теперь вы можете отбирать столько раз, сколько нужно, из строк template
(что будет очень быстро), и это будет то же самое, что и случайная выборка из указанного пространства.
Ответ 8
Введение
Я нашел этот вопрос настолько интересным, что я чувствовал себя вынужденным подумать об этом и в конечном итоге предоставить свой собственный ответ. Поскольку алгоритм, к которому я пришел, не сразу следует из описания проблемы, я начну с объяснения того, как я пришел к этому решению, а затем предоставил примерную реализацию на С++ (я никогда не писал R).
Разработка решения
На первый взгляд
Чтение описания проблемы изначально было запутанным, но как только я увидел редактирование с изображениями деревьев, я сразу понял проблему, и моя интуиция предположила, что бинарное дерево также является решением: построить дерево (сборник деревьев размером 1) и разбить дерево на коллекцию меньших деревьев после устранения ветвей и предков при выборе.
Хотя это изначально казалось хорошим, процесс отбора и обслуживания коллекции был бы болью. Тем не менее, дерево казалось, что оно должно играть большую роль в любом решении.
Версия 1
Не разбивайте дерево. Вместо этого на каждом node имеется полезная нагрузка булевых данных, указывающая, была ли она уже устранена. Это оставляет только одно дерево, которое поддерживает форму.
Обратите внимание, что это не просто бинарное дерево, это фактически полное двоичное дерево глубины max_len-1
.
Версия 2
Полное двоичное дерево может быть очень хорошо представлено в виде массива. Типичное представление массива использовало поиск по дереву в ширину, имеет следующие свойства:
Let x be the array index.
x = 0 is the root of the entire tree
left_child(x) = 2x + 1
right_child(x) = 2x + 2
parent(x) = floor((n-1)/2)
На рисунке ниже каждый node помечен индексом массива:
![breadth-first]()
В качестве массива это занимает меньше памяти (не более указателей), используемая память является непрерывной (хорошо подходит для кэширования) и может полностью переходить в стек, а не в кучу (при условии, что ваш язык дает вам выбор). Конечно, здесь применяются некоторые условия, в частности размер массива. Я вернусь к этому позже.
Как и в редакции 1, данные, хранящиеся в массиве, будут иметь значение boolean: true для доступных, false для недоступных. Поскольку корень node на самом деле не является допустимым выбором, индекс 0 должен быть инициализирован равным false. Проблема в том, как сделать выбор по-прежнему остается:
По мере устранения указаний, тривиально отслеживать, сколько из них устранено, и, следовательно, сколько осталось. Выберите случайное число в этом диапазоне, затем пройдите массив до тех пор, пока не увидите, что многие указатели установлены в true (включая текущий индекс). Прибывший индекс - это выбор.
Выберите до тех пор, пока не будут выбраны n указателей, или нет ничего, что можно было бы выбрать.
Это полный алгоритм, который будет работать, однако есть возможности для улучшения процесса выбора, а также существует практическая проблема с размером, который не был рассмотрен: размер массива в котором был бы равен O (2 ^ n). По мере увеличения размера n сначала извлекается преимущество кэширования, затем данные начинают выгружаться на диск, и в какой-то момент становится невозможным сохранить его.
Версия 3
Я решил сначала решить более легкую проблему: улучшить процесс выбора.
Сканирование массива слева направо напрасно. Возможно, более эффективно отслеживать диапазоны, которые были устранены, а не проверять и находить несколько фальши в строке. Однако наше древовидное представление не является идеальным для этого, так как немногие из узлов, которые будут устранены каждый раунд, будут смежными в массиве.
Переупорядочивая, как массив сопоставляется с деревом, можно, однако, лучше использовать это. В частности, мы будем использовать предварительный поиск по глубине, а не поиск по ширине. Для этого дерево необходимо фиксировать по размеру, что и относится к этой проблеме. Это также менее очевидно, так как математически математически связаны индексы дочерних и родительских узлов.
![depth-first]()
Используя эту компоновку, каждый выбор, который не является листом, гарантированно устраняет смежный диапазон: его поддерево.
Версия 4
Отслеживая удаленные диапазоны, больше нет необходимости в истинных/ложных данных и, следовательно, нет необходимости в массиве или дереве вообще.
При каждом случайном розыгрыше исключенные диапазоны можно использовать для быстрого поиска node для выбора. Все предки и все поддеревья удаляются и могут быть представлены как диапазоны, которые легко могут быть объединены с другими.
Последняя задача - превратить выбранные узлы в строковое представление, которое требуется OP. Это довольно просто, так как это двоичное дерево по-прежнему поддерживает строгий порядок: переходя от корня, все элементы >= правый ребенок справа, а остальные - слева. Таким образом, поиск дерева даст как список предков, так и двоичную строку, добавив "0" при прохождении влево; или "1" при правильном перемещении.
Пример реализации
#include <stdint.h>
#include <algorithm>
#include <cmath>
#include <list>
#include <deque>
#include <ctime>
#include <cstdlib>
#include <iostream>
/*
* A range of values of the form (a, b), where a <= b, and is inclusive.
* Ex (1,1) is the range from 1 to 1 (ie: just 1)
*/
class Range
{
private:
friend bool operator< (const Range& lhs, const Range& rhs);
friend std::ostream& operator<<(std::ostream& os, const Range& obj);
int64_t m_start;
int64_t m_end;
public:
Range(int64_t start, int64_t end) : m_start(start), m_end(end) {}
int64_t getStart() const { return m_start; }
int64_t getEnd() const { return m_end; }
int64_t size() const { return m_end - m_start + 1; }
bool canMerge(const Range& other) const {
return !((other.m_start > m_end + 1) || (m_start > other.m_end + 1));
}
int64_t merge(const Range& other) {
int64_t change = 0;
if (m_start > other.m_start) {
change += m_start - other.m_start;
m_start = other.m_start;
}
if (other.m_end > m_end) {
change += other.m_end - m_end;
m_end = other.m_end;
}
return change;
}
};
inline bool operator< (const Range& lhs, const Range& rhs){return lhs.m_start < rhs.m_start;}
std::ostream& operator<<(std::ostream& os, const Range& obj) {
os << '(' << obj.m_start << ',' << obj.m_end << ')';
return os;
}
/*
* Stuct to allow returning of multiple values
*/
struct NodeInfo {
int64_t subTreeSize;
int64_t depth;
std::list<int64_t> ancestors;
std::string representation;
};
/*
* Collection of functions representing a complete binary tree
* as an array created using pre-order depth-first search,
* with 0 as the root.
* Depth of the root is defined as 0.
*/
class Tree
{
private:
int64_t m_depth;
public:
Tree(int64_t depth) : m_depth(depth) {}
int64_t size() const {
return (int64_t(1) << (m_depth+1))-1;
}
int64_t getDepthOf(int64_t node) const{
if (node == 0) { return 0; }
int64_t searchDepth = m_depth;
int64_t currentDepth = 1;
while (true) {
int64_t rightChild = int64_t(1) << searchDepth;
if (node == 1 || node == rightChild) {
break;
} else if (node > rightChild) {
node -= rightChild;
} else {
node -= 1;
}
currentDepth += 1;
searchDepth -= 1;
}
return currentDepth;
}
int64_t getSubtreeSizeOf(int64_t node, int64_t nodeDepth = -1) const {
if (node == 0) {
return size();
}
if (nodeDepth == -1) {
nodeDepth = getDepthOf(node);
}
return (int64_t(1) << (m_depth + 1 - nodeDepth)) - 1;
}
int64_t getLeftChildOf(int64_t node, int64_t nodeDepth = -1) const {
if (nodeDepth == -1) {
nodeDepth = getDepthOf(node);
}
if (nodeDepth == m_depth) { return -1; }
return node + 1;
}
int64_t getRightChildOf(int64_t node, int64_t nodeDepth = -1) const {
if (nodeDepth == -1) {
nodeDepth = getDepthOf(node);
}
if (nodeDepth == m_depth) { return -1; }
return node + 1 + ((getSubtreeSizeOf(node, nodeDepth) - 1) / 2);
}
NodeInfo getNodeInfo(int64_t node) const {
NodeInfo info;
int64_t depth = 0;
int64_t currentNode = 0;
while (currentNode != node) {
if (currentNode != 0) {
info.ancestors.push_back(currentNode);
}
int64_t rightChild = getRightChildOf(currentNode, depth);
if (rightChild == -1) {
break;
} else if (node >= rightChild) {
info.representation += '1';
currentNode = rightChild;
} else {
info.representation += '0';
currentNode = getLeftChildOf(currentNode, depth);
}
depth++;
}
info.depth = depth;
info.subTreeSize = getSubtreeSizeOf(node, depth);
return info;
}
};
// random selection amongst remaining allowed nodes
int64_t selectNode(const std::deque<Range>& eliminationList, int64_t poolSize, std::mt19937_64& randomGenerator)
{
std::uniform_int_distribution<> randomDistribution(1, poolSize);
int64_t selection = randomDistribution(randomGenerator);
for (auto const& range : eliminationList) {
if (selection >= range.getStart()) { selection += range.size(); }
else { break; }
}
return selection;
}
// determin how many nodes have been elimintated
int64_t countEliminated(const std::deque<Range>& eliminationList)
{
int64_t count = 0;
for (auto const& range : eliminationList) {
count += range.size();
}
return count;
}
// merge all the elimination ranges to listA, and return the number of new elimintations
int64_t mergeEliminations(std::deque<Range>& listA, std::deque<Range>& listB) {
if(listB.empty()) { return 0; }
if(listA.empty()) {
listA.swap(listB);
return countEliminated(listA);
}
int64_t newEliminations = 0;
int64_t x = 0;
auto listA_iter = listA.begin();
auto listB_iter = listB.begin();
while (listB_iter != listB.end()) {
if (listA_iter == listA.end()) {
listA_iter = listA.insert(listA_iter, *listB_iter);
x = listB_iter->size();
assert(x >= 0);
newEliminations += x;
++listB_iter;
} else if (listA_iter->canMerge(*listB_iter)) {
x = listA_iter->merge(*listB_iter);
assert(x >= 0);
newEliminations += x;
++listB_iter;
} else if (*listB_iter < *listA_iter) {
listA_iter = listA.insert(listA_iter, *listB_iter) + 1;
x = listB_iter->size();
assert(x >= 0);
newEliminations += x;
++listB_iter;
} else if ((listA_iter+1) != listA.end() && listA_iter->canMerge(*(listA_iter+1))) {
listA_iter->merge(*(listA_iter+1));
listA_iter = listA.erase(listA_iter+1);
} else {
++listA_iter;
}
}
while (listA_iter != listA.end()) {
if ((listA_iter+1) != listA.end() && listA_iter->canMerge(*(listA_iter+1))) {
listA_iter->merge(*(listA_iter+1));
listA_iter = listA.erase(listA_iter+1);
} else {
++listA_iter;
}
}
return newEliminations;
}
int main (int argc, char** argv)
{
std::random_device rd;
std::mt19937_64 randomGenerator(rd());
int64_t max_len = std::stoll(argv[1]);
int64_t num_samples = std::stoll(argv[2]);
int64_t samplesRemaining = num_samples;
Tree tree(max_len);
int64_t poolSize = tree.size() - 1;
std::deque<Range> eliminationList;
std::deque<Range> eliminated;
std::list<std::string> foundList;
while (samplesRemaining > 0 && poolSize > 0) {
// find a valid node
int64_t selectedNode = selectNode(eliminationList, poolSize, randomGenerator);
NodeInfo info = tree.getNodeInfo(selectedNode);
foundList.push_back(info.representation);
samplesRemaining--;
// determine which nodes this choice eliminates
eliminated.clear();
for( auto const& ancestor : info.ancestors) {
Range r(ancestor, ancestor);
if(eliminated.empty() || !eliminated.back().canMerge(r)) {
eliminated.push_back(r);
} else {
eliminated.back().merge(r);
}
}
Range r(selectedNode, selectedNode + info.subTreeSize - 1);
if(eliminated.empty() || !eliminated.back().canMerge(r)) {
eliminated.push_back(r);
} else {
eliminated.back().merge(r);
}
// add the eliminated nodes to the existing list
poolSize -= mergeEliminations(eliminationList, eliminated);
}
// Print some stats
// std::cout << "tree: " << tree.size() << " samplesRemaining: "
// << samplesRemaining << " poolSize: "
// << poolSize << " samples: " << foundList.size()
// << " eliminated: "
// << countEliminated(eliminationList) << std::endl;
// Print list of binary strings
// std::cout << "list:";
// for (auto const& s : foundList) {
// std::cout << " " << s;
// }
// std::cout << std::endl;
}
Дополнительные мысли
Этот алгоритм будет очень хорошо масштабироваться для max_len. Масштабирование с n не очень хорошо, но, основываясь на моем собственном профилировании, похоже, это лучше, чем другие решения.
Этот алгоритм может быть изменен без особых усилий, чтобы строки, содержащие больше, чем просто "0" и "1". Более возможные буквы в словах увеличивают разветвление дерева и устраняют более широкий диапазон с каждым выбором - все еще со всеми узлами в каждом поддереве, остающемся смежным.