Comptez le nombre de matrices Hankelable

12

Contexte

Une matrice de Hankel binaire est une matrice à diagonales asymétriques constantes (diagonales à pente positive) contenant uniquement 0s et 1s. Par exemple, une matrice de Hankel binaire 5x5 ressemble

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

a, b, c, d, e, f, g, h, isont soit 0ou 1.

Définissons une matrice M comme Hankelable s'il y a une permutation de l'ordre des lignes et des colonnes de M afin que M soit une matrice de Hankel. Cela signifie que l'on peut appliquer une permutation à l'ordre des lignes et éventuellement une autre aux colonnes.

Le défi

Le défi est de compter le nombre de Hankelable n par nmatrices pour nune valeur aussi élevée que possible.

Production

Pour chaque entier n à partir de 1, affichez le nombre de Hankelablen par des nmatrices avec des entrées qui sont 0ou 1.

Car n = 1,2,3,4,5les réponses devraient être 2,12,230,12076,1446672. (Merci à orlp pour le code pour les produire.)

Limite de temps

Je vais exécuter votre code sur ma machine et l'arrêter après 1 minute. Le code qui génère les bonnes réponses jusqu'à la plus grande valeur de n victoires. Le délai est pour tout, de n = 1la plus grande valeur npour laquelle vous donnez une réponse.

Le gagnant sera la meilleure réponse d'ici la fin du samedi 18 avril.

Tie break

Dans le cas d'une égalité pour un maximum, nje chronomètrerai combien de temps il faut pour abandonner les sorties n+1et la plus rapide gagne. Dans le cas où ils s'exécutent en même temps dans une seconde jusqu'à n+1, la première soumission gagne.

Langues et bibliothèques

Vous pouvez utiliser n'importe quelle langue disposant d'un compilateur / interprète / etc disponible gratuitement. pour Linux et toutes les bibliothèques qui sont également disponibles gratuitement pour Linux.

Ma machine

Les horaires seront exécutés sur ma machine. Il s'agit d'une installation ubuntu standard sur un processeur AMD FX-8350 à huit cœurs sur une carte mère Asus M5A78L-M / USB3 (Socket AM3 +, 8 Go DDR3). Cela signifie également que je dois pouvoir exécuter votre code. Par conséquent, n'utilisez que des logiciels gratuits facilement disponibles et veuillez inclure des instructions complètes sur la façon de compiler et d'exécuter votre code.

Remarques

Je recommande de ne pas itérer sur toutes les matrices n par n et d'essayer de détecter si chacune a la propriété que je décris. Premièrement, il y en a trop et deuxièmement, il ne semble pas y avoir de moyen rapide de faire cette détection .

Principales entrées à ce jour

  • n = 8 par Peter Taylor. Java
  • n = 5 par orlp. Python
Communauté
la source
4
J'aime vraiment le mot «Hankelable».
Alex A.
3
Pour n=6le total, c'est 260357434. Je pense que la pression de la mémoire est un problème plus important que le temps CPU.
Peter Taylor
C'est une question géniale. J'ai été complètement snipé.
alexander-brett

Réponses:

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;
        }
    }
}

Enregistrer sous HankelCombinatorics.java, compiler sous javac HankelCombinatorics.java, exécuter sous java -Xmx2G HankelCombinatorics.

Avec NUM_THREADS = 4ma machine quad-core , il obtient 20420819767436pour n=8dans 50 à 55 secondes se sont écoulées, avec une quantité de juste de la variabilité entre les pistes; J'espère qu'il devrait facilement gérer la même chose sur votre machine octa-core, mais cela prendra une heure ou plus pour l'obtenir n=9.

Comment ça fonctionne

Étant donné n, il existe des matrices 2^(2n-1)binaires nx nHankel. Les lignes peuvent être permutées de n!différentes manières et les colonnes de n!différentes manières. Tout ce que nous devons faire est d'éviter le double comptage ...

Si vous calculez la somme de chaque ligne, alors ni permuter les lignes ni permuter les colonnes ne modifie le multi-ensemble de sommes. Par exemple

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

possède un multi-ensemble de somme de lignes {3, 3, 2, 2, 2}, tout comme toutes les matrices Hankelable qui en dérivent. Cela signifie que nous pouvons regrouper les matrices Hankel par ces multisets de somme de lignes, puis gérer chaque groupe indépendamment, en exploitant plusieurs cœurs de processeur.

Il y a aussi une symétrie exploitable: les matrices avec plus de zéros que de uns sont en bijection avec les matrices avec plus de zéros.

Double comptage se produit lorsque la matrice Hankel M_1avec permutation de ligne r_1et permutation de colonne c_1correspond à la matrice Hankel M_2avec permutation de ligne r_2et permutation de colonne c_2(avec un maximum de deux , mais pas tous les trois M_1 = M_2, r_1 = r_2, c_1 = c_2). Les permutations de ligne et de colonne sont indépendantes, donc si nous appliquons la permutation de ligne r_1à M_1et la permutation de ligne r_2à M_2, les colonnes en tant que multisets doivent être égales. Donc, pour chaque groupe, je calcule tous les multisets de colonne obtenus en appliquant une permutation de ligne à une matrice du groupe. Le moyen facile d'obtenir une représentation canonique des multisets est de trier les colonnes, ce qui est également utile à l'étape suivante.

Après avoir obtenu les multisets de colonnes distinctes, nous devons trouver combien de n!permutations de chacune sont uniques. À ce stade, le double comptage ne peut se produire que si un multiset de colonne donné a des colonnes en double: ce que nous devons faire est de compter le nombre d'occurrences de chaque colonne distincte dans le multiset, puis de calculer le coefficient multinomial correspondant. Comme les colonnes sont triées, il est facile de faire le décompte.

Enfin, nous les additionnons tous.

La complexité asymptotique n'est pas triviale à calculer avec une précision totale, car nous devons faire quelques hypothèses sur les ensembles. Nous évaluons sur l'ordre des 2^(2n-2) n!multisets de colonne, en prenant du n^2 ln ntemps pour chacun (y compris le tri); si le regroupement ne prend pas plus d'un ln nfacteur, nous avons une complexité temporelle Theta(4^n n! n^2 ln n). Mais puisque les facteurs exponentiels dominent complètement les polynômes, c'est le cas Theta(4^n n!) = Theta((4n/e)^n).

Peter Taylor
la source
C'est très impressionnant. Pourriez-vous dire quelque chose sur l'algorithme que vous avez utilisé?
3

Python2 / 3

Approche assez naïve, dans un langage lent:

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

Exécutez en tapant python script.py.

orlp
la source
Vous avez le langage répertorié comme Python 2/3, mais pour qu'il fonctionne en Python 2, n'avez-vous pas besoin from __future__ import print_function(ou quelque chose comme ça)?
Alex A.
2
@AlexA. Normalement, oui, mais pas dans ce cas. Tenez compte du comportement de Python2 lorsque vous tapez return(1). Remplacez maintenant returnpar print:)
orlp
Cool! J'apprends quelque chose de nouveau tous les jours. :)
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

Nulle part aussi vite que Peter - c'est une configuration assez impressionnante qu'il a là! Maintenant, avec beaucoup plus de code copié à partir d'Internet. Usage:

$ ghc -threaded hankell.hs
$ ./hankell
alexander-brett
la source
Une réponse Haskell est toujours la bienvenue. Je vous remercie.
@Lembik - comment va le mien sur votre machine?
alexander-brett