Подтвердить что ты не робот

Случайная выборка символьного вектора без элементов, префиксных друг другу

Рассмотрим вектор символов, 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 будет отменено на этом этапе.

4b9b3361

Ответ 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". Более возможные буквы в словах увеличивают разветвление дерева и устраняют более широкий диапазон с каждым выбором - все еще со всеми узлами в каждом поддереве, остающемся смежным.