usar - leer cadena de caracteres en c
Muestra aleatoria de vector de caracteres, sin elementos prefijándose entre sí. (8)
Introducción
Encontré esta pregunta tan interesante que me sentí obligada a reflexionar y, en última instancia, a dar mi propia respuesta. Dado que el algoritmo al que llegué no sigue inmediatamente de la descripción del problema, comenzaré explicando cómo llegué a esta solución y luego proporcionaré una implementación de ejemplo en C ++ (nunca he escrito R).
Desarrollo de la Solución.
A primera vista
La lectura de la descripción del problema fue inicialmente confusa, pero tan pronto como vi la edición con las imágenes de los árboles, comprendí de inmediato el problema y mi intuición sugirió que un árbol binario también era una solución: construir un árbol tamaño 1), y rompa el árbol en una colección de árboles más pequeños después de eliminar las ramas y los antepasados al hacer una selección.
Aunque inicialmente esto parecía bueno, el proceso de selección y el mantenimiento de la colección serían difíciles. Aún así, el árbol parecía que debería jugar un papel importante en cualquier solución.
Revisión 1
No rompas el árbol. En su lugar, tenga una carga útil de datos booleanos en cada nodo que indique si ya se ha eliminado. Esto deja solo un árbol que mantiene la forma.
Sin embargo, tenga en cuenta que esto no es un árbol binario, en realidad es un árbol binario completo de profundidad max_len-1
.
Revisión 2
Un árbol binario completo se puede representar muy bien como una matriz. La representación de matriz típica utilizada en la búsqueda en primer lugar del árbol, tiene las siguientes propiedades:
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)
En la siguiente ilustración, cada nodo está etiquetado con su índice de matriz:
Como una matriz, ocupa menos memoria (no más punteros), la memoria utilizada es contigua (buena para el almacenamiento en caché) y puede ir completamente en la pila en lugar del montón (suponiendo que su idioma le dé una opción). Por supuesto, algunas condiciones se aplican aquí, en particular el tamaño de la matriz. Volveré a eso más tarde.
Al igual que en la Revisión 1, los datos almacenados en la matriz serán booleanos: verdadero para disponible, falso para no disponible. Dado que el nodo raíz no es realmente una opción válida, el índice 0 debe inicializarse en falso. El problema de cómo hacer una selección aún permanece:
A medida que se eliminan los indicios, es trivial hacer un seguimiento de cuántos se eliminan y, por lo tanto, cuántos quedan. Elija un número aleatorio en ese rango, luego recorra la matriz hasta que haya visto que muchos de los indicadores están establecidos en verdadero (se incluye el índice actual). El índice al que llega es entonces la selección a realizar. Seleccione hasta que se seleccionen n indicaciones, o no quede nada para seleccionar.
Este es un algoritmo completo que funcionará, sin embargo, hay espacio para mejorar en el proceso de selección, y también existe un problema práctico de tamaño que no se ha abordado: el tamaño del arreglo sería O (2 ^ n). A medida que n aumenta, primero desaparece el beneficio de almacenamiento en caché, luego los datos comienzan a paginarse en el disco y, en algún momento, resulta imposible almacenarlos.
Revisión 3
Decidí abordar el problema más fácil primero: mejorar el proceso de selección.
Escanear la matriz de izquierda a derecha es un desperdicio. Puede ser más eficiente rastrear los rangos que se han eliminado en lugar de verificar y encontrar varias falsificaciones seguidas. Sin embargo, nuestra representación de árbol no es ideal para esto, ya que pocos de los nodos que serían eliminados en cada ronda serían contiguos en la matriz.
Al reorganizar la forma en que la matriz se asigna al árbol, es posible, sin embargo, explotar mejor esto. En particular, utilicemos la búsqueda de primer orden de profundidad en lugar de la búsqueda de primer alcance. Para hacerlo, el árbol debe ser de tamaño fijo, como es el caso de este problema. También es menos obvio cómo se conectan matemáticamente los índices de los nodos hijo y padre.
Al utilizar este arreglo, se garantiza que cualquier elección que no sea una hoja elimine un rango contiguo: su subárbol.
Revisión 4
Al rastrear los rangos eliminados, ya no hay necesidad de los datos verdaderos / falsos, y por lo tanto, no es necesaria la matriz o el árbol. En cada sorteo aleatorio, los rangos eliminados se pueden usar para encontrar rápidamente el nodo a seleccionar. Todos los antepasados y todo el subárbol se eliminan, y pueden representarse como rangos que pueden combinarse fácilmente con los demás.
La tarea final es convertir los nodos seleccionados en la representación de cadena que el OP deseaba. Eso es bastante fácil, ya que este árbol binario aún mantiene un orden estricto: atravesar desde la raíz, todos los elementos> = el hijo derecho está a la derecha, y los otros están a la izquierda. Por lo tanto, buscar en el árbol proporcionará tanto la lista de ancestros como la cadena binaria al agregar ''0'' cuando se recorra hacia la izquierda; o un ''1'' cuando se cruza a la derecha.
Implementación de la muestra
#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;
}
Pensamientos adicionales
Este algoritmo se escalará extremadamente bien para max_len. La escala con n no es muy buena, pero según mi propio perfil parece funcionar mejor que las otras soluciones.
Este algoritmo podría modificarse con poco esfuerzo para permitir cadenas que contengan más que solo ''0'' y ''1''. Más letras posibles en las palabras aumentan el abanico de salida del árbol y eliminarían un rango más amplio con cada selección, aún con todos los nodos en cada subárbol que permanecen contiguos.
Considere un vector de caracteres, pool
, cuyos elementos son números binarios (con relleno cero) con hasta max_len
dígitos.
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"
Me gustaría muestrear n
de estos elementos, con la restricción de que ninguno de los elementos muestreados son prefixes de ninguno de los otros elementos muestreados (es decir, si muestreamos 1101
, prohibimos 1
, 11
y 110
, mientras que si muestreamos 1
, prohibimos aquellos elementos que comienzan con 1
, como 10
, 11
, 100
, etc.).
El siguiente es mi intento de usar while
, pero por supuesto esto es lento cuando n
es grande (o cuando se acerca a 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"
Esto se puede mejorar ligeramente eliminando inicialmente del pool
aquellos elementos cuya inclusión significaría que no hay elementos suficientes en el pool
para tomar una muestra total de tamaño n
. Por ejemplo, cuando max_len = 4
y n > 9
, podemos eliminar inmediatamente 0
y 1
de la pool
, ya que al incluir cualquiera de ellos, la muestra máxima sería 9 (ya sea 0
y los ocho elementos de 4 caracteres que comienzan con 1
, o 1
y Ocho elementos de 4 caracteres que comienzan con 0
).
Según esta lógica, podríamos omitir elementos de la pool
, antes de tomar la muestra inicial, con algo como:
pool <- pool[
nchar(pool) > tail(which(n > (2^max_len - rev(2^(0:max_len))[-1] + 1)), 1)]
¿Alguien puede pensar en un mejor enfoque? Siento que estoy pasando por alto algo mucho más simple.
EDITAR
Para aclarar mis intenciones, retrataré el conjunto como un conjunto de ramas, donde las uniones y las puntas son los nodos (elementos del pool
). Supongamos que se dibujó el nodo amarillo en la siguiente figura (es decir, 010). Ahora, toda la "rama" roja, que consta de los nodos 0, 01 y 010, se elimina de la agrupación. Esto es lo que quise decir al prohibir el muestreo de nodos que "prefijan" los nodos que ya están en nuestra muestra (así como los que ya están prefijados por los nodos en nuestra muestra).
Si el nodo muestreado estaba a mitad de camino a lo largo de una rama, como 01 en la siguiente figura, entonces todos los nodos rojos (0, 01, 010 y 011) no están permitidos, ya que 0 prefijos 01 y 01 prefijos tanto 010 como 011.
No quiero probar 1 o 0 en cada cruce (es decir, caminar a lo largo de las ramas lanzando monedas en los tenedores); está bien tener ambos en la muestra, siempre y cuando: (1) los padres (o los abuelos, etc.) .) o los hijos (nietos, etc.) del nodo no se han muestreado; y (2) al muestrear el nodo, habrá suficientes nodos restantes para lograr la muestra deseada de tamaño n
.
En la segunda figura anterior, si 010 fue la primera selección, todos los nodos en los nodos negros siguen siendo (actualmente) válidos, suponiendo que n <= 4
. Por ejemplo, si n==4
y muestreamos el nodo 1 a continuación (y por lo tanto, nuestras selecciones ahora incluían 01 y 1), posteriormente rechazaremos el nodo 00 (debido a la regla 2 anterior), pero aún podríamos seleccionar 000 y 001, lo que nos da nuestra Muestra de 4 elementos. Si n==5
, por otro lado, el nodo 1 se rechazaría en esta etapa.
Introducción
Esta es una variación numérica del algoritmo de cadena que implementamos en la otra respuesta. Es más rápido y no requiere crear ni ordenar la agrupación.
Esquema del algoritmo
Podemos usar números enteros para representar sus cadenas binarias, lo que simplifica enormemente el problema de la generación de agrupaciones y la eliminación secuencial de valores. Por ejemplo, con max_len==3
, podemos tomar el número 1--
(donde -
representa el relleno) para significar 4
en decimal. Además, podemos determinar que los números que requieren eliminación si seleccionamos este número son aquellos entre 4
y 4 + 2 ^ x - 1
. Aquí x
es el número de elementos de relleno (2 en este caso), por lo que los números a eliminar están entre 4
y 4 + 2 ^ 2 - 1
(o entre 4
y 7
, expresados como 100
, 110
y 111
).
Para poder resolver exactamente su problema, necesitamos un poco de arrugas, ya que los números que posiblemente serían los mismos en binario podrían ser distintos para algunas partes de su algoritmo. Por ejemplo, 100
, 10-
y 1--
son todos el mismo número, pero deben tratarse de manera diferente en su esquema. En un mundo max_len==3
, tenemos 8 números posibles, pero 14 representaciones posibles:
0 - 000: 0--, 00-
1 - 001:
2 - 010: 01-
3 - 011:
4 - 100: 1--, 10-
5 - 101:
6 - 110: 11-
7 - 111:
Así que 0 y 4 tienen tres codificaciones posibles, 2 y 6 tienen dos, y todas las demás son solo una. Necesitamos generar un grupo de enteros que represente la mayor probabilidad de selección para los números con representación múltiple, así como un mecanismo para rastrear cuántos espacios en blanco incluye el número. Podemos hacer esto agregando algunos bits al final del número para indicar las ponderaciones que queremos. Así que nuestros números se vuelven (usamos dos bits aquí):
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
Algunas propiedades útiles:
-
enc.bits
representa la cantidad de bits que necesitamos para la codificación (dos en este caso) -
int.enc %% enc.bits
nos dice cuántos de los números están explícitamente especificados -
int.enc %/% enc.bits
devuelveint
-
int * 2 ^ enc.bits + explicitly.specified
devuelveint.enc
Tenga en cuenta que explicitly.specified
especifica aquí entre 0
y max_len - 1
en nuestra implementación, ya que siempre hay al menos un dígito especificado. Ahora tenemos una codificación que representa completamente su estructura de datos utilizando solo números enteros. Podemos muestrear desde enteros y reproducir los resultados deseados, con los pesos correctos, etc. Una limitación de este enfoque es que usamos enteros de 32 bits en R, y debemos reservar algunos bits para la codificación, por lo que nos limitamos a grupos de max_len==25
o menos. Usted podría ir más grande si usó números enteros especificados por puntos flotantes de doble precisión, pero no lo hicimos aquí.
Evitando Selecciones Duplicadas
Hay dos formas básicas para asegurarnos de que no elegimos el mismo valor dos veces
- Haga un seguimiento de los valores que quedan disponibles para seleccionar y muestre de forma aleatoria entre ellos
- Muestra aleatoriamente de todos los valores posibles, y luego verifica si el valor ha sido seleccionado previamente, y si lo ha hecho, muestra nuevamente
Si bien la primera opción parece la más limpia, en realidad es muy costosa computacionalmente. Se requiere realizar una exploración vectorial de todos los valores posibles para cada selección para descalificar previamente los valores seleccionados o crear un vector de reducción que contenga valores no descalificados. La opción de reducción solo es más eficiente que el escaneo vectorial si se hace que el vector se reduzca por referencia a través del código C, pero incluso entonces requiere traducciones repetidas de porciones potencialmente grandes del vector, y requiere C.
Aquí usamos el método # 2. Esto nos permite mezclar aleatoriamente el universo de valores posibles una vez, y luego seleccionar cada valor en secuencia, verificar que no haya sido descalificado, y si lo ha hecho, elegir otro, etc. Esto funciona porque es trivial verificar si un el valor ha sido elegido como resultado de nuestra codificación de valor; podemos inferir la ubicación de un valor en una tabla ordenada en base únicamente al valor . Por lo tanto, registramos el estado de cada valor en una tabla ordenada, y podemos actualizar o buscar ese estado a través del acceso directo al índice (no se requiere escaneo).
Ejemplos
La implementación de este algoritmo en base R está disponible como una esencia . Esta implementación particular solo extrae los sorteos completos. Aquí hay una muestra de 10 dibujos de 8 elementos de un grupo 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"
Originalmente teníamos dos implementaciones que dependían del método # 1 para evitar duplicados, una en la base R y otra en C, pero incluso la versión C no es tan rápida como la nueva versión de la base R cuando n
es grande. Estas funciones implementan la capacidad de dibujar sorteos incompletos, por lo que les ofrecemos aquí como referencia:
Puntos de referencia comparativos
Aquí hay un conjunto de puntos de referencia que comparan varias de las funciones que se muestran en esta Q / A. Tiempos en milisegundos. La versión brodie.b
es la que se describe en esta respuesta. brodie
es la implementación original, brodie.C
es el original con algo de C. Todos estos hacen cumplir el requisito de muestras completas. brodie.str
es la versión basada en cadena en la otra respuesta.
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
Esto se adapta relativamente bien a piscinas más grandes.
system.time(sample0110b(18, 100000))
user system elapsed
8.441 0.079 8.527
Notas de referencia:
- frank, y brodie (menos brodie.str) no requieren ninguna agrupación previa, lo que afectará las comparaciones (ver más abajo)
- Josilber es la versión LP
- jbaum es el ejemplo OP
- tensibai se modifica ligeramente para salir en lugar de fallar si el grupo está vacío
- no está configurado para ejecutar Python, por lo que no puedo comparar / tener en cuenta el almacenamiento en búfer
-
-
Representa opciones inviables o demasiado lentas para calcular el tiempo razonablemente.
Los tiempos no incluyen dibujar las agrupaciones ( 0.8
, 2.5
, 401
milisegundos respectivamente para el tamaño 4
, 8
y 16
), lo que es necesario para las jbaum
, josilber
y brodie.str
, ni clasificarlas ( 0.1
, 2.7
, 3700
milisegundos para el tamaño 4
, 8
y 16
), que es necesario para brodie.str
además del sorteo. Si desea incluirlos o no, depende de cuántas veces ejecute la función para un grupo en particular. Además, es casi seguro que hay mejores formas de generar / clasificar las agrupaciones.
Estos son el tiempo medio de tres carreras con microbenchmark
. El código está disponible como una esencia , aunque tenga en cuenta que deberá cargar las sample0110b
, sample0110
, sample01101
y sample01
antemano.
Encontré el problema interesante, así que intenté y obtuve esto con habilidades realmente bajas en R (para que esto pueda mejorarse):
Frank versión editada, gracias al consejo de Frank :
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)
La idea principal es utilizar la comparación de máscara de bits para saber si una muestra es padre de la otra (parte de bit común), si es así, suprima este elemento de la agrupación.
Ahora se necesitan 1,4 s para extraer 100 muestras de un grupo de longitud 16 en mi máquina.
Puede ordenar el grupo para ayudar a decidir qué elementos descalificar. Por ejemplo, mirando un grupo ordenado de tres elementos:
[1] "0" "00" "000" "001" "01" "010" "011" "1" "10" "100" "101" "11"
[13] "110" "111"
Puedo decir que puedo descalificar cualquier cosa que siga a mi elemento seleccionado que tenga más caracteres que mi elemento hasta el primer elemento que tenga el mismo número o menos caracteres. Por ejemplo, si selecciono "01", puedo ver inmediatamente que los siguientes dos elementos ("010", "011") deben eliminarse, pero no el siguiente después de eso porque "1" tiene menos caracteres. Eliminar el "0" después es fácil. Aquí hay una implementación:
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
}
Y una función para hacer grupos ordenados (exactamente como su código, pero los resultados ordenados):
make_pool <- function(size)
sort(
unlist(
lapply(
seq_len(size),
function(x) do.call(paste0, expand.grid(rep(list(c(''0'', ''1'')), x)))
) ) )
Luego, utilizando el grupo max_len
3 (útil para inspeccionar visualmente las cosas se comportan como se espera):
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"
Observe que en el último caso obtenemos todas las combinaciones binarias de 3 dígitos (2 ^ 3) porque por casualidad seguimos muestreando de los de 3 dígitos. Además, con solo un grupo de 3 tamaños hay muchos muestreos que impiden un total de 8 sorteos; usted podría abordar esto con su sugerencia de eliminar las combinaciones que impiden el sorteo completo de la piscina.
Esto es bastante rápido. Mirando el ejemplo de max_len==9
que tomó 2 segundos con la solución 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
Alrededor de medio milisegundo. También puedes intentar razonablemente piscinas bastante grandes:
pool16 <- make_pool(16) # 131K entries
system.time(sample01(pool16, 100))
# user system elapsed
# 3.407 0.146 3.552
Esto no es muy rápido, pero estamos hablando de un grupo con elementos de 130K. También hay potencial para la optimización adicional.
Tenga en cuenta que el paso de clasificación es relativamente lento para los grupos grandes, pero no lo estoy contando porque solo necesita hacerlo una vez, y probablemente pueda encontrar un algoritmo razonable para producir el grupo ordenado previamente.
También existe la posibilidad de un enfoque mucho más rápido de enteros a binarios que exploré en una respuesta ahora eliminada, pero eso requiere un poco más de trabajo para vincular exactamente lo que está buscando.
Está en python y no en r, pero jbaums dijo que estaría bien.
Así que aquí está mi contribución, vea los comentarios en la fuente para obtener una explicación de las partes cruciales.
Todavía estoy trabajando en la solución analítica para determinar la cantidad de posibles combinaciones c
para un árbol de muestras t
y S
de profundidad, por lo que puedo mejorar los combs
la función. Tal vez alguien lo tiene? Este es realmente el cuello de botella ahora.
El muestreo de 100 nodos de un árbol con profundidad 16 toma alrededor de 8 ms en mi computadora portátil. No es la primera vez, pero se acelera hasta cierto punto cuanto más se muestrea porque se llena 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''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()
La salida de referencia:
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 muestras de tamaño 10 de un árbol con profundidad 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'']
Mapeo de ids a cadenas. Puede asignar números a sus vectores 0/1, como mencionó @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
}
Encontrar ids para soltar. Vamos a soltar secuencialmente id
desde el grupo de muestreo:
# 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)
}
Los índices están agrupados por g
, el número de caracteres nchar(get_01str(id))
,. Debido a que los índices están ordenados por g
, g=findInterval(id,cuts)
es una ruta más rápida.
Un índice en el grupo g
, 1 < g < max_len
tiene un índice de "niño" de tamaño g-1
y dos índices de matrices de tamaño g+1
. Para cada nodo hijo, tomamos su nodo hijo hasta que golpeamos g==1
; y para cada nodo principal, tomamos su par de nodos principales hasta que pulsamos g==max_len
.
La estructura del árbol es la más simple en términos del identificador dentro del grupo gid
,. gid
mapas a los dos padres, gid
y gid+2^g
; e invirtiendo este mapeo uno encuentra al niño.
Muestreo
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
}
La oklens
parte integra la idea del OP para omitir cadenas que están garantizadas para hacer imposible el muestreo. Sin embargo, incluso haciendo esto, podemos seguir un camino de muestreo que nos deja sin más opciones. Tomando el ejemplo de la OP de max_len=4
y n=10
, sabemos que hay que bajar 0
y 1
de la consideración, pero ¿Qué pasa si nuestro primer cuatro empates son 00
, 01
, 11
y 10
? Oh bueno, supongo que estamos fuera de suerte. Esta es la razón por la que debería definir las probabilidades de muestreo. (El OP tiene otra idea, para determinar, en cada paso, qué nodos conducirán a un estado imposible, pero eso parece una tarea difícil).
Ilustración
# 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" "" ""
Puntos de referencia (más antiguos que los de la respuesta de @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
Otras respuestas. Hay otras dos respuestas:
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)
Dejé de lado el método de enumeración de @josilber porque tomó demasiado tiempo; y el método numérico / índice de @ BrodieG porque no funcionó en ese momento, pero ahora sí. Vea la respuesta actualizada de @BrodieG para más evaluaciones comparativas.
Velocidad vs corrección. Si bien las respuestas de @josilber son mucho más lentas (y para el método de enumeración, al parecer mucho más intensivas en memoria), están garantizadas para obtener una muestra de tamaño n
en el primer intento. Con el método de cuerdas de @BrodieG o esta respuesta, tendrás que volver a muestrear una y otra vez con la esperanza de dibujar un completo n
. Con grandes max_len
, eso no debería ser un problema, supongo.
Esta respuesta se escala mejor que bro_string
porque no requiere la construcción del pool
frente.
Un enfoque sería simplemente generar todas las tuplas posibles del tamaño apropiado con un enfoque iterativo:
- Construye todas las tuplas de tamaño 1 (todos los elementos en
pool
) - Tome un producto cruzado con elementos en
pool
- Remueve cualquier tupla usando el mismo elemento de
pool
más de una vez. - Eliminar cualquier duplicado exacto de otra tupla
- Retire cualquier tupla con un par que no pueda usarse juntos
- Enjuague y repita hasta que obtenga el tamaño de tupla apropiado
Esto es ejecutable para los tamaños dados ( pool
de longitud 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
Ahora puede muestrear tantas veces como desee de las filas de template
(que será muy rápido), y esto será lo mismo que muestrear al azar del espacio definido.
Si no desea generar el conjunto de todas las tuplas posibles y luego muestrearlas de forma aleatoria (lo que, como puede observar, puede ser imposible para tamaños de entrada grandes), otra opción es dibujar una muestra única con programación de enteros. Básicamente, podría asignar cada elemento en pool
un valor aleatorio y luego seleccionar la tupla factible con la mayor suma de valores. Esto debería dar a cada tupla la misma probabilidad de ser seleccionado porque todos son del mismo tamaño y sus valores fueron seleccionados al azar. Las restricciones del modelo garantizarían que ninguno de los pares de tuplas no permitidos estén seleccionados y que se seleccione el número correcto de elementos.
Aquí hay una solución con el lpSolve
paquete:
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"
Esto parece escalar a piscinas bastante grandes. Por ejemplo, se tarda un poco más de 2 segundos en obtener una muestra de longitud 20 de un grupo de tamaño 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
Si realmente necesita resolver los problemas de los problemas realmente grandes, puede pasar de los solucionadores de código no abierto que se envían lpSolve
a un solucionador comercial como gurobi o cplex (no es gratis en general, pero es gratuito para uso académico).