Domanda Esempio casuale di vettore di caratteri, senza elementi che si prefissano l'un l'altro


Consideri un vettore di caratteri, pool, i cui elementi sono numeri binari (azzerati) fino a max_len cifre.

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"

Mi piacerebbe provare n di questi elementi, con il vincolo che nessuno degli elementi campionati è prefissi di uno qualsiasi degli altri elementi campionati (cioè se campioniamo 1101, vietiamo 1, 11, e 110, mentre se proviamo 1, vietiamo quegli elementi che iniziano con 1, ad esempio 10, 11, 100, eccetera.).

Quello che segue è il mio tentativo di usare while, ma ovviamente questo è lento quando n è grande (o quando si avvicina 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"

Questo può essere leggermente migliorato rimuovendo inizialmente da pool quegli elementi la cui inclusione significherebbe che non sono rimasti elementi sufficienti pool per prendere un campione totale di dimensioni n. Ad esempio, quando max_len = 4 e n > 9, possiamo rimuovere immediatamente 0 e 1 a partire dal pool, poiché includendo entrambi, il campione massimo sarebbe 9 (o 0 e gli otto elementi di 4 caratteri che iniziano con 1, o 1 e gli otto elementi di 4 caratteri che iniziano con 0).

Sulla base di questa logica, potremmo quindi omettere gli elementi pool, prima di prendere il campione iniziale, con qualcosa di simile:

pool <- pool[
  nchar(pool) > tail(which(n > (2^max_len - rev(2^(0:max_len))[-1] + 1)), 1)]

Qualcuno può pensare ad un approccio migliore? Mi sento come se stessi trascurando qualcosa di molto più semplice.


MODIFICARE

Per chiarire le mie intenzioni, ritroverò il pool come un insieme di rami, in cui i nodi e i suggerimenti sono i nodi (elementi di pool). Supponiamo che il nodo giallo nella figura seguente (vale a dire 010) sia stato disegnato. Ora, l'intero "ramo" rosso, costituito dai nodi 0, 01 e 010, viene rimosso dal pool. Questo è quello che intendevo proibendo il campionamento dei nodi che già "prefissano" i nodi nel nostro campione (così come quelli già prefissata dai nodi nel nostro esempio).

enter image description here

Se il nodo campionato era a metà strada lungo un ramo, ad esempio 01 nella figura seguente, tutti i nodi rossi (0, 01, 010 e 011) non sono consentiti, poiché 0 prefissi 01 e 01 prefissi entrambi 010 e 011.

enter image description here

Non intendo campionare o 1 o 0 a ogni incrocio (cioè camminando lungo i rami che lanciano monete nelle forcelle) - è bene avere entrambi nel campione, purché: (1) genitori (o nonni, ecc.) O figli (nipoti, ecc.) del nodo non sono già campionati; e (2) al momento del campionamento del nodo ci saranno dei nodi sufficienti per ottenere il campione desiderato di dimensione n.

Nella seconda figura sopra, se 010 è stata la prima scelta, tutti i nodi a nodi neri sono ancora (attualmente) validi, assumendo n <= 4. Ad esempio, se n==4 e abbiamo campionato il nodo 1 successivo (e quindi le nostre scelte ora includevano 01 e 1), successivamente non avremmo consentito il nodo 00 (dovuto alla regola 2 sopra), ma potremmo comunque selezionare 000 e 001, dandoci il nostro campione di 4 elementi. Se n==5, d'altra parte, il nodo 1 non sarebbe consentito in questa fase.


44
2018-06-10 23:42


origine


risposte:


Intro

Questa è una variazione numerica dell'algoritmo di stringa che abbiamo implementato nell'altra risposta. È più veloce e non richiede né la creazione né l'ordinamento del pool.

Profilo dell'algoritmo

Possiamo usare numeri interi per rappresentare le stringhe binarie, il che semplifica enormemente il problema della generazione del pool e dell'eliminazione sequenziale dei valori. Ad esempio, con max_len==3, possiamo prendere il numero 1-- (dove - rappresenta padding) per significare 4 in decimale. Inoltre, possiamo determinare che i numeri che richiedono l'eliminazione se prendiamo questo numero sono quelli in mezzo 4 e 4 + 2 ^ x - 1. Qui x è il numero di elementi di riempimento (2 in questo caso), quindi i numeri da eliminare sono tra 4 e 4 + 2 ^ 2 - 1 (o tra 4 e 7, espresso come 100, 110, e 111).

Per far combaciare esattamente il tuo problema, abbiamo bisogno di un po 'di rughe dato che tratti i numeri che sarebbero potenzialmente gli stessi in binario come distinti per alcune parti del tuo algoritmo. Per esempio, 100, 10-, e 1-- sono tutti lo stesso numero, ma devono essere trattati in modo diverso nel tuo schema. In un max_len==3 mondo, abbiamo 8 numeri possibili, ma 14 possibili rappresentazioni:

0 - 000: 0--, 00-
1 - 001:
2 - 010: 01-
3 - 011:
4 - 100: 1--, 10-
5 - 101:
6 - 110: 11-
7 - 111:

Quindi 0 e 4 hanno tre possibili codifiche, 2 e 6 ne hanno due, e tutte le altre solo una. Abbiamo bisogno di generare un pool intero che rappresenti la maggiore probabilità di selezione per i numeri con rappresentazione multipla, oltre al meccanismo per tenere traccia di quanti spazi il numero include. Possiamo farlo aggiungendo alcuni bit alla fine del numero per indicare le ponderazioni che vogliamo. Quindi i nostri numeri diventano (usiamo due bit qui):

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

Alcune proprietà utili:

  • enc.bits rappresenta quanti bit abbiamo bisogno per la codifica (due in questo caso)
  • int.enc %% enc.bits ci dice quanti dei numeri sono esplicitamente specificati
  • int.enc %/% enc.bits ritorna int
  • int * 2 ^ enc.bits + explicitly.specified ritorna int.enc

Nota che explicitly.specified qui è tra 0 e max_len - 1 nella nostra implementazione poiché c'è sempre almeno una cifra specificata. Ora abbiamo una codifica che rappresenta completamente la struttura dei dati utilizzando solo gli interi. Possiamo campionare da interi e riprodurre i risultati desiderati, con pesi corretti, ecc. Una limitazione di questo approccio è che usiamo numeri interi a 32 bit in R, e dobbiamo riservare alcuni bit per la codifica, quindi ci limitiamo a pool di max_len==25 o così. Potresti andare più grande se usassi interi specificati da punti mobili a precisione doppia, ma non lo abbiamo fatto qui.

Evitare le selezioni duplicate

Esistono due modi per garantire che non selezioniamo lo stesso valore due volte

  1. Tieni traccia di quali valori rimangono disponibili per il prelievo e campionali casualmente da quelli
  2. Campiona casualmente da tutti i possibili valori, quindi controlla se il valore è stato selezionato in precedenza e, se lo ha, prova nuovamente

Mentre la prima opzione sembra la più pulita, in realtà è molto dispendiosa dal punto di vista computazionale. Richiede sia una scansione vettoriale di tutti i possibili valori per ciascun prelievo per pre-squalificare i valori selezionati, sia creare un vettore di restringimento contenente valori non squalificati. L'opzione di restringimento è più efficiente della scansione vettoriale se il vettore viene fatto restringere per riferimento tramite codice C, ma anche allora richiede ripetute traduzioni di porzioni potenzialmente grandi del vettore, e richiede C.

Qui usiamo il metodo # 2. Questo ci permette di mescolare casualmente l'universo dei possibili valori una volta, e quindi selezionare ciascun valore in sequenza, controllare che non sia stato squalificato, e se lo ha, sceglierne un altro, ecc. Funziona perché è banale controllare se un il valore è stato scelto come risultato della nostra codifica del valore; possiamo dedurre la posizione di un valore in una tabella ordinata in base al solo valore. Quindi registriamo lo stato di ogni valore in una tabella ordinata e possiamo aggiornare o cercare quello stato tramite accesso diretto all'indice (nessuna scansione richiesta).

Esempi

L'implementazione di questo algoritmo nella base R è disponibile come un succo. Questa particolare implementazione tira solo i progetti completi. Ecco un esempio di 10 estrazioni di 8 elementi di a max_len==4 piscina:

# 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"

Inizialmente avevamo anche due implementazioni che si basavano sul metodo # 1 per evitare duplicati, uno in base R e uno in C, ma anche la versione C non è veloce come la nuova versione di base R quando n è grande. Queste funzioni implementano la capacità di disegnare disegni incompleti, quindi li forniamo qui per riferimento:

Benchmark comparativi

Ecco una serie di benchmark che confrontano molte delle funzioni visualizzate in questa Q / A. Tempi in millisecondi. Il brodie.b la versione è quella descritta in questa risposta. brodie è l'implementazione originale, brodie.C è l'originale con alcuni C. Tutti questi fanno rispettare il requisito per i campioni completi. brodie.str è la versione basata su stringa nell'altra risposta.

   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

Questo scala relativamente bene per le piscine più grandi

system.time(sample0110b(18, 100000))
   user  system elapsed 
  8.441   0.079   8.527 

Note di riferimento:

  • frank e brodie (meno brodie.str) non richiedono alcuna pre-generazione di pool, il che influenzerà i confronti (vedi sotto)
  • Josilber è la versione LP
  • jbaum è l'esempio OP
  • tensibai è leggermente modificato per uscire invece di fallire se il pool è vuoto
  • non è impostato per eseguire python, quindi non è possibile confrontare / account per il buffering
  • - rappresentano o opzioni non fattibili o troppo lente nel tempo ragionevolmente

I tempi non includono il disegno delle piscine (0.8, 2.5, 401 millisecondi rispettivamente per dimensione 4, 8, e 16), che è necessario per il jbaum, josilber, e brodie.str corre, o ordinandoli (0.1, 2.7, 3700 millisecondi per dimensione 4, 8, e 16), che è necessario per brodie.str oltre al sorteggio. Se si desidera includerli o meno dipende da quante volte si esegue la funzione per un determinato pool. Inoltre, ci sono quasi certamente modi migliori per generare / ordinare i pool.

Questi sono il tempo medio di tre esecuzioni con microbenchmark. Il codice è disponibile come sintesi, sebbene si noti che dovrete caricare il file sample0110b, sample0110, sample01101, e sample01 funzioni in anticipo.


19
2018-06-11 12:40



Ho trovato il problema interessante, quindi ho provato e ottenere questo con competenze veramente basse in R (quindi questo può essere migliorato):

Versione modificata più recente, grazie a @Franck consigli:

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's 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's 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)

L'idea principale è usare il confronto maschera di bitper sapere se un campione è padre dell'altro (parte bit comune), se si sopprime questo elemento dal pool.

Occorrono 1,4 secondi per prelevare 100 campioni da un pool di 16 sulla mia macchina.


15
2018-06-17 08:44



È possibile ordinare il pool per decidere quali elementi da squalificare. Ad esempio, guardando un pool ordinato di tre elementi:

 [1] "0"   "00"  "000" "001" "01"  "010" "011" "1"   "10"  "100" "101" "11" 
[13] "110" "111"

Posso dire che posso squalificare tutto ciò che segue il mio elemento selezionato che ha più caratteri del mio oggetto fino al primo oggetto che ha lo stesso numero o meno caratteri. Ad esempio, se seleziono "01", posso immediatamente vedere che i due elementi successivi ("010", "011") devono essere rimossi, ma non quello successivo perché "1" ha meno caratteri. Rimuovere lo "0" in seguito è facile. Ecco una implementazione:

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  
}

E una funzione per creare pool ordinati (esattamente come il tuo codice, ma restituisce ordinati):

make_pool <- function(size)
  sort(
    unlist(
      lapply(
        seq_len(size), 
        function(x) do.call(paste0, expand.grid(rep(list(c('0', '1')), x))) 
  ) ) )

Quindi, usando il max_len 3 pool (utile per ispezionare visivamente le cose comportarsi come previsto):

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"    

Notate nell'ultimo caso che otteniamo tutte le combinazioni binarie a 3 cifre (2 ^ 3) perché per caso abbiamo continuato a campionare da quelle a 3 cifre. Inoltre, con solo una piscina di dimensioni 3 ci sono molti campionamenti che impediscono un totale di 8 pareggi; puoi affrontare questo problema con il tuo suggerimento di eliminare le combinazioni che impediscono la piena estrazione dal pool.

Questo è abbastanza veloce. Guardando il max_len==9 esempio che ha impiegato 2 secondi con la soluzione alternativa:

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    

Circa mezzo millisecondo. Puoi ragionevolmente provare anche piscine abbastanza grandi:

pool16 <- make_pool(16)  # 131K entries
system.time(sample01(pool16, 100))
#  user  system elapsed 
# 3.407   0.146   3.552 

Questo non è molto veloce, ma stiamo parlando di un pool con 130K elementi. C'è anche il potenziale per l'ottimizzazione aggiuntiva.

Nota che il passaggio di ordinamento diventa relativamente lento per i grandi pool, ma non lo sto contando perché hai sempre bisogno di farlo una volta, e potresti probabilmente trovare un algoritmo ragionevole per produrre il pool pre ordinato.

C'è anche la possibilità di un approccio molto più veloce all'approccio binario che ho esplorato in una risposta ora cancellata, ma ciò richiede un bel po 'di lavoro in più per legare esattamente quello che stai cercando.


13
2018-06-11 02:39



È in python e non in r ma jbaums ha detto che andrebbe bene.

Quindi, ecco il mio contributo, vedere i commenti nella fonte per la spiegazione delle parti cruciali.
Sto ancora lavorando alla soluzione analitica per determinare la quantità di combinazioni possibili c per un albero di profondità t e S campioni, quindi posso migliorare la funzione combs. Forse qualcuno lo ha? Questo è davvero il collo di bottiglia ora.

Campionare 100 nodi da un albero con profondità 16 richiede circa 8ms sul mio portatile. Non è la prima volta, ma diventa più veloce fino a un certo punto più si campiona perché combBuffer viene riempito.

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's 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's 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's depth. (maxdepth-level)
        """
        return self.maxdepth-self.level

    @property
    def id(self):
        """
        This tree's id, string of '0's and '1's 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's 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's ids is
        a prefix for another one's.

        For an example see Tree's 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's 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's 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()

L'output di riferimento:

  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 campioni di dimensione 10 da un albero con profondità 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']

13
2018-06-14 15:06



Mappatura degli ID alle stringhe. Puoi mappare i numeri ai tuoi vettori 0/1, come menzionato da @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      
} 

Trovare gli ID da abbandonare. Cadiamo sequenzialmente ids dal pool di campionamento:

 # 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)
}

Gli indici sono raggruppati per g, il numero di caratteri, nchar(get_01str(id)). Perché gli indici sono ordinati per g, g=findInterval(id,cuts) è un percorso più veloce.

Un indice in gruppo g, 1 < g < max_len ha un indice di dimensioni "bambino" g-1 e due indici di dimensioni principali g+1. Per ogni nodo figlio, prendiamo il suo nodo figlio fino a quando non colpiamo g==1; e per ogni nodo genitore, prendiamo la loro coppia di nodi genitore fino a quando non colpiamo g==max_len.

La struttura dell'albero è la più semplice in termini di identificatore all'interno del gruppo, gid. gid mappe per i due genitori, gid e gid+2^g; e invertendo questa mappatura si trova il bambino.

campionatura

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
}

Il oklens parte integra l'idea dell'OP per omettere le stringhe che sono garantite per rendere impossibile il campionamento. Tuttavia, anche facendo questo, potremmo seguire un percorso di campionamento che ci lascia senza ulteriori opzioni. Prendendo l'esempio del PO di max_len=4 e n=10, sappiamo che dobbiamo lasciare 0 e 1 dalla considerazione, ma cosa succede se i nostri primi quattro progetti sono 00, 01, 11 e 10? Oh bene, immagino che siamo sfortunati. Questo è il motivo per cui dovresti effettivamente definire le probabilità di campionamento. (L'OP ha un'altra idea, per determinare, ad ogni passo, quali nodi porteranno ad uno stato impossibile, ma sembra un compito arduo).

Illustrazione

# 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"  ""     ""  

benchmark (più vecchio di quelli nella risposta di @ 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

Altre risposte Ci sono altre due risposte:

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)    

Ho omesso il metodo di enumerazione di @josilber perché ci è voluto troppo tempo; e @ il metodo numerico / indice di BrodieG perché non funzionava al momento, ma ora lo fa. Vedi la risposta aggiornata di @ BrodieG per ulteriori analisi comparative.

Velocità vs correttezza. Mentre le risposte di @josilber sono molto più lente (e per il metodo di enumerazione, apparentemente molto più ricco di memoria), sono garantite per disegnare un campione di dimensioni n al primo tentativo. Con il metodo di stringa di @ BrodieG o questa risposta, dovrai ricampionare ancora e ancora nella speranza di ottenere un risultato completo n. Con grande max_len, non dovrebbe essere un tale problema, immagino.

Questa risposta è migliore di bro_string perché non richiede la costruzione di pool in anticipo.


13
2018-06-11 16:32



Se non si desidera generare l'insieme di tutte le possibili tuple e quindi campionare casualmente (che come si nota potrebbe non essere fattibile per grandi dimensioni di input), un'altra opzione consiste nel disegnare un singolo campione con la programmazione di interi. Fondamentalmente puoi assegnare ogni elemento in pool un valore casuale e quindi selezionare la tupla fattibile con la più grande somma di valori. Questo dovrebbe dare ad ogni tupla una uguale probabilità di essere selezionata perché sono tutte della stessa dimensione e i loro valori sono stati selezionati casualmente. I vincoli del modello assicurano che nessuna delle coppie di tuple non consentite sia selezionata e che sia selezionato il numero corretto di elementi.

Ecco una soluzione con il lpSolve pacchetto:

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"

Questo sembra ridimensionarsi a piscine abbastanza grandi. Ad esempio, ci vuole poco più di 2 secondi per ottenere un campione di lunghezza 20 da un pool di dimensioni 510:

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 

Se avevi bisogno di risolvere problemi davvero enormi, allora potevi passare dai solutori non open source con cui navighi lpSolvead un risolutore commerciale come gurobi o cplex (non libero in generale, ma gratuito per uso accademico).


11
2018-06-11 02:00



Un approccio sarebbe semplicemente generare tutte le possibili tuple delle dimensioni appropriate con un approccio iterativo:

  1. Crea tutte le tuple di dimensione 1 (tutti gli elementi in pool)
  2. Prendi un prodotto incrociato con elementi in pool
  3. Rimuovi qualsiasi tupla usando lo stesso elemento di pool più di una volta
  4. Rimuovi qualsiasi duplicato esatto di un'altra tupla
  5. Rimuovi qualsiasi tupla con una coppia che non può essere utilizzata insieme
  6. Risciacquare e ripetere fino ad ottenere la giusta dimensione della tupla

Questo è eseguibile per le dimensioni indicate (pool di lunghezza 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 

Ora puoi campionare tutte le volte che vuoi dalle file di template (che sarà molto veloce), e questo sarà lo stesso campionamento casuale dallo spazio definito.


9
2018-06-11 01:40