Cuenta el número de matrices Hankelable

12

Antecedentes

Una matriz binaria de Hankel es una matriz con diagonales oblicuas constantes (diagonales de pendiente positiva) que contiene solo 0sy 1s. Por ejemplo, una matriz de Hankel binaria 5x5 se parece a

a b c d e
b c d e f
c d e f g
d e f g h
e f g h i

donde a, b, c, d, e, f, g, h, ison o 0o 1.

Definamos una matriz M como Hankelable si hay una permutación del orden de las filas y columnas de M para que M sea ​​una matriz de Hankel. Esto significa que uno puede aplicar una permutación al orden de las filas y una posiblemente diferente a las columnas.

El reto

El desafío es contar cuántos Hankelable n por nmatrices hay para todos ncon el mayor valor posible.

Salida

Para cada número entero n desde 1 hacia arriba, genere el número de Hankelablen por nmatrices con entradas que son 0o 1.

Para n = 1,2,3,4,5las respuestas debe ser 2,12,230,12076,1446672. (Gracias a orlp por el código para producir estos).

Límite de tiempo

Ejecutaré su código en mi máquina y lo detendré después de 1 minuto. El código que genera las respuestas correctas hasta el mayor valor de n gana. El límite de tiempo es para todo, desde n = 1el valor más grande npara el que da una respuesta.

El ganador será la mejor respuesta para finales del sábado 18 de abril.

Desempate

En el caso de un empate por un máximo n, calcularé el tiempo que tarda en dar los resultados n+1y gana el más rápido. En el caso de que se ejecuten al mismo tiempo hasta dentro de un segundo n+1, la primera presentación gana.

Idiomas y bibliotecas

Puede usar cualquier idioma que tenga un compilador / intérprete / etc. para Linux y cualquier biblioteca que también esté disponible gratuitamente para Linux.

Mi maquina

Los tiempos se ejecutarán en mi máquina. Esta es una instalación estándar de ubuntu en un procesador AMD FX-8350 de ocho núcleos en una placa base Asus M5A78L-M / USB3 (Socket AM3 +, 8GB DDR3). Esto también significa que necesito poder ejecutar su código. Como consecuencia, solo use software gratuito fácilmente disponible e incluya instrucciones completas sobre cómo compilar y ejecutar su código.

Notas

Recomiendo no iterar sobre todas las matrices n por n e intentar detectar si cada una tiene la propiedad que describo. Primero, hay demasiados y segundo, parece que no hay una forma rápida de hacer esta detección .

Entradas principales hasta ahora

  • n = 8 por Peter Taylor. Java
  • n = 5 por orlp. Pitón
Comunidad
fuente
44
Realmente disfruto la palabra "Hankelable".
Alex A.
3
Por n=6el total es 260357434. Creo que la presión de la memoria es un problema mayor que el tiempo de CPU.
Peter Taylor
Esta es una pregunta asombrosa. He sido completamente francotirador.
alexander-brett

Respuestas:

7

Java (n = 8)

import java.util.*;
import java.util.concurrent.*;

public class HankelCombinatorics {
    public static final int NUM_THREADS = 8;

    private static final int[] FACT = new int[13];
    static {
        FACT[0] = 1;
        for (int i = 1; i < FACT.length; i++) FACT[i] = i * FACT[i-1];
    }

    public static void main(String[] args) {
        long prevElapsed = 0, start = System.nanoTime();
        for (int i = 1; i < 12; i++) {
            long count = count(i), elapsed = System.nanoTime() - start;
            System.out.format("%d in %dms, total elapsed %dms\n", count, (elapsed - prevElapsed) / 1000000, elapsed / 1000000);
            prevElapsed = elapsed;
        }
    }

    @SuppressWarnings("unchecked")
    private static long count(int n) {
        int[][] perms = new int[FACT[n]][];
        genPermsInner(0, 0, new int[n], perms, 0);

        // We partition by canonical representation of the row sum multiset, discarding any with a density > 50%.
        Map<CanonicalMatrix, Map<CanonicalMatrix, Integer>> part = new HashMap<CanonicalMatrix, Map<CanonicalMatrix, Integer>>();
        for (int m = 0; m < 1 << (2*n-1); m++) {
            int density = 0;
            int[] key = new int[n];
            for (int i = 0; i < n; i++) {
                key[i] = Integer.bitCount((m >> i) & ((1 << n) - 1));
                density += key[i];
            }
            if (2 * density <= n * n) {
                CanonicalMatrix _key = new CanonicalMatrix(key);
                Map<CanonicalMatrix, Integer> map = part.get(_key);
                if (map == null) part.put(_key, map = new HashMap<CanonicalMatrix, Integer>());
                map.put(new CanonicalMatrix(m, perms[0]), m);
            }
        }

        List<Job> jobs = new ArrayList<Job>();
        ExecutorService pool = Executors.newFixedThreadPool(NUM_THREADS);

        for (Map.Entry<CanonicalMatrix, Map<CanonicalMatrix, Integer>> e : part.entrySet()) {
            Job job = new Job(n, perms, e.getKey().sum() << 1 == n * n ? 0 : 1, e.getValue());
            jobs.add(job);
            pool.execute(job);
        }

        pool.shutdown();
        try {
            pool.awaitTermination(1, TimeUnit.DAYS); // i.e. until it's finished - inaccurate results are useless
        }
        catch (InterruptedException ie) {
            throw new IllegalStateException(ie);
        }

        long total = 0;
        for (Job job : jobs) total += job.subtotal;
        return total;
    }

    private static int genPermsInner(int idx, int usedMask, int[] a, int[][] perms, int off) {
        if (idx == a.length) perms[off++] = a.clone();
        else for (int i = 0; i < a.length; i++) {
            int m = 1 << (a[idx] = i);
            if ((usedMask & m) == 0) off = genPermsInner(idx+1, usedMask | m, a, perms, off);
        }
        return off;
    }

    static class Job implements Runnable {
        private volatile long subtotal = 0;
        private final int n;
        private final int[][] perms;
        private final int shift;
        private final Map<CanonicalMatrix, Integer> unseen;

        public Job(int n, int[][] perms, int shift, Map<CanonicalMatrix, Integer> unseen) {
            this.n = n;
            this.perms = perms;
            this.shift = shift;
            this.unseen = unseen;
        }

        public void run() {
            long result = 0;
            int[][] perms = this.perms;
            Map<CanonicalMatrix, Integer> unseen = this.unseen;
            while (!unseen.isEmpty()) {
                int m = unseen.values().iterator().next();
                Set<CanonicalMatrix> equiv = new HashSet<CanonicalMatrix>();
                for (int[] perm : perms) {
                    CanonicalMatrix canonical = new CanonicalMatrix(m, perm);
                    if (equiv.add(canonical)) {
                        result += canonical.weight() << shift;
                        unseen.remove(canonical);
                    }
                }
            }

            subtotal = result;
        }
    }

    static class CanonicalMatrix {
        private final int[] a;
        private final int hash;

        public CanonicalMatrix(int m, int[] r) {
            this(permuteRows(m, r));
        }

        public CanonicalMatrix(int[] a) {
            this.a = a;
            Arrays.sort(a);

            int h = 0;
            for (int i : a) h = h * 37 + i;
            hash = h;
        }

        private static int[] permuteRows(int m, int[] perm) {
            int[] cols = new int[perm.length];
            for (int i = 0; i < perm.length; i++) {
                for (int j = 0; j < cols.length; j++) cols[j] |= ((m >> (perm[i] + j)) & 1L) << i;
            }
            return cols;
        }

        public int sum() {
            int sum = 0;
            for (int i : a) sum += i;
            return sum;
        }

        public int weight() {
            int prev = -1, count = 0, weight = FACT[a.length];
            for (int col : a) {
                if (col == prev) weight /= ++count;
                else {
                    prev = col;
                    count = 1;
                }
            }
            return weight;
        }

        @Override public boolean equals(Object obj) {
            // Deliberately unsuitable for general-purpose use, but helps catch bugs faster.
            CanonicalMatrix that = (CanonicalMatrix)obj;
            for (int i = 0; i < a.length; i++) {
                if (a[i] != that.a[i]) return false;
            }
            return true;
        }

        @Override public int hashCode() {
            return hash;
        }
    }
}

Guardar como HankelCombinatorics.java, compilar como javac HankelCombinatorics.java, ejecutar como java -Xmx2G HankelCombinatorics.

Con NUM_THREADS = 4en mi máquina de cuatro núcleos se pone 20420819767436para n=8de 50 a 55 segundos han pasado, con una buena cantidad de variabilidad entre las carreras; Espero que pueda manejar fácilmente lo mismo en su máquina de ocho núcleos, pero tardará una hora o más en llegar n=9.

Cómo funciona

Dado n, hay matrices 2^(2n-1)binarias nx nHankel. Las filas se pueden permutar en n!formas, y las columnas en n!formas. Todo lo que tenemos que hacer es evitar el doble conteo ...

Si calcula la suma de cada fila, entonces ni permutar las filas ni permutar las columnas cambia el conjunto múltiple de sumas. P.ej

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

tiene multiset de suma de filas {3, 3, 2, 2, 2}, y también lo hacen todas las matrices Hankelable derivadas de él. Esto significa que podemos agrupar las matrices de Hankel por estos conjuntos de suma de filas y luego manejar cada grupo de forma independiente, explotando múltiples núcleos de procesador.

También hay una simetría explotable: las matrices con más ceros que unos están en biyección con las matrices con más unos que ceros.

Doble conteo se produce cuando matriz de Hankel M_1con permutación fila r_1y la permutación de columna c_1coincide con matriz de Hankel M_2con permutación fila r_2y la permutación de columna c_2(con un máximo de dos, pero no los tres M_1 = M_2, r_1 = r_2, c_1 = c_2). Las permutaciones de fila y columna son independientes, por lo que si aplicamos permutación r_1de M_1fila y permutación de fila r_2a M_2, las columnas como multisets deben ser iguales. Entonces, para cada grupo, calculo todos los conjuntos de columnas múltiples obtenidos aplicando una permutación de fila a una matriz en el grupo. La manera fácil de obtener una representación canónica de los conjuntos múltiples es ordenar las columnas, y eso también es útil en el siguiente paso.

Una vez obtenidos los distintos conjuntos de columnas, necesitamos encontrar cuántas n!permutaciones de cada uno son únicos. En este punto, el recuento doble solo puede ocurrir si un conjunto múltiple de columnas dado tiene columnas duplicadas: lo que debemos hacer es contar el número de ocurrencias de cada columna distinta en el conjunto múltiple y luego calcular el coeficiente multinomial correspondiente. Como las columnas están ordenadas, es fácil hacer el recuento.

Finalmente los sumamos todos.

La complejidad asintótica no es trivial para calcular con total precisión, porque necesitamos hacer algunas suposiciones sobre los conjuntos. Evaluamos el orden de los conjuntos de 2^(2n-2) n!columnas múltiples, tomando n^2 ln ntiempo para cada uno (incluida la clasificación); Si la agrupación no toma más que un ln nfactor, tenemos complejidad de tiempo Theta(4^n n! n^2 ln n). Pero como los factores exponenciales dominan por completo a los polinomios, lo es Theta(4^n n!) = Theta((4n/e)^n).

Peter Taylor
fuente
Esto es muy impresionante ¿Podría decir algo sobre el algoritmo que ha utilizado?
3

Python2 / 3

Enfoque bastante ingenuo, en un lenguaje lento:

import itertools

def permute_rows(m):
    for perm in itertools.permutations(m):
        yield perm

def permute_columns(m):
    T = zip(*m)
    for perm in itertools.permutations(T):
        yield zip(*perm)

N = 1
while True:
    base_template = ["abcdefghijklmnopqrstuvwxyz"[i:i+N] for i in range(N)]

    templates = set()
    for c in permute_rows(base_template):
        for m in permute_columns(c):
            templates.add("".join("".join(row) for row in m))

    def possibs(free, templates):
        if free == 2*N - 1:
            return set(int(t, 2) for t in templates)

        s = set()
        for b in "01":
            new_templates = set(t.replace("abcdefghijklmnopqrstuvwxyz"[free], b) for t in templates)
            s |= possibs(free + 1, new_templates)

        return s

    print(len(possibs(0, templates)))
    N += 1

Ejecutar escribiendo python script.py.

orlp
fuente
Tiene el lenguaje listado como Python 2/3, pero para que funcione en Python 2, ¿no lo necesita from __future__ import print_function(o algo así)?
Alex A.
2
@AlexA. Normalmente sí, pero no en este caso. Considere el comportamiento de Python2 cuando escriba return(1). Ahora reemplace returncon print:)
orlp
¡Frio! Aprendo algo nuevo cada día. :)
Alex A.
2

Haskell

import Data.List
import Data.Hashable
import Control.Parallel.Strategies
import Control.Parallel
import qualified Data.HashSet as S

main = mapM putStrLn $ map (show.countHankellable) [1..]

a§b=[a!!i|i<-b]

hashNub :: (Hashable a, Eq a) => [a] -> [a]
hashNub l = go S.empty l
    where
      go _ []     = []
      go s (x:xs) = if x `S.member` s then go s xs
                                    else x : go (S.insert x s) xs

pmap = parMap rseq

makeMatrix :: Int->[Bool]->[[Bool]]
makeMatrix n vars = [vars§[i..i+n-1]|i<-[0..n-1]]

countHankellable :: Int -> Int
countHankellable n = let
    s = permutations [0..n-1]
    conjugates m = concat[permutations[r§q|r<-m]|q<-s]
    variableSets = sequence [[True,False]|x<-[0..2*(n-1)]]
 in
    length.hashNub.concat.pmap (conjugates.makeMatrix n ) $ variableSets

No es tan rápido como el de Peter: ¡es una configuración bastante impresionante que tiene allí! Ahora con mucho más código copiado de internet. Uso:

$ ghc -threaded hankell.hs
$ ./hankell
alexander-brett
fuente
Una respuesta de Haskell siempre es bienvenida. Gracias.
@Lembik: ¿cómo va el mío en tu máquina?
alexander-brett