algorithm - valor - Algoritmo para matrices binarias cuadradas aleatorias equiprobables con dos no ceros no adyacentes en cada fila y columna

suma de filas y columnas de una matriz en java

Question

Sería genial si alguien pudiera apuntarme hacia un algoritmo que me permitiera:

1. crear una matriz cuadrada aleatoria, con las entradas 0 y 1, de modo que
2. Cada fila y cada columna contienen exactamente dos entradas distintas de cero,
3. dos entradas distintas de cero no pueden ser adyacentes,
4. Todas las matrices posibles son equiprobables.

En este momento, logro alcanzar los puntos 1 y 2 haciendo lo siguiente: una matriz de este tipo puede transformarse, utilizando permutaciones adecuadas de filas y columnas, en una matriz de bloques diagonal con bloques de la forma

1 1 0 0 ... 0
0 1 1 0 ... 0
0 0 1 1 ... 0
.............
1 0 0 0 ... 1

Así que empiezo desde una matriz de este tipo utilizando una partición de [0, ..., n-1] y la aleatorizo ​​permutando filas y columnas al azar. Desafortunadamente, no puedo encontrar una manera de integrar la condición de adyacencia, y estoy bastante seguro de que mi algoritmo no tratará todas las matrices por igual.

Actualizar

Me las arreglé para lograr el punto 3. La respuesta fue realmente directa: la matriz de bloques que estoy creando contiene toda la información necesaria para tener en cuenta la condición de adyacencia. Primero algunas propiedades y definiciones:

• una matriz adecuada define las permutaciones de [1, ..., n] que se pueden construir así: seleccione un 1 en la fila 1 . La columna que contiene esta entrada contiene exactamente otra entrada igual a 1 en una fila diferente de 1. De nuevo, la fila a contiene otra entrada 1 en una columna que contiene una segunda entrada 1 en una fila b , y así sucesivamente. Esto inicia una permutación 1 -> a -> b ...

Por ejemplo, con la siguiente matriz, comenzando con la entrada marcada

v
1 0 1 0 0 0 | 1
0 1 0 0 0 1 | 2
1 0 0 1 0 0 | 3
0 0 1 0 1 0 | 4
0 0 0 1 0 1 | 5
0 1 0 0 1 0 | 6
------------+--
1 2 3 4 5 6 |

obtenemos la permutación 1 -> 3 -> 5 -> 2 -> 6 -> 4 -> 1 .

• los ciclos de tal permutación llevan a la matriz de bloques que mencioné anteriormente. También mencioné mezclar la matriz de bloques utilizando permutaciones arbitrarias en las filas y columnas para reconstruir una matriz compatible con los requisitos.

Pero estaba usando cualquier permutación, lo que llevó a algunas entradas adyacentes distintas de cero. Para evitar eso, tengo que elegir permutaciones que separen filas (y columnas) adyacentes en la matriz de bloques. En realidad, para ser más precisos, si dos filas pertenecen a un mismo bloque y son cíclicamente consecutivas (la primera y la última fila de un bloque también se consideran consecutivas), entonces la permutación que quiero aplicar tiene que mover estas filas a no consecutivas. filas de la matriz final (llamaré dos filas incompatibles en ese caso).

Así que la pregunta es: ¿Cómo construir todas esas permutaciones?

La idea más simple es construir una permutación progresivamente mediante la adición aleatoria de filas que sean compatibles con la anterior. Como ejemplo, considere el caso n = 6 usando la partición 6 = 3 + 3 y la matriz de bloques correspondiente

1 1 0 0 0 0 | 1
0 1 1 0 0 0 | 2
1 0 1 0 0 0 | 3
0 0 0 1 1 0 | 4
0 0 0 0 1 1 | 5
0 0 0 1 0 1 | 6
------------+--
1 2 3 4 5 6 |

Aquí las filas 1 , 2 y 3 son mutuamente incompatibles, como lo son 4 , 5 y 6 . Elija una fila aleatoria, digamos 3 .

Escribiremos una permutación como una matriz: [2, 5, 6, 4, 3, 1] significa 1 -> 2 , 2 -> 5 , 3 -> 6 , ... Esto significa que la fila 2 de la matriz de bloques se convertirá en la primera fila de la matriz final, la fila 5 se convertirá en la segunda fila, y así sucesivamente.

Ahora construyamos una permutación adecuada eligiendo aleatoriamente una fila, digamos 3 :

• p = [3, ...]

La siguiente fila se elegirá al azar entre las filas restantes que sean compatibles con 3 : 4 , 5 y 6 . Digamos que elegimos 4 :

• p = [3, 4, ...]

La siguiente elección debe hacerse entre 1 y 2 , por ejemplo 1 :

• p = [3, 4, 1, ...]

Y así sucesivamente: p = [3, 4, 1, 5, 2, 6] .

Aplicando esta permutación a la matriz de bloques, obtenemos:

1 0 1 0 0 0 | 3
0 0 0 1 1 0 | 4
1 1 0 0 0 0 | 1
0 0 0 0 1 1 | 5
0 1 1 0 0 0 | 2
0 0 0 1 0 1 | 6
------------+--
1 2 3 4 5 6 |

Al hacerlo, conseguimos aislar verticalmente todas las entradas que no son cero. Lo mismo se debe hacer con las columnas, por ejemplo, utilizando la permutación p'' = [6, 3, 5, 1, 4, 2] para obtener finalmente

0 1 0 1 0 0 | 3
0 0 1 0 1 0 | 4
0 0 0 1 0 1 | 1
1 0 1 0 0 0 | 5
0 1 0 0 0 1 | 2
1 0 0 0 1 0 | 6
------------+--
6 3 5 1 4 2 |

Así que esto parece funcionar de manera bastante eficiente, pero la construcción de estas permutaciones debe hacerse con precaución, porque se puede bloquear fácilmente: por ejemplo, con n=6 y partición 6 = 2 + 2 + 2 , siguiendo las reglas de construcción establecidas anteriormente puede llevar a p = [1, 3, 2, 4, ...] . Desafortunadamente, 5 y 6 son incompatibles, por lo que elegir uno u otro hace que la última opción sea imposible. Creo que he encontrado todas las situaciones que llevan a un callejón sin salida. Indicaré por r el conjunto de opciones restantes:

1. p = [..., x, ?] , r = {y} con x e y incompatibles
2. p = [..., x, ?, ?] , r = {y, z} con y y z siendo ambos incompatibles con x (no se puede elegir)
3. p = [..., ?, ?] , r = {x, y} con x e y incompatibles (cualquier elección llevaría a la situación 1)
4. p = [..., ?, ?, ?] , r = {x, y, z} con x , y z siendo cíclicamente consecutivos (elegir x o z llevaría a la situación 2, eligiendo y a la situación 3)
5. p = [..., w, ?, ?, ?] , r = {x, y, z} siendo xwy un ciclo de 3 ciclos (no se puede elegir ni x ni y , elegir z llevaría a la situación 3)
6. p = [..., ?, ?, ?, ?] , r = {w, x, y, z} siendo wxyz un ciclo de 4 (cualquier elección llevaría a la situación 4)
7. p = [..., ?, ?, ?, ?] , r = {w, x, y, z} con xyz como un ciclo de 3 (elegir w llevaría a la situación 4, elegir cualquier otro conduciría a la situación 4)

Ahora parece que el siguiente algoritmo da todas las permutaciones adecuadas:

• Mientras haya estrictamente más de 5 números para elegir, elija al azar entre los compatibles.
• Si quedan 5 números para elegir: si los números restantes contienen un ciclo de 3 o un ciclo de 4, rompa ese ciclo (es decir, elija un número que pertenezca a ese ciclo).
• Si quedan 4 números para elegir: si los números restantes contienen tres números cíclicamente consecutivos, elija uno de ellos.
• Si quedan 3 números para elegir: si los números restantes contienen dos números cíclicamente consecutivos, elija uno de ellos.

Estoy bastante seguro de que esto me permite generar todas las permutaciones adecuadas y, por lo tanto, todas las matrices adecuadas.

Desafortunadamente, cada matriz se obtendrá varias veces, dependiendo de la partición elegida.

Answer 1

Introducción

Aquí hay un enfoque prototipo, tratando de resolver la tarea más general del muestreo combinatorio uniforme , que para nuestro enfoque aquí significa: podemos usar este enfoque para todo lo que podamos formular como SAT-problem .

No está explotando su problema directamente y toma un gran desvío. Este desvío al problema SAT puede ayudar en lo que respecta a la teoría (resultados teóricos generales más poderosos) y la eficiencia (solucionadores SAT).

Dicho esto, no es un enfoque si desea muestrear en segundos o menos (en mis experimentos), al menos mientras esté preocupado por la uniformidad.

Teoría

El enfoque, basado en los resultados de la teoría de la complejidad, sigue este trabajo :

GOMES, Carla P .; SABHARWAL, Ashish; SELMAN, Bart. Muestreo casi uniforme de espacios combinatorios usando restricciones XOR. En: Avances En Sistemas De Procesamiento De La Información Neural. 2007. S. 481-488.

La idea básica:

• formular el problema como SAT-problema
• agregue xors generados al azar al problema (¡actuando solo sobre las variables de decisión! eso es importante en la práctica)
  • esto reducirá el número de soluciones (algunas soluciones serán imposibles)
  • haga eso en un bucle (con parámetros ajustados) hasta que solo quede una solución!
    • la búsqueda de alguna solución se está realizando mediante solucionadores SAT o # solucionadores SAT (= conteo de modelos)
    • si hay más de una solución: no se agregarán xors, pero se realizará un reinicio completo: ¡agregue random-xors al problema de inicio!

Las garantías:

• Al ajustar los parámetros correctamente, este enfoque logra un muestreo casi uniforme.
  • este ajuste puede ser costoso, ya que se basa en aproximar la cantidad de soluciones posibles
  • Empíricamente esto también puede ser costoso!
  • La respuesta de Ante, al mencionar la secuencia de números A001499 en realidad da un buen límite superior en el espacio de la solución (¡ya que simplemente está ignorando las restricciones de adyacencia!)

Los inconvenientes:

• ineficiente para grandes problemas (en general; no necesariamente en comparación con las alternativas como MCMC y compañía)
  • Necesidad de cambiar / reducir parámetros para producir muestras.
  • Esos parámetros reducidos pierden las garantías teóricas.
  • pero empíricamente: ¡los buenos resultados son todavía posibles!

Parámetros:

En la práctica, los parámetros son:

• N: número de xors añadidos
• L: número mínimo de variables que forman parte de una restricción xor
• U: número máximo de variables que forman parte de una restricción xor

N es importante para reducir el número de soluciones posibles. Dada la constante de N, las otras variables, por supuesto, también tienen algún efecto sobre eso.

La teoría dice (si interpreto correctamente), que deberíamos usar L = R = 0.5 * # dec-vars.

Esto es imposible en la práctica aquí, ya que las xor restricciones afectan mucho a las personas que resuelven SAT.

Here algunas diapositivas más científicas sobre el impacto de L y U.

Llaman xors de tamaño 8-20 corto-XORS, ¡mientras que necesitaremos usar aún más cortos más tarde!

Implementación

Versión definitiva

Aquí hay una implementación bastante pirata en Python, usando los scripts XorSample de here .

El solucionador de SAT subyacente en uso es Cryptominisat .

El código básicamente se reduce a:

• Transformar el problema en forma normal conjuntiva.
  • como DIMACS-CNF
• Implementar el enfoque de muestreo:
  • Pide XorSample (basado en tuberías + basado en archivos)
  • Llame a SAT-solver (basado en archivo)
• Añadir muestras a algún archivo para su posterior análisis.

Código: (espero que ya te haya advertido sobre la calidad del código)

from itertools import count
from time import time
import subprocess
import numpy as np
import os
import shelve
import uuid
import pickle
from random import SystemRandom
cryptogen = SystemRandom()
""" Helper functions """
# K-ARY CONSTRAINT GENERATION
# ###########################
# SINZ, Carsten. Towards an optimal CNF encoding of boolean cardinality constraints.
# CP, 2005, 3709. Jg., S. 827-831.
def next_var_index(start):
next_var = start
while(True):
yield next_var
next_var += 1
class s_index():
def __init__(self, start_index):
self.firstEnvVar = start_index
def next(self,i,j,k):
return self.firstEnvVar + i*k +j
def gen_seq_circuit(k, input_indices, next_var_index_gen):
cnf_string = ''''
s_index_gen = s_index(next_var_index_gen.next())
# write clauses of first partial sum (i.e. i=0)
cnf_string += (str(-input_indices[0]) + '' '' + str(s_index_gen.next(0,0,k)) + '' 0/n'')
for i in range(1, k):
cnf_string += (str(-s_index_gen.next(0, i, k)) + '' 0/n'')
# write clauses for general case (i.e. 0 < i < n-1)
for i in range(1, len(input_indices)-1):
cnf_string += (str(-input_indices[i]) + '' '' + str(s_index_gen.next(i, 0, k)) + '' 0/n'')
cnf_string += (str(-s_index_gen.next(i-1, 0, k)) + '' '' + str(s_index_gen.next(i, 0, k)) + '' 0/n'')
for u in range(1, k):
cnf_string += (str(-input_indices[i]) + '' '' + str(-s_index_gen.next(i-1, u-1, k)) + '' '' + str(s_index_gen.next(i, u, k)) + '' 0/n'')
cnf_string += (str(-s_index_gen.next(i-1, u, k)) + '' '' + str(s_index_gen.next(i, u, k)) + '' 0/n'')
cnf_string += (str(-input_indices[i]) + '' '' + str(-s_index_gen.next(i-1, k-1, k)) + '' 0/n'')
# last clause for last variable
cnf_string += (str(-input_indices[-1]) + '' '' + str(-s_index_gen.next(len(input_indices)-2, k-1, k)) + '' 0/n'')
return (cnf_string, (len(input_indices)-1)*k, 2*len(input_indices)*k + len(input_indices) - 3*k - 1)
# K=2 clause GENERATION
# #####################
def gen_at_most_2_constraints(vars, start_var):
constraint_string = ''''
used_clauses = 0
used_vars = 0
index_gen = next_var_index(start_var)
circuit = gen_seq_circuit(2, vars, index_gen)
constraint_string += circuit[0]
used_clauses += circuit[2]
used_vars += circuit[1]
start_var += circuit[1]
return [constraint_string, used_clauses, used_vars, start_var]
def gen_at_least_2_constraints(vars, start_var):
k = len(vars) - 2
vars = [-var for var in vars]
constraint_string = ''''
used_clauses = 0
used_vars = 0
index_gen = next_var_index(start_var)
circuit = gen_seq_circuit(k, vars, index_gen)
constraint_string += circuit[0]
used_clauses += circuit[2]
used_vars += circuit[1]
start_var += circuit[1]
return [constraint_string, used_clauses, used_vars, start_var]
# Adjacency conflicts
# ###################
def get_all_adjacency_conflicts_4_neighborhood(N, X):
conflicts = set()
for x in range(N):
for y in range(N):
if x < (N-1):
conflicts.add(((x,y),(x+1,y)))
if y < (N-1):
conflicts.add(((x,y),(x,y+1)))
cnf = '''' # slow string appends
for (var_a, var_b) in conflicts:
var_a_ = X[var_a]
var_b_ = X[var_b]
cnf += ''-'' + var_a_ + '' '' + ''-'' + var_b_ + '' 0 /n''
return cnf, len(conflicts)
# Build SAT-CNF
#############
def build_cnf(N, verbose=False):
var_counter = count(1)
N_CLAUSES = 0
X = np.zeros((N, N), dtype=object)
for a in range(N):
for b in range(N):
X[a,b] = str(next(var_counter))
# Adjacency constraints
CNF, N_CLAUSES = get_all_adjacency_conflicts_4_neighborhood(N, X)
# k=2 constraints
NEXT_VAR = N*N+1
for row in range(N):
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_most_2_constraints(X[row, :].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_least_2_constraints(X[row, :].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
for col in range(N):
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_most_2_constraints(X[:, col].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_least_2_constraints(X[:, col].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
# build final cnf
CNF = ''p cnf '' + str(NEXT_VAR-1) + '' '' + str(N_CLAUSES) + ''/n'' + CNF
return X, CNF, NEXT_VAR-1
# External tools
# ##############
def get_random_xor_problem(CNF_IN_fp, N_DEC_VARS, N_ALL_VARS, s, min_l, max_l):
# .cnf not part of arg!
p = subprocess.Popen([''./gen-wff'', CNF_IN_fp,
str(N_DEC_VARS), str(N_ALL_VARS),
str(s), str(min_l), str(max_l), ''xored''],
stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
result = p.communicate()
os.remove(CNF_IN_fp + ''-str-xored.xor'') # file not needed
return CNF_IN_fp + ''-str-xored.cnf''
def solve(CNF_IN_fp, N_DEC_VARS):
seed = cryptogen.randint(0, 2147483647) # actually no reason to do it; but can''t hurt either
p = subprocess.Popen(["./cryptominisat5", ''-t'', ''4'', ''-r'', str(seed), CNF_IN_fp], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
result = p.communicate()[0]
sat_line = result.find(''s SATISFIABLE'')
if sat_line != -1:
# solution found!
vars = parse_solution(result)[:N_DEC_VARS]
# forbid solution (DeMorgan)
negated_vars = list(map(lambda x: x*(-1), vars))
with open(CNF_IN_fp, ''a'') as f:
f.write( (str(negated_vars)[1:-1] + '' 0/n'').replace('','', ''''))
# assume solve is treating last constraint despite not changing header!
# solve again
seed = cryptogen.randint(0, 2147483647)
p = subprocess.Popen(["./cryptominisat5", ''-t'', ''4'', ''-r'', str(seed), CNF_IN_fp], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
result = p.communicate()[0]
sat_line = result.find(''s SATISFIABLE'')
if sat_line != -1:
os.remove(CNF_IN_fp) # not needed anymore
return True, False, None
else:
return True, True, vars
else:
return False, False, None
def parse_solution(output):
# assumes there is one
vars = []
for line in output.split("/n"):
if line:
if line[0] == ''v'':
line_vars = list(map(lambda x: int(x), line.split()[1:]))
vars.extend(line_vars)
return vars
# Core-algorithm
# ##############
def xorsample(X, CNF_IN_fp, N_DEC_VARS, N_VARS, s, min_l, max_l):
start_time = time()
while True:
# add s random XOR constraints to F
xored_cnf_fp = get_random_xor_problem(CNF_IN_fp, N_DEC_VARS, N_VARS, s, min_l, max_l)
state_lvl1, state_lvl2, var_sol = solve(xored_cnf_fp, N_DEC_VARS)
print(''------------'')
if state_lvl1 and state_lvl2:
print(''FOUND'')
d = shelve.open(''N_15_70_4_6_TO_PLOT'')
d[str(uuid.uuid4())] = (pickle.dumps(var_sol), time() - start_time)
d.close()
return True
else:
if state_lvl1:
print(''sol not unique'')
else:
print(''no sol found'')
print(''------------'')
""" Run """
N = 15
N_DEC_VARS = N*N
X, CNF, N_VARS = build_cnf(N)
with open(''my_problem.cnf'', ''w'') as f:
f.write(CNF)
counter = 0
while True:
print(''sample: '', counter)
xorsample(X, ''my_problem'', N_DEC_VARS, N_VARS, 70, 4, 6)
counter += 1

La salida se verá como (se eliminaron algunas advertencias):

------------
no sol found
------------
------------
no sol found
------------
------------
no sol found
------------
------------
sol not unique
------------
------------
FOUND

Núcleo: formulación de CNF

Introducimos una variable para cada celda de la matriz. N = 20 significa 400 variables binarias.

Adjancencia:

Precalcular todos los conflictos de reducción de simetría y agregar cláusulas de conflicto.

Teoría básica:

a -> !b
<->
!a v !b (propositional logic)

Cardinalidad Fila / Col:

Esto es difícil de expresar en CNF y los enfoques ingenuos necesitan un número exponencial de restricciones.

Utilizamos alguna codificación basada en el circuito del sumador (SINZ, Carsten. Hacia una codificación CNF óptima de restricciones de cardinalidad booleana) que introduce nuevas variables auxiliares.

Observación:

sum(var_set) <= k
<->
sum(negated(var_set)) >= len(var_set) - k

Estas codificaciones SAT se pueden colocar en contadores de modelos exactos (para N pequeña, por ejemplo, <9). El número de soluciones es igual a los resultados de Ante, ¡lo cual es una fuerte indicación de una transformación correcta!

También hay interesantes contadores de modelos aproximados (también basados ​​en gran medida en las restricciones xor) como approxMC que muestra una cosa más que podemos hacer con la formulación SAT. Pero en la práctica no he podido usarlos (approxMC = autoconf; sin comentarios).

Otros experimentos

También construí una versión usando pblib , para usar formulaciones de cardinalidad más potentes para la formulación SAT-CNF. No intenté usar la API basada en C ++, sino solo el pbencoder reducido, que selecciona automáticamente la mejor codificación, que fue mucho peor que la codificación que usé anteriormente (lo que es mejor todavía es un problema de investigación; a menudo incluso restricciones redundantes poder ayudar).

Análisis empírico

Para obtener un tamaño de muestra (dada mi paciencia ), solo calculé muestras para N = 15 . En este caso utilizamos:

• N = 70 xors
• L, U = 4,6

También calculé algunas muestras para N = 20 con (100,3,6), ¡pero esto toma unos minutos y redujimos el límite inferior!

Visualización

Aquí un poco de animación (fortaleciendo mi relación de amor-odio con matplotlib ):

Edición: Y una comparación (reducida) con el muestreo uniforme de fuerza bruta con N = 5 (NXOR, L, U = 4, 10, 30):

(Todavía no me he decidido sobre la adición del código de trazado. Es tan feo como el anterior y la gente podría ver demasiado en mi confusión estadística; normalizaciones y co).

Teoría

El análisis estadístico es probablemente difícil de realizar ya que el problema subyacente es de naturaleza combinatoria. Incluso no es del todo obvio cómo debe verse el PDF final de la celda. En el caso de N = impar, es probable que no sea uniforme y se vea como un tablero de ajedrez (hice una comprobación de fuerza bruta N = 5 para observar esto).

De una cosa podemos estar seguros (imho): ¡la simetría!

Dada una matriz de célula-PDF, debemos esperar que la matriz sea simétrica (A = AT). Esto se verifica en la visualización y se grafica la norma euclidiana de las diferencias en el tiempo.

Podemos hacer lo mismo en otra observación: pares observados.

Para N = 3, podemos observar los siguientes pares:

• 0,1
• 0,2
• 1,2

¡Ahora podemos hacer esto por fila y por columna y también debemos esperar simetría!

Lamentablemente, probablemente no sea fácil decir algo sobre la variación y, por lo tanto, ¡las muestras necesarias para hablar de confianza!

Observación

De acuerdo con mi percepción simplificada, las muestras actuales y el PDF de la celda se ven bien, aunque la convergencia aún no se alcanza (o estamos muy lejos de la uniformidad).

El aspecto más importante son probablemente las dos normas, que disminuyen agradablemente hacia 0. (sí; se podría ajustar un algoritmo para eso al transponer con prob = 0.5; pero esto no se hace aquí ya que no serviría para nada).

Posibles próximos pasos

• Parámetros de sintonía
• Echa un vistazo al enfoque utilizando # solversores / contadores de modelos en lugar de solversores SAT
• Pruebe diferentes formulaciones de CNF, especialmente en lo que respecta a codificaciones de cardinalidad y codificaciones xor
  • XorSample es por defecto el uso de la codificación similar a tseitin para crecer exponencialmente
    • para los xors más pequeños (como se usan) podría ser una buena idea usar codificación ingenua (que se propaga más rápido)
      • XorSample apoya eso en teoría; Pero el guión funciona de manera diferente en la práctica.
      • Cryptominisat es conocido por su manejo XOR dedicado (ya que fue creado para analizar criptografía incluyendo muchos Xors) y podría ganar algo por codificación ingenua (como inferir Xors de CNFs inflados es mucho más difícil)
• Más análisis estadístico
• Deshazte de los scripts XorSample (shell + perl ...)

Resumen

• El enfoque es muy general.
• Este código produce muestras factibles.
• No debería ser difícil de probar, que todas las soluciones posibles pueden ser muestreadas
• Otros tienen probadas garantías teóricas de uniformidad para algunos parámetros.
  • no se mantiene para nuestros params
• Otros han analizado empíricamente / teóricamente parámetros más pequeños (en uso aquí)

Answer 2

(Resultados de prueba actualizados, ejemplo de ejecución y fragmentos de código a continuación).

Puede usar la programación dinámica para calcular el número de soluciones que resultan de cada estado (de una manera mucho más eficiente que un algoritmo de fuerza bruta) y usar esos valores (precalculados) para crear soluciones aleatorias equiprobables.

Considere el ejemplo de una matriz de 7x7; Al inicio, el estado es:

0,0,0,0,0,0,0 

lo que significa que hay siete columnas adyacentes no utilizadas. Después de agregar dos unidades a la primera fila, el estado podría ser, por ejemplo:

0,1,0,0,1,0,0 

Con dos columnas que ahora tienen una. Después de agregar unos a la segunda fila, el estado podría ser, por ejemplo:

0,1,1,0,1,0,1 

Después de que se llenen tres filas, existe la posibilidad de que una columna tenga su máximo de dos unidades; esto divide efectivamente la matriz en dos zonas independientes:

1,1,1,0,2,0,1 -> 1,1,1,0 + 0,1 

Estas zonas son independientes en el sentido de que la regla de no-adyacentes no tiene ningún efecto cuando se agregan unas a diferentes zonas, y el orden de las zonas no tiene ningún efecto en el número de soluciones.

Para utilizar estos estados como firmas para tipos de soluciones, tenemos que transformarlos en una notación canónica. Primero, debemos tener en cuenta el hecho de que las columnas con solo 1 una en ellas pueden ser inutilizables en la siguiente fila, porque contienen una en la fila actual. Entonces, en lugar de una notación binaria, tenemos que usar una notación ternaria, por ejemplo:

2,1,1,0 + 0,1 

donde 2 significa que esta columna se usó en la fila actual (y no que hay 2 unos en la columna). En el siguiente paso, deberíamos convertir los dos de nuevo en unos.

Además, también podemos reflejar los grupos separados para colocarlos en su notación lexicográficamente más pequeña:

2,1,1,0 + 0,1 -> 0,1,1,2 + 0,1 

Por último, clasificamos los grupos separados de pequeños a grandes, y luego lexicográficamente, de modo que un estado en una matriz más grande puede ser, por ejemplo:

0,0 + 0,1 + 0,0,2 + 0,1,0 + 0,1,0,1 

Luego, al calcular el número de soluciones que resultan de cada estado, podemos usar la memoria mediante la notación canónica de cada estado como clave.

Crear un diccionario de los estados y la cantidad de soluciones para cada uno de ellos solo debe hacerse una vez, y una tabla para matrices más grandes probablemente también se puede usar para matrices más pequeñas.

En la práctica, generaría un número aleatorio entre 0 y el número total de soluciones, y luego, para cada fila, observaría los diferentes estados que podría crear a partir del estado actual, la cantidad de soluciones únicas que cada uno genere y vea qué opción lleva a la solución que corresponde con su número generado aleatoriamente.

Tenga en cuenta que cada estado y la clave correspondiente solo pueden aparecer en una fila en particular, por lo que puede almacenar las claves en diccionarios separados por fila.

RESULTADOS DE LA PRUEBA

Una primera prueba con JavaScript no optimizado dio resultados muy prometedores. Con la programación dinámica, calcular el número de soluciones para una matriz de 10x10 ahora toma un segundo, donde un algoritmo de fuerza bruta tomó varias horas (y esta es la parte del algoritmo que solo se necesita hacer una vez). El tamaño del diccionario con las firmas y el número de soluciones crece con un factor decreciente que se aproxima a 2.5 para cada paso en tamaño; El tiempo para generarlo crece con un factor de alrededor de 3.

Estos son el número de soluciones, estados, firmas (tamaño total de los diccionarios) y el número máximo de firmas por fila (el diccionario más grande por fila) que se crean:

size unique solutions states signatures max/row
4x4 2 9 6 2
5x5 16 73 26 8
6x6 722 514 107 40
7x7 33,988 2,870 411 152
8x8 2,215,764 13,485 1,411 596
9x9 179,431,924 56,375 4,510 1,983
10x10 17,849,077,140 218,038 13,453 5,672
11x11 2,138,979,146,276 801,266 38,314 14,491
12x12 304,243,884,374,412 2,847,885 104,764 35,803
13x13 50,702,643,217,809,908 9,901,431 278,561 96,414
14x14 9,789,567,606,147,948,364 33,911,578 723,306 238,359
15x15 2,168,538,331,223,656,364,084 114,897,838 1,845,861 548,409
16x16 546,386,962,452,256,865,969,596 4,952,501 1,444,487
17x17 155,420,047,516,794,379,573,558,433 12,837,870 3,754,040
18x18 48,614,566,676,379,251,956,711,945,475 31,452,747 8,992,972
19x19 17,139,174,923,928,277,182,879,888,254,495 74,818,773 20,929,008
20x20 6,688,262,914,418,168,812,086,412,204,858,650 175,678,000 50,094,203

(Resultados adicionales obtenidos con C ++, utilizando una implementación simple de enteros de 128 bits).

EJEMPLO

El diccionario para una matriz de 5x5 se ve así:

row 0: 00000 -> 16 row 3: 101 -> 0
1112 -> 1
row 1: 20002 -> 2 1121 -> 1
00202 -> 4 1+01 -> 0
02002 -> 2 11+12 -> 2
02020 -> 2 1+121 -> 1
0+1+1 -> 0
row 2: 10212 -> 1 1+112 -> 1
12012 -> 1
12021 -> 2 row 4: 0 -> 0
12102 -> 1 11 -> 0
21012 -> 0 12 -> 0
02121 -> 3 1+1 -> 1
01212 -> 1 1+2 -> 0

El número total de soluciones es 16; Si elegimos aleatoriamente un número del 0 al 15, por ejemplo, 13, podemos encontrar la solución correspondiente (es decir, la 14ª) de esta manera:

state: 00000
options: 10100 10010 10001 01010 01001 00101
signature: 00202 02002 20002 02020 02002 00202
solutions: 4 2 2 2 2 4 

Esto nos dice que la 14ª solución es la segunda solución de la opción 00101. El siguiente paso es:

state: 00101
options: 10010 01010
signature: 12102 02121
solutions: 1 3 

Esto nos dice que la segunda solución es la primera solución de la opción 01010. El siguiente paso es:

state: 01111
options: 10100 10001 00101
signature: 11+12 1112 1+01
solutions: 2 1 0 

Esto nos dice que la primera solución es la primera solución de la opción 10100. El siguiente paso es:

state: 11211
options: 01010 01001
signature: 1+1 1+1
solutions: 1 1 

Esto nos dice que la primera solución es la primera solución de la opción 01010. El último paso es:

state: 12221
options: 10001 

Y la matriz 5x5 correspondiente al número 13 elegido al azar es:

0 0 1 0 1
0 1 0 1 0
1 0 1 0 0
0 1 0 1 0
1 0 0 0 1 

Y aquí hay un ejemplo de código rápido; ejecute el fragmento para generar la firma y el diccionario de recuento de soluciones, y genere una matriz aleatoria de 10x10 (se tarda un segundo en generar el diccionario; una vez hecho esto, genera soluciones aleatorias en medio milisegundo):

function signature(state, prev) {
var zones = [], zone = [];
for (var i = 0; i < state.length; i++) {
if (state[i] == 2) {
if (zone.length) zones.push(mirror(zone));
zone = [];
}
else if (prev[i]) zone.push(3);
else zone.push(state[i]);
}
if (zone.length) zones.push(mirror(zone));
zones.sort(function(a,b) {return a.length - b.length || a - b;});
return zones.length ? zones.join("2") : "2";
function mirror(zone) {
var ltr = zone.join('''');
zone.reverse();
var rtl = zone.join('''');
return (ltr < rtl) ? ltr : rtl;
}
}
function memoize(n) {
var memo = [], empty = [];
for (var i = 0; i <= n; i++) memo[i] = [];
for (var i = 0; i < n; i++) empty[i] = 0;
memo[0][signature(empty, empty)] = next_row(empty, empty, 1);
return memo;
function next_row(state, prev, row) {
if (row > n) return 1;
var solutions = 0;
for (var i = 0; i < n - 2; i++) {
if (state[i] == 2 || prev[i] == 1) continue;
for (var j = i + 2; j < n; j++) {
if (state[j] == 2 || prev[j] == 1) continue;
var s = state.slice(), p = empty.slice();
++s[i]; ++s[j]; ++p[i]; ++p[j];
var sig = signature(s, p);
var sol = memo[row][sig];
if (sol == undefined)
memo[row][sig] = sol = next_row(s, p, row + 1);
solutions += sol;
}
}
return solutions;
}
}
function random_matrix(n, memo) {
var matrix = [], empty = [], state = [], prev = [];
for (var i = 0; i < n; i++) empty[i] = state[i] = prev[i] = 0;
var total = memo[0][signature(empty, empty)];
var pick = Math.floor(Math.random() * total);
document.write("solution " + pick.toLocaleString(''en-US'') +
" from a total of " + total.toLocaleString(''en-US'') + "<br>");
for (var row = 1; row <= n; row++) {
var options = find_options(state, prev);
for (var i in options) {
var state_copy = state.slice();
for (var j in state_copy) state_copy[j] += options[i][j];
var sig = signature(state_copy, options[i]);
var solutions = memo[row][sig];
if (pick < solutions) {
matrix.push(options[i].slice());
prev = options[i].slice();
state = state_copy.slice();
break;
}
else pick -= solutions;
}
}
return matrix;
function find_options(state, prev) {
var options = [];
for (var i = 0; i < n - 2; i++) {
if (state[i] == 2 || prev[i] == 1) continue;
for (var j = i + 2; j < n; j++) {
if (state[j] == 2 || prev[j] == 1) continue;
var option = empty.slice();
++option[i]; ++option[j];
options.push(option);
}
}
return options;
}
}
var size = 10;
var memo = memoize(size);
var matrix = random_matrix(size, memo);
for (var row in matrix) document.write(matrix[row] + "<br>");

El siguiente fragmento de código muestra el diccionario de firmas y los recuentos de soluciones para una matriz de tamaño 10x10. He utilizado un formato de firma ligeramente diferente de la explicación anterior: las zonas están delimitadas por un ''2'' en lugar de un signo más, y una columna que tiene un uno en la fila anterior está marcada con un ''3'' en lugar de un ''2''. Esto muestra cómo las claves podrían almacenarse en un archivo como enteros con 2 × N bits (rellenados con 2).

function signature(state, prev) {
var zones = [], zone = [];
for (var i = 0; i < state.length; i++) {
if (state[i] == 2) {
if (zone.length) zones.push(mirror(zone));
zone = [];
}
else if (prev[i]) zone.push(3);
else zone.push(state[i]);
}
if (zone.length) zones.push(mirror(zone));
zones.sort(function(a,b) {return a.length - b.length || a - b;});
return zones.length ? zones.join("2") : "2";
function mirror(zone) {
var ltr = zone.join('''');
zone.reverse();
var rtl = zone.join('''');
return (ltr < rtl) ? ltr : rtl;
}
}
function memoize(n) {
var memo = [], empty = [];
for (var i = 0; i <= n; i++) memo[i] = [];
for (var i = 0; i < n; i++) empty[i] = 0;
memo[0][signature(empty, empty)] = next_row(empty, empty, 1);
return memo;
function next_row(state, prev, row) {
if (row > n) return 1;
var solutions = 0;
for (var i = 0; i < n - 2; i++) {
if (state[i] == 2 || prev[i] == 1) continue;
for (var j = i + 2; j < n; j++) {
if (state[j] == 2 || prev[j] == 1) continue;
var s = state.slice(), p = empty.slice();
++s[i]; ++s[j]; ++p[i]; ++p[j];
var sig = signature(s, p);
var sol = memo[row][sig];
if (sol == undefined)
memo[row][sig] = sol = next_row(s, p, row + 1);
solutions += sol;
}
}
return solutions;
}
}
var memo = memoize(10);
for (var i in memo) {
document.write("row " + i + ":<br>");
for (var j in memo[i]) {
document.write("&quot;" + j + "&quot;: " + memo[i][j] + "<br>");
}
}

Answer 3

Aquí hay un enfoque muy rápido de generar la matriz fila por fila, escrito en Java:

public static void main(String[] args) throws Exception {
int n = 100;
Random rnd = new Random();
byte[] mat = new byte[n*n];
byte[] colCount = new byte[n];
//generate row by row
for (int x = 0; x < n; x++) {
//generate a random first bit
int b1 = rnd.nextInt(n);
while ( (x > 0 && mat[(x-1)*n + b1] == 1) || //not adjacent to the one above
(colCount[b1] == 2) //not in a column which has 2
) b1 = rnd.nextInt(n);
//generate a second bit, not equal to the first one
int b2 = rnd.nextInt(n);
while ( (b2 == b1) || //not the same as bit 1
(x > 0 && mat[(x-1)*n + b2] == 1) || //not adjacent to the one above
(colCount[b2] == 2) || //not in a column which has 2
(b2 == b1 - 1) || //not adjacent to b1
(b2 == b1 + 1)
) b2 = rnd.nextInt(n);
//fill the matrix values and increment column counts
mat[x*n + b1] = 1;
mat[x*n + b2] = 1;
colCount[b1]++;
colCount[b2]++;
}
String arr = Arrays.toString(mat).substring(1, n*n*3 - 1);
System.out.println(arr.replaceAll("(.{" + n*3 + "})", "$1/n"));
}

Esencialmente genera cada una fila aleatoria a la vez. Si la fila viola cualquiera de las condiciones, se genera de nuevo (de nuevo al azar). Creo que esto va a satisfacer la condición 4 también.

Agregando una nota rápida de que girará para siempre en Ns donde no hay soluciones (como N = 3).

Answer 4

Esta solución utiliza la recursión para establecer la celda de la matriz una por una. Si la caminata aleatoria termina con una solución imposible, entonces retrocedemos un paso en el árbol y continuamos la caminata aleatoria.

El algoritmo es eficiente y creo que los datos generados son altamente equiprobable.

package rndsqmatrix;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.IntStream;
public class RndSqMatrix {
/**
* Generate a random matrix
* @param size the size of the matrix
* @return the matrix encoded in 1d array i=(x+y*size)
*/
public static int[] generate(final int size) {
return generate(size, new int[size * size], new int[size],
new int[size]);
}
/**
* Build a matrix recursivly with a random walk
* @param size the size of the matrix
* @param matrix the matrix encoded in 1d array i=(x+y*size)
* @param rowSum
* @param colSum
* @return
*/
private static int[] generate(final int size, final int[] matrix,
final int[] rowSum, final int[] colSum) {
// generate list of valid positions
final List<Integer> positions = new ArrayList();
for (int y = 0; y < size; y++) {
if (rowSum[y] < 2) {
for (int x = 0; x < size; x++) {
if (colSum[x] < 2) {
final int p = x + y * size;
if (matrix[p] == 0
&& (x == 0 || matrix[p - 1] == 0)
&& (x == size - 1 || matrix[p + 1] == 0)
&& (y == 0 || matrix[p - size] == 0)
&& (y == size - 1 || matrix[p + size] == 0)) {
positions.add(p);
}
}
}
}
}
// no valid positions ?
if (positions.isEmpty()) {
// if the matrix is incomplete => return null
for (int i = 0; i < size; i++) {
if (rowSum[i] != 2 || colSum[i] != 2) {
return null;
}
}
// the matrix is complete => return it
return matrix;
}
// random walk
Collections.shuffle(positions);
for (int p : positions) {
// set ''1'' and continue recursivly the exploration
matrix[p] = 1;
rowSum[p / size]++;
colSum[p % size]++;
final int[] solMatrix = generate(size, matrix, rowSum, colSum);
if (solMatrix != null) {
return solMatrix;
}
// rollback
matrix[p] = 0;
rowSum[p / size]--;
colSum[p % size]--;
}
// we can''t find a valid matrix from here => return null
return null;
}
public static void printMatrix(final int size, final int[] matrix) {
for (int y = 0; y < size; y++) {
for (int x = 0; x < size; x++) {
System.out.print(matrix[x + y * size]);
System.out.print(" ");
}
System.out.println();
}
}
public static void printStatistics(final int size, final int count) {
final int sumMatrix[] = new int[size * size];
for (int i = 0; i < count; i++) {
final int[] matrix = generate(size);
for (int j = 0; j < sumMatrix.length; j++) {
sumMatrix[j] += matrix[j];
}
}
printMatrix(size, sumMatrix);
}
public static void checkAlgorithm() {
final int size = 8;
final int count = 2215764;
final int divisor = 122;
final int sumMatrix[] = new int[size * size];
for (int i = 0; i < count/divisor ; i++) {
final int[] matrix = generate(size);
for (int j = 0; j < sumMatrix.length; j++) {
sumMatrix[j] += matrix[j];
}
}
int total = 0;
for(int i=0; i < sumMatrix.length; i++) {
total += sumMatrix[i];
}
final double factor = (double)total / (count/divisor);
System.out.println("Factor=" + factor + " (theory=16.0)");
}
public static void benchmark(final int size, final int count,
final boolean parallel) {
final long begin = System.currentTimeMillis();
if (!parallel) {
for (int i = 0; i < count; i++) {
generate(size);
}
} else {
IntStream.range(0, count).parallel().forEach(i -> generate(size));
}
final long end = System.currentTimeMillis();
System.out.println("rate="
+ (double) (end - begin) / count + "ms/matrix");
}
public static void main(String[] args) {
checkAlgorithm();
benchmark(8, 10000, true);
//printStatistics(8, 2215764/36);
printStatistics(8, 2215764);
}
}

La salida es:

Factor=16.0 (theory=16.0)
rate=0.2835ms/matrix
552969 554643 552895 554632 555680 552753 554567 553389
554071 554847 553441 553315 553425 553883 554485 554061
554272 552633 555130 553699 553604 554298 553864 554028
554118 554299 553565 552986 553786 554473 553530 554771
554474 553604 554473 554231 553617 553556 553581 553992
554960 554572 552861 552732 553782 554039 553921 554661
553578 553253 555721 554235 554107 553676 553776 553182
553086 553677 553442 555698 553527 554850 553804 553444

Answer 5

Sólo unos pocos pensamientos. Número de matrices que satisfacen las condiciones para n <= 10:

3 0
4 2
5 16
6 722
7 33988
8 2215764
9 179431924
10 17849077140

Desafortunadamente no hay una secuencia con estos números en OEIS .

Hay uno similar ( A001499 ), sin condición para el vecino. El número de matrices nxn en este caso es ''de orden'' como el número de matrices (n-1) x (n-1) de A001499. Eso es lo que se espera, ya que hay varias formas de llenar una fila en este caso, la posición 2 está en n lugares con al menos un cero entre ellas ((n-1) elige 2). Igual que la posición 2 de uno en (n-1) lugares sin la restricción.

No creo que haya una conexión fácil entre estas matrices de orden n y A001499 matriz de orden n-1, lo que significa que si tenemos matriz A001499 podemos construir algunas de estas matrices.

Con esto, para n = 20, el número de matrices es> 10 ^ 30. Bastante :-/