Calcule el Hafnian lo más rápido posible

12

El desafío es escribir el código más rápido posible para calcular el hafniano de una matriz .

El hafniano de una matriz simétrica 2n-por- se define como:2nA

Aquí S 2n representa el conjunto de todas las permutaciones de los enteros de 1a 2n, es decir [1, 2n].

El enlace de wikipedia también ofrece una fórmula de aspecto diferente que puede ser de interés (e incluso existen métodos más rápidos si busca más en la web). La misma página wiki habla de matrices de adyacencia, pero su código también debería funcionar para otras matrices. Puede suponer que todos los valores serán enteros, pero no que sean todos positivos.

También hay un algoritmo más rápido, pero parece difícil de entender. y Christian Sievers fue el primero en implementarlo (en Haskell).

En esta pregunta, todas las matrices son cuadradas y simétricas con una dimensión par.

Implementación de referencia (tenga en cuenta que esto está utilizando el método más lento posible).

Aquí hay un ejemplo de código de Python del Sr. Xcoder.

from itertools import permutations
from math import factorial

def hafnian(matrix):
    my_sum = 0
    n = len(matrix) // 2
    for sigma in permutations(range(n*2)):
        prod = 1
        for j in range(n):
            prod *= matrix[sigma[2*j]][sigma[2*j+1]]
        my_sum += prod
    return my_sum / (factorial(n) * 2 ** n)

print(hafnian([[-1, 1, 1, -1, 0, 0, 1, -1], [1, 0, 1, 0, -1, 0, -1, -1], [1, 1, -1, 1, -1, -1, 0, -1], [-1, 0, 1, -1, -1, 1, -1, 0], [0, -1, -1, -1, -1, 0, 0, -1], [0, 0, -1, 1, 0, 0, 1, 1], [1, -1, 0, -1, 0, 1, 1, 0], [-1, -1, -1, 0, -1, 1, 0, 1]]))
4

M = [[1, 1, 0, 0, 0, 0, 0, 1, 0, 0], [1, 1, -1, 0, -1, 1, 1, 1, 0, -1], [0, -1, -1, -1, 0, -1, -1, 0, -1, 1], [0, 0, -1, 1, -1, 1, -1, 0, 1, -1], [0, -1, 0, -1, -1, -1, -1, 1, -1, 1], [0, 1, -1, 1, -1, 1, -1, -1, 1, -1], [0, 1, -1, -1, -1, -1, 1, 0, 0, 0], [1, 1, 0, 0, 1, -1, 0, 1, 1, -1], [0, 0, -1, 1, -1, 1, 0, 1, 1, 1], [0, -1, 1, -1, 1, -1, 0, -1, 1, 1]]

print(hafnian(M))
-13

M = [[-1, 0, -1, -1, 0, -1, 0, 1, -1, 0, 0, 0], [0, 0, 0, 0, 0, -1, 0, 1, -1, -1, -1, -1], [-1, 0, 0, 1, 0, 0, 0, 1, -1, 1, -1, 0], [-1, 0, 1, -1, 1, -1, -1, -1, 0, -1, -1, -1], [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, -1, 0], [-1, -1, 0, -1, 0, 0, 1, 1, 1, 1, 1, 0], [0, 0, 0, -1, 0, 1, 1, -1, -1, 0, 1, 0], [1, 1, 1, -1, 0, 1, -1, 1, -1, -1, -1, -1], [-1, -1, -1, 0, 0, 1, -1, -1, -1, 1, -1, 0], [0, -1, 1, -1, 1, 1, 0, -1, 1, -1, 1, 1], [0, -1, -1, -1, -1, 1, 1, -1, -1, 1, 0, -1], [0, -1, 0, -1, 0, 0, 0, -1, 0, 1, -1, 1]]

print(hafnian(M))
13

M = [[-1, 1, 0, 1, 0, -1, 0, 0, -1, 1, -1, 1, 0, -1], [1, -1, 1, -1, 1, 1, -1, 0, -1, 1, 1, 0, 0, -1], [0, 1, 1, 1, -1, 1, -1, -1, 0, 0, -1, 0, -1, -1], [1, -1, 1, -1, 1, 0, 1, 1, -1, -1, 0, 0, 1, 1], [0, 1, -1, 1, 0, 1, 0, 1, -1, -1, 1, 1, 0, -1], [-1, 1, 1, 0, 1, 1, -1, 0, 1, -1, -1, -1, 1, -1], [0, -1, -1, 1, 0, -1, -1, -1, 0, 1, -1, 0, 1, -1], [0, 0, -1, 1, 1, 0, -1, 0, 0, -1, 0, 0, 0, 1], [-1, -1, 0, -1, -1, 1, 0, 0, 1, 1, 0, 1, -1, 0], [1, 1, 0, -1, -1, -1, 1, -1, 1, 1, 1, 0, 1, 0], [-1, 1, -1, 0, 1, -1, -1, 0, 0, 1, -1, 0, -1, 0], [1, 0, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 1], [0, 0, -1, 1, 0, 1, 1, 0, -1, 1, -1, 1, 1, -1], [-1, -1, -1, 1, -1, -1, -1, 1, 0, 0, 0, 1, -1, -1]]

print(hafnian(M))
83

La tarea

Debería escribir un código que, dado 2npor una 2nmatriz, genere su hafniano.

Como tendré que probar su código, sería útil si pudiera darme una forma simple de dar una matriz como entrada a su código, por ejemplo, leyendo desde el estándar. Probaré su código en matrices elegidas al azar con elementos seleccionado de {-1, 0, 1}. El propósito de pruebas como esta es reducir la posibilidad de que el Hafnian sea un valor muy grande.

Idealmente, su código podrá leerse en matrices exactamente como las tengo en los ejemplos de esta pregunta directamente desde el estándar. Así es como se vería la entrada, [[1,-1],[-1,-1]]por ejemplo. Si desea utilizar otro formato de entrada, solicítelo y haré todo lo posible para adaptarlo.

Puntuaciones y lazos

Probaré su código en matrices aleatorias de tamaño creciente y me detendré la primera vez que su código tarde más de 1 minuto en mi computadora. Las matrices de puntuación serán consistentes para todas las presentaciones a fin de garantizar la equidad.

Si dos personas obtienen el mismo puntaje, entonces el ganador es el que tiene el valor más rápido n. Si están dentro de 1 segundo el uno del otro, entonces es el publicado primero.

Idiomas y bibliotecas

Puede usar cualquier idioma y bibliotecas disponibles que desee, pero ninguna función preexistente para calcular el hafniano. Siempre que sea posible, sería bueno poder ejecutar su código, así que incluya una explicación completa sobre cómo ejecutar / compilar su código en Linux, si es posible.

Mi máquina Los tiempos se ejecutarán en mi máquina de 64 bits. Esta es una instalación estándar de ubuntu con 8 GB de RAM, procesador AMD FX-8350 de ocho núcleos y Radeon HD 4250. Esto también significa que necesito poder ejecutar su código.

Solicite respuestas en más idiomas.

Sería genial obtener respuestas en su lenguaje de programación súper rápido favorito. Para empezar, ¿qué tal fortran , nim y óxido ?

Tabla de clasificación

  • 52 por millas usando C ++ . 30 segundos.
  • 50 por ngn usando C . 50 segundos
  • 46 por Christian Sievers usando Haskell . 40 segundos
  • 40 por millas usando Python 2 + pypy . 41 segundos
  • 34 por ngn usando Python 3 + pypy . 29 segundos
  • 28 por Dennis usando Python 3 . 35 segundos (Pypy es más lento)

fuente
¿Existe un límite para los valores absolutos de las entradas de la matriz? ¿Podemos devolver una aproximación de coma flotante? ¿Tenemos que usar enteros de precisión arbitrarios?
Dennis
@Dennis En la práctica solo usaré -1,0,1 para probar (elegido al azar). No quiero que sea un gran desafío internacional. Honestamente, no sé si alcanzaremos los límites de las entradas de 64 bits antes de que el código sea demasiado lento para ejecutarse, pero supongo que no lo haremos. Actualmente no estamos cerca de eso.
Si las entradas están limitadas a -1,0,1 , eso debería mencionarse en la pregunta. ¿Nuestro código tiene que funcionar para otras matrices?
Dennis
@ Dennis Una versión antigua solía decir eso, pero debo haber escrito sobre ella. Preferiría que el código no estuviera especializado para -1,0,1 entradas, pero supongo que no puedo detenerlo.
¿Tienes más casos de prueba? tal vez para mayor n ?
millas

Respuestas:

14

Haskell

import Control.Parallel.Strategies
import qualified Data.Vector.Unboxed as V
import qualified Data.Vector as VB

type Poly = V.Vector Int

type Matrix = VB.Vector ( VB.Vector Poly )

constpoly :: Int -> Int -> Poly
constpoly n c = V.generate (n+1) (\i -> if i==0 then c else 0)

add :: Poly -> Poly -> Poly
add = V.zipWith (+)

shiftmult :: Poly -> Poly -> Poly
shiftmult a b = V.generate (V.length a) 
                           (\i -> sum [ a!j * b!(i-1-j) | j<-[0..i-1] ])
  where (!) = V.unsafeIndex

x :: Matrix -> Int -> Int -> Int -> Poly -> Int
x  _    0  _ m p = m * V.last p
x mat n c m p =
  let mat' = VB.generate (2*n-2) $ \i ->
             VB.generate i       $ \j ->
                 shiftmult (mat!(2*n-1)!i) (mat!(2*n-2)!j) `add`
                 shiftmult (mat!(2*n-1)!j) (mat!(2*n-2)!i) `add`
                 (mat!i!j)
      p' = p `add` shiftmult (mat!(2*n-1)!(2*n-2)) p
      (!) = VB.unsafeIndex
      r = if c>0 then parTuple2 rseq rseq else r0
      (a,b) = (x mat (n-1) (c-1) m p, x mat' (n-1) (c-1) (-m) p')
              `using` r
  in a+b

haf :: [[Int]] -> Int
haf m = let n=length m `div` 2
        in x (VB.fromList $ map (VB.fromList . map (constpoly n)) m) 
             n  5  ((-1)^n)  (constpoly n 1) 

main = getContents >>= print . haf . read

Esto implementa una variación del Algoritmo 2 de Andreas Björklund: Contando combinaciones perfectas tan rápido como Ryser .

Compile usando ghcopciones de tiempo de compilación -O3 -threadedy use opciones de tiempo de ejecución +RTS -Npara la paralelización. Toma entrada de stdin.

Christian Sievers
fuente
2
Tal vez tenga en cuenta que parallely vectordebe ser instalado?
H.PWiz
@ H.PWiz Nadie se quejó aquí , pero claro, señalando que no hará daño. Bueno, ahora lo hiciste.
Christian Sievers
@ChristianSievers No creo que se quejen. Es posible que el OP no esté familiarizado con Haskell, por lo que es una buena idea indicar explícitamente qué se debe instalar para poder sincronizar el código.
Dennis
@ Dennis No quise decir "te quejaste" sino "lo notaste". Y no pensé en quejarme como algo negativo. El OP es el mismo que en el desafío al que me vinculé, así que no debería haber ningún problema.
Christian Sievers
N = 40 en 7.5 segundos en TIO ... ¡Hombre, esto es rápido!
Dennis
6

Python 3

from functools import lru_cache

@lru_cache(maxsize = None)
def haf(matrix):
	n = len(matrix)
	if n == 2: return matrix[0][1]
	h = 0
	for j in range(1, n):
		if matrix[0][j] == 0: continue
		copy = list(matrix)
		del copy[:j+1:j]
		copy = list(zip(*copy))
		del copy[:j+1:j]
		h += matrix[0][j] * haf(tuple(copy))
	return h

print(haf(tuple(map(tuple, eval(open(0).read())))))

Pruébalo en línea!

Dennis
fuente
6

C ++ (gcc)

#define T(x) ((x)*((x)-1)/2)
#define S 1
#define J (1<<S)
#define TYPE int

#include <iostream>
#include <vector>
#include <string>
#include <pthread.h>

using namespace std;

struct H {
    int s, w, t;
    TYPE *b, *g;
};

void *solve(void *a);
void hafnian(TYPE *b, int s, TYPE *g, int w, int t);

int n, m, ti = 0;
TYPE r[J] = {0};
pthread_t pool[J];

int main(void) {
    vector<int> a;
    string s;
    getline(cin, s);

    for (int i = 0; i < s.size(); i++)
        if (s[i] == '0' || s[i] == '1')
            a.push_back((s[i-1] == '-' ? -1 : 1)*(s[i] - '0'));

    for (n = 1; 4*n*n < a.size(); n++);
    m = n+1;

    TYPE z[T(2*n)*m] = {0}, g[m] = {0};

    for (int j = 1; j < 2*n; j++)
        for (int k = 0; k < j; k++)
            z[(T(j)+k)*m] = a[j*2*n+k];
    g[0] = 1;

    hafnian(z, 2*n, g, 1, -1);

    TYPE h = 0;
    for (int t = 0; t < ti; t++) {
        pthread_join(pool[t], NULL);
        h += r[t];
    }

    cout << h << endl;

    return 0;
}

void *solve(void *a) {
    H *p = reinterpret_cast<H*>(a);
    hafnian(p->b, p->s, p->g, p->w, p->t);
    delete[] p->b;
    delete[] p->g;
    delete p;
    return NULL;
}

void hafnian(TYPE *b, int s, TYPE *g, int w, int t) {
    if (t == -1 && (n < S || s/2 == n-S)) {
        H *p = new H;
        TYPE *c = new TYPE[T(s)*m], *e = new TYPE[m];
        copy(b, b+T(s)*m, c);
        copy(g, g+m, e);
        p->b = c;
        p->s = s;
        p->g = e;
        p->w = w;
        p->t = ti;
        pthread_create(pool+ti, NULL, solve, p);
        ti++;
    }
    else if (s > 0) {
        TYPE c[T(s-2)*m], e[m];
        copy(b, b+T(s-2)*m, c);
        hafnian(c, s-2, g, -w, t);
        copy(g, g+m, e);

        for (int u = 0; u < n; u++) {
            TYPE *d = e+u+1,
                  p = g[u], *x = b+(T(s)-1)*m;
            for (int v = 0; v < n-u; v++)
                d[v] += p*x[v];
        }

        for (int j = 1; j < s-2; j++)
            for (int k = 0; k < j; k++)
                for (int u = 0; u < n; u++) {
                    TYPE *d = c+(T(j)+k)*m+u+1,
                          p = b[(T(s-2)+j)*m+u], *x = b+(T(s-1)+k)*m,
                          q = b[(T(s-2)+k)*m+u], *y = b+(T(s-1)+j)*m;
                    for (int v = 0; v < n-u; v++)
                        d[v] += p*x[v] + q*y[v];
                }

        hafnian(c, s-2, e, w, t);
    }
    else
        r[t] += w*g[n];
}

Pruébalo en línea! (13s para n = 24)

Basado en la implementación más rápida de Python en mi otra publicación. Edite la segunda línea #define S 3en su máquina de 8 núcleos y compile con g++ -pthread -march=native -O2 -ftree-vectorize.

Divide el trabajo por la mitad, por lo que el valor de Sdebería ser log2(#threads). Los tipos se pueden cambiar fácilmente entre int, long, floaty doublemodificando el valor de #define TYPE.

millas
fuente
Esta es la respuesta principal hasta ahora. Su código en realidad no lee en la entrada como se especifica, ya que no puede hacer frente a espacios. Tuve que hacer, por ejemplotr -d \ < matrix52.txt > matrix52s.txt
@Lembik Lo siento, solo lo usé contra la matriz sin espacio de tamaño 24. Sin embargo, ahora se corrigió para trabajar con espacios.
millas
4

Python 3

Esto calcula haf (A) como una suma memorizada (A [i] [j] * haf (A sin filas y cols i y j)).

#!/usr/bin/env python3
import json,sys
a=json.loads(sys.stdin.read())
n=len(a)//2
b={0:1}
def haf(x):
 if x not in b:
  i=0
  while not x&(1<<i):i+=1
  x1=x&~(1<<i)
  b[x]=sum(a[i][j]*haf(x1&~(1<<j))for j in range(2*n)if x1&(1<<j)and a[i][j])
 return b[x]
print(haf((1<<2*n)-1))
ngn
fuente
3

C

Otra implicación del artículo de Andreas Björklund , que es mucho más fácil de entender si también observa el código Haskell de Christian Sievers . Para los primeros niveles de la recursividad, distribuye subprocesos round-robin sobre las CPU disponibles. El último nivel de la recursión, que representa la mitad de las invocaciones, se optimiza a mano.

Compilar con: gcc -O3 -pthread -march=native; gracias @Dennis por una aceleración 2x

n = 24 en 24s en TIO

#define _GNU_SOURCE
#include<sched.h>
#include<stdio.h>
#include<stdlib.h>
#include<memory.h>
#include<unistd.h>
#include<pthread.h>
#define W while
#define R return
#define S static
#define U (1<<31)
#define T(i)((i)*((i)-1)/2)
typedef int I;typedef long L;typedef char C;typedef void V;
I n,ncpu,icpu;
S V f(I*x,I*y,I*z){I i=n,*z1=z+n;W(i){I s=0,*x2=x,*y2=y+--i;W(y2>=y)s+=*x2++**y2--;*z1--+=s;}}
typedef struct{I m;V*a;V*p;pthread_barrier_t*bar;I r;}A;S V*(h1)(V*);
I h(I m,I a[][n+1],I*p){
 m-=2;I i,j,k=0,u=T(m),v=u+m,b[u][n+1],q[n+1];
 if(!m){I*x=a[v+m],*y=p+n-1,s=0;W(y>=p)s-=*x++**y--;R s;}
 memcpy(b,a,sizeof(b));memcpy(q,p,sizeof(q));f(a[v+m],p,q);
 for(i=1;i<m;i++)for(j=0;j<i;j++){f(a[u+i],a[v+j],b[k]);f(a[u+j],a[v+i],b[k]);k++;}
 if(2*n-m>8)R h(m,a,p)-h(m,b,q);
 pthread_barrier_t bar;pthread_barrier_init(&bar,0,2);pthread_t th;
 cpu_set_t cpus;CPU_ZERO(&cpus);CPU_SET(icpu++%ncpu,&cpus);
 pthread_attr_t attr;pthread_attr_init(&attr);
 pthread_attr_setaffinity_np(&attr,sizeof(cpu_set_t),&cpus);
 A arg={m,a,p,&bar};pthread_create(&th,&attr,h1,&arg);
 I r=h(m,b,q);pthread_barrier_wait(&bar);pthread_join(th,0);pthread_barrier_destroy(&bar);
 R arg.r-r;
}
S V*h1(V*x0){A*x=(A*)x0;x->r=h(x->m,x->a,x->p);pthread_barrier_wait(x->bar);R 0;}
I main(){
 ncpu=sysconf(_SC_NPROCESSORS_ONLN);
 S C s[200000];I i=0,j=0,k,l=0;W((k=read(0,s+l,sizeof(s)-l))>0)l+=k;
 n=1;W(s[i]!=']')n+=s[i++]==',';n/=2;
 I a[T(2*n)][n+1];memset(a,0,sizeof(a));k=0;
 for(i=0;i<2*n;i++)for(j=0;j<2*n;j++){
  W(s[k]!='-'&&(s[k]<'0'||s[k]>'9'))k++;
  I v=0,m=s[k]=='-';k+=m;W(k<l&&('0'<=s[k]&&s[k]<='9'))v=10*v+s[k++]-'0';
  if(i>j)*a[T(i)+j]=v*(1-2*m);
 }
 I p[n+1];memset(p,0,sizeof(p));*p=1;
 printf("%d\n",(1-2*(n&1))*h(2*n,a,p));
 R 0;
}

Algoritmo:

La matriz, que es simétrica, se almacena en forma triangular inferior izquierda. Los índices de triángulo i,jcorresponden al índice lineal T(max(i,j))+min(i,j)donde Tes una macro para i*(i-1)/2. Los elementos de la matriz son polinomios de grado n. Un polinomio se representa como una matriz de coeficientes ordenados desde el término constante ( p[0]) hasta el coeficiente de x n ( p[n]). Los valores iniciales de la matriz -1,0,1 se convierten primero en polinomios constantes.

Realizamos un paso recursivo con dos argumentos: la media matriz (es decir, el triángulo) ade polinomios y un polinomio separado p(denominado beta en el documento). Reducimos el mproblema de tamaño (inicialmente m=2*n) recursivamente a dos problemas de tamaño m-2y devolvemos la diferencia de sus hafnianos. Una de ellas es usar lo mismo asin sus dos últimas filas, y lo mismo p. Otra es usar el triángulo b[i][j] = a[i][j] + shmul(a[m-1][i],a[m-2][j]) + shmul(a[m-1][j],a[m-2][i])(donde shmulestá la operación shift-multiply en polinomios; es como un producto polinomial como de costumbre, adicionalmente multiplicado por la variable "x"; las potencias superiores a x ^ n son ignoradas), y el polinomio separado q = p + shmul(p,a[m-1][m-2]). Cuando la repetición realiza un tamaño-0 a, volvemos el mayor coeficiente de p: p[n].

La operación shift-and-multiply se implementa en función f(x,y,z). Se modifica zen el lugar. Hablando libremente, lo hace z += shmul(x,y). Esta parece ser la parte más crítica para el rendimiento.

Una vez que la recursión ha finalizado, debemos corregir el signo del resultado multiplicando por (-1) n .

ngn
fuente
¿Podría mostrar un ejemplo explícito de la entrada que acepta su código? Digamos para una matriz de 2 por 2. Además, ¡parece que has jugado tu código! (Este es un desafío de código más rápido, no un desafío de golf.)
@Lembik Para el registro, como dije en el chat, la entrada está en el mismo formato que los ejemplos: json (en realidad, solo lee los números y usa n = sqrt (len (input)) / 2). Por lo general, escribo un código corto, incluso cuando el golf no es un requisito.
ngn
¿Cuál es la matriz de mayor tamaño que debe admitir este nuevo código?
1
-march=nativeHará una gran diferencia aquí. Al menos en TIO, casi reduce el tiempo de la pared a la mitad.
Dennis
1
Además, al menos en TIO, el ejecutable producido por gcc será aún más rápido.
Dennis
3

Pitón

Esta es una implementación de referencia directa del Algoritmo 2 del documento mencionado . Los únicos cambios fueron mantener solo el valor actual de B , eliminando los valores de β solo actualizando g cuando iX , y multiplicando el polinomio truncado calculando solo los valores hasta el grado n .

from itertools import chain,combinations

def powerset(s):
    return chain.from_iterable(combinations(s, k) for k in range(len(s)+1))

def padd(a, b):
    return [a[i]+b[i] for i in range(len(a))]

def pmul(a, b):
    n = len(a)
    c = [0]*n
    for i in range(n):
        for j in range(n):
            if i+j < n:
                c[i+j] += a[i]*b[j]
    return c

def hafnian(m):
    n = len(m) / 2
    z = [[[c]+[0]*n for c in r] for r in m]
    h = 0
    for x in powerset(range(1, n+1)):
        b = z
        g = [1] + [0]*n
        for i in range(1, n+1):
            if i in x:
                g = pmul(g, [1] + b[0][1][:n])
                b = [[padd(b[j+2][k+2], [0] + padd(pmul(b[0][j+2], b[1][k+2]), pmul(b[0][k+2], b[1][j+2]))[:n]) if j != k else 0 for k in range(2*n-2*i)] for j in range(2*n-2*i)]
            else:
                b = [r[2:] for r in b[2:]]
        h += (-1)**(n - len(x)) * g[n]
    return h

Pruébalo en línea!

Aquí hay una versión más rápida con algunas de las optimizaciones fáciles.

def hafnian(m):
  n = len(m)/2
  z = [[0]*(n+1) for _ in range(n*(2*n-1))]
  for j in range(1, 2*n):
    for k in range(j):
      z[j*(j-1)/2+k][0] = m[j][k]
  return solve(z, 2*n, 1, [1] + [0]*n, n)

def solve(b, s, w, g, n):
  if s == 0:
    return w*g[n]
  c = [b[(j+1)*(j+2)/2+k+2][:] for j in range(1, s-2) for k in range(j)]
  h = solve(c, s-2, -w, g, n)
  e = g[:]
  for u in range(n):
    for v in range(n-u):
      e[u+v+1] += g[u]*b[0][v]
  for j in range(1, s-2):
    for k in range(j):
      for u in range(n):
        for v in range(n-u):
          c[j*(j-1)/2+k][u+v+1] += b[(j+1)*(j+2)/2][u]*b[(k+1)*(k+2)/2+1][v] + b[(k+1)*(k+2)/2][u]*b[(j+1)*(j+2)/2+1][v]
  return h + solve(c, s-2, w, e, n)

Pruébalo en línea!

Para mayor diversión, aquí hay una implementación de referencia en J.

millas
fuente
Esto es bastante lento de todas las comprensiones de la lista y de calcular valores equivalentes en la diagonal, por lo que no hay necesidad de comparar esto.
millas
¡Bastante impresionante!
¡Muy agradable! Intenté algo similar con sympy, que fue sorprendentemente lento y, aunque fue correcto para los pequeños ejemplos, arrojó, después de mucho tiempo, un resultado incorrecto para la matriz 24 * 24. No tengo idea de lo que está pasando allí. - Según la descripción anterior del Algoritmo 2, la multiplicación polinómica ya está destinada a ser truncada.
Christian Sievers
2
En pmul, use for j in range(n-i):y evite elif
Christian Sievers
1
@Lembik Calcula toda la matriz; por un factor de aproximadamente dos, reemplace j != kcon j < k. Copia una submatriz en el caso else que se puede evitar cuando manejamos y eliminamos las dos últimas en lugar de las dos primeras filas y columnas. Y cuando calcula con x={1,2,4}y más tarde con x={1,2,4,6}, repite sus cálculos hasta i=5. Reemplacé la estructura de los dos bucles externos con el primer bucle activado iy luego asumiendo recursivamente i in Xy i not in X. - Por cierto, puede ser interesante observar el crecimiento del tiempo necesario en comparación con los otros programas más lentos.
Christian Sievers
1

Octava

Esto es básicamente una copia de la entrada de Dennis , pero optimizado para Octave. La optimización principal se realiza utilizando la matriz de entrada completa (y su transposición) y la recursión utilizando solo índices de matriz, en lugar de crear matrices reducidas.

La principal ventaja es la copia reducida de matrices. Mientras que Octave no tiene una diferencia entre punteros / referencias y valores y funcionalmente solo pasa por valor, es una historia diferente detrás de escena. Allí, se utiliza la copia en escritura (copia diferida). Eso significa que, para el código a=1;b=a;b=b+1, la variable bsolo se copia en una nueva ubicación en la última instrucción, cuando se cambia. Dado matiny matranspno se cambian, ellos nunca por copiado. La desventaja es que la función pasa más tiempo calculando los índices correctos. Puede que tenga que probar diferentes variaciones entre índices numéricos y lógicos para optimizar esto.

Nota importante: la matriz de entrada debe ser int32! Guarde la función en un archivo llamadohaf.m

function h=haf(matin,indices,matransp,transp)

    if nargin-4
        indices=int32(1:length(matin));
        matransp=matin';
        transp=false;
    end
    if(transp)
        matrix=matransp;
    else
        matrix=matin;
    end
    ind1=indices(1);
    n=length(indices);
    if n==2
        h=matrix(ind1,indices(2));
        return
    end
    h=0*matrix(1); 
    for j=1:(n-1)
        indj=indices(j+1);
        k=matrix(ind1,indj);
        if logical(k)
            indicestemp=true(n,1);
            indicestemp(1:j:j+1)=false;
            h=h+k.*haf(matin,indices(indicestemp),matransp,~transp);
        end
    end
end

Ejemplo de script de prueba:

matrix = int32([0 0 1 -1 1 0 -1 -1 -1 0 -1 1 0 1 1 0 0 1 0 0 1 0 1 1;0 0 1 0 0 -1 -1 -1 -1 0 1 1 1 1 0 -1 -1 0 0 1 1 -1 0 0;-1 -1 0 1 0 1 -1 1 -1 1 0 0 1 -1 0 0 0 -1 0 -1 1 0 0 0;1 0 -1 0 1 1 0 1 1 0 0 0 1 0 0 0 1 -1 -1 -1 -1 1 0 -1;-1 0 0 -1 0 0 1 -1 0 1 -1 -1 -1 1 1 0 1 1 1 0 -1 1 -1 -1;0 1 -1 -1 0 0 1 -1 -1 -1 0 -1 1 0 0 0 -1 0 0 1 0 0 0 -1;1 1 1 0 -1 -1 0 -1 -1 0 1 1 -1 0 1 -1 0 0 1 -1 0 0 0 -1;1 1 -1 -1 1 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 -1 1 0 0;1 1 1 -1 0 1 1 0 0 -1 1 -1 1 1 1 0 -1 -1 -1 -1 0 1 1 -1;0 0 -1 0 -1 1 0 -1 1 0 1 0 0 0 0 0 1 -1 0 0 0 1 -1 -1;1 -1 0 0 1 0 -1 0 -1 -1 0 0 1 0 0 -1 0 -1 -1 -1 -1 -1 1 -1;-1 -1 0 0 1 1 -1 -1 1 0 0 0 -1 0 0 -1 0 -1 -1 0 1 -1 0 0;0 -1 -1 -1 1 -1 1 0 -1 0 -1 1 0 1 -1 -1 1 -1 1 0 1 -1 1 -1;-1 -1 1 0 -1 0 0 0 -1 0 0 0 -1 0 0 -1 1 -1 -1 0 1 0 -1 -1;-1 0 0 0 -1 0 -1 0 -1 0 0 0 1 0 0 1 1 1 1 -1 -1 0 -1 -1;0 1 0 0 0 0 1 0 0 0 1 1 1 1 -1 0 0 1 -1 -1 -1 0 -1 -1;0 1 0 -1 -1 1 0 -1 1 -1 0 0 -1 -1 -1 0 0 -1 1 0 0 -1 -1 1;-1 0 1 1 -1 0 0 0 1 1 1 1 1 1 -1 -1 1 0 1 1 -1 -1 -1 1;0 0 0 1 -1 0 -1 -1 1 0 1 1 -1 1 -1 1 -1 -1 0 1 1 0 0 -1;0 -1 1 1 0 -1 1 0 1 0 1 0 0 0 1 1 0 -1 -1 0 0 0 1 0;-1 -1 -1 1 1 0 0 1 0 0 1 -1 -1 -1 1 1 0 1 -1 0 0 0 0 0;0 1 0 -1 -1 0 0 -1 -1 -1 1 1 1 0 0 0 1 1 0 0 0 0 1 0;-1 0 0 0 1 0 0 0 -1 1 -1 0 -1 1 1 1 1 1 0 -1 0 -1 0 1;-1 0 0 1 1 1 1 0 1 1 1 0 1 1 1 1 -1 -1 1 0 0 0 -1 0])

tic
i=1;
while(toc<60)
    tic
    haf(matrix(1:i,1:i));
    i=i+1;
end

He probado esto usando TIO y MATLAB (en realidad nunca he instalado Octave). Me imagino que hacerlo funcionar es tan simple como sudo apt-get install octave. El comando octavecargará la interfaz gráfica de usuario de Octave. Si es más complicado que esto, eliminaré esta respuesta hasta que haya proporcionado instrucciones de instalación más detalladas.

Sanchises
fuente
0

Recientemente Andreas Bjorklund, Brajesh Gupt y yo publicamos un nuevo algoritmo para hafnianos de matrices complejas: https://arxiv.org/pdf/1805.12498.pdf . Para una matriz n \ times n se escala como n ^ 3 2 ^ {n / 2}.

Si entiendo correctamente el algoritmo original de Andreas de https://arxiv.org/pdf/1107.4466.pdf se escala como n ^ 4 2 ^ {n / 2} o n ^ 3 log (n) 2 ^ {n / 2} si usaste transformadas de Fourier para hacer multiplicaciones polinómicas.
Nuestro algoritmo está específicamente diseñado para matrices complejas, por lo que no será tan rápido como los desarrollados aquí para matrices {-1,0,1}. Sin embargo, me pregunto si uno puede usar algunos de los trucos que usó para mejorar nuestra implementación. Además, si las personas están interesadas, me gustaría ver cómo funciona su implementación cuando se les dan números complejos en lugar de números enteros. Finalmente, cualquier comentario, crítica, mejoras, errores, mejoras son bienvenidos en nuestro repositorio https://github.com/XanaduAI/hafnian/

¡Salud!

Nicolás Quesada
fuente
Bienvenido al sitio! Sin embargo, las respuestas a esta pregunta deben contener código. Sería mejor dejarlo como un comentario (que lamentablemente no tiene el representante para hacer).
Ad Hoc Garf Hunter
Bienvenido a PPCG. Si bien su respuesta puede hacer un buen comentario, el sitio no es para control de calidad. Este es un sitio de desafío y la respuesta a un desafío debe tener código y no una explicación de otra cosa. Actualice o elimine amablemente (si no lo hace, los mods lo harán)
Muhammad Salman
Bueno, el código está en github, pero supongo que no tiene sentido simplemente copiarlo y pegarlo aquí.
Nicolás Quesada
2
Si cabe en una respuesta, especialmente si usted es uno de los autores, no creo que haya nada malo en publicar una solución competitiva debidamente atribuida que se haya publicado en otro lugar.
Dennis
@ NicolásQuesada Las respuestas en este sitio deberían ser independientes si es posible, lo que significa que no deberíamos tener que ir a otro sitio para ver su respuesta / código.
mbomb007