Codegolf le Hafnian

22

Le défi est d'écrire du codegolf pour le Hafnian d'une matrice . Le Hafnien d'une matrice 2n-par- 2nsymétrique Aest défini comme:

entrez la description de l'image ici

Ici, S 2n représente l'ensemble de toutes les permutations des entiers de 1à 2n, c'est-à-dire [1, 2n].

Le lien wikipedia parle des matrices d'adjacence mais votre code devrait fonctionner pour toutes les matrices d'entrée symétriques à valeur réelle.

Pour ceux qui s'intéressent aux applications du Hafnian, le lien mathoverflow en discute plus.

Votre code peut prendre une entrée comme il le souhaite et donner une sortie dans n'importe quel format raisonnable, mais veuillez inclure dans votre réponse un exemple complet comprenant des instructions claires sur la façon de fournir une entrée à votre code.

La matrice d'entrée est toujours carrée et sera au maximum de 16 sur 16. Il n'est pas nécessaire de pouvoir gérer la ou les matrices vides de dimension impaire.

Implémentation de référence

Voici un exemple de code python de M. 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([[0, 4.5], [4.5, 0]]))
4.5
print(hafnian([[0, 4.7, 4.6, 4.5], [4.7, 0, 2.1, 0.4], [4.6, 2.1, 0, 1.2], [4.5, 0.4, 1.2, 0]])
16.93
print(hafnian([[1.3, 4.1, 1.2, 0.0, 0.9, 4.4], [4.1, 4.2, 2.7, 1.2, 0.4, 1.7], [1.2, 2.7, 4.9, 4.7, 4.0, 3.7], [0.0, 1.2, 4.7, 2.2, 3.3, 1.8], [0.9, 0.4, 4.0, 3.3, 0.5, 4.4], [4.4, 1.7, 3.7, 1.8, 4.4, 3.2]])
262.458

La page wiki a maintenant (2 mars 2018) été mise à jour par ShreevatsaR pour inclure une manière différente de calculer le Hafnian. Il serait très intéressant de voir ce golf.

user202729
la source
5
Je pense que ce serait plus facile à digérer avec une explication informelle du hafnien. Quelque chose comme, prenez tous les sous-ensembles de n entrées de matrice où leurs n indices de ligne et n indices de colonne forment une partition de 1..2n, prenez le produit de chacun, ajoutez-les et mettez la somme à l'échelle.
xnor

Réponses:

9

R , 150 142 127 127 119 octets

function(A,N=nrow(A),k=1:(N/2)*2)sum(apply(gtools::permutations(N,N),1,function(r)prod(A[cbind(r[k-1],r[k])])))/prod(k)

Essayez-le en ligne!

Utilise la même astuce que j'ai découverte en descendant cette réponse pour indexer la matrice P, et @Vlo a suggéré une approche pour supprimer entièrement la forboucle pour -6 octets!

Pour créer un nouveau scénario de test, vous pouvez le faire matrix(c(values,separated,by,commas,going,across,rows),nrow=2n,ncol=2n,byrow=T).

Explication: (le code est le même; il utilise une boucle applyplutôt qu'une forboucle mais la logique est par ailleurs identique).

function(A){
N <- nrow(A)                   #N = 2*n
k <- 1:(N/2) * 2               #k = c(2,4,...,N) -- aka 2*j in the formula
P <- gtools::permutations(N,N) #all N-length permutations of 1:N
for(i in 1:nrow(P))
 F <- F + prod(A[cbind(P[i,k-1],P[i,k])]) # takes the product of all A_sigma(2j-1)sigma(2j) for fixed i and adds it to F (initialized to 0)
F / prod(k)                    #return value; prod(k) == n! * 2^n
}

Giuseppe
la source
Appliquer est moins cher de 2 octets, ce qui permet d'économiser 4 octets supplémentaires en entassant les autres lignes ensemble. tio.run/##PY6xDoIwEIZ3nsLxzpxiS4ymkYEXYHIjDFDEEKBtSokS47PX4sDw5/… Il est également très intéressant de voir comment la base R n'a pas de fonction de permutation pour un langage de programmation statistique
Vlo
@Vlo très sympa! nous pouvons déplacer Net kdans les arguments de fonction pour obtenir une instruction, en supprimant {}et en enregistrant deux autres octets.
Giuseppe
@Giuseppe Darn oublie toujours que vous pouvez les définir dans la fonction args. J'ai passé quelques minutes à essayer de braquer ces variables ...
Vlo
8

Pyth , 24 octets

sm*Fmc@@Qhkek2d{mScd2.pU

Essayez-le ici!


Ancienne version, 35 octets

*c1**FK/lQ2^2Ksm*Fm@@Q@[email protected]

Essayez-le ici!

M. Xcoder
la source
3
Actuellement en tête mais il faut craindre les réponses de Jelly à venir .... :)
Eh Jelly va certainement me battre d'environ 10 octets. Pyth n'est pas le meilleur outil pour le travail
M. Xcoder
05AB1E semble même pouvoir lier Pyth (croyez-le ou non, enfin un défi matriciel où il a[b]suffit de rivaliser).
Urne de poulpe magique
@MagicOctopusUrn J'ai déjà une solution 05AB1E qui bat Pyth :-) Je ne vais pas la publier (pour l'instant, au moins)
M. Xcoder
Est-ce quelque chose du xÍysè<¹sès·<ysè<ègenre lmao? PS Mine fait 40 octets et ne fonctionne pas si bien, alors n'hésitez pas à le poster, je ne suis pas sûr de pouvoir terminer avant de rentrer chez moi.
Urne de poulpe magique
6

Stax , 23 22 19 17 octets

ü;Y╙◘▌Φq↓ê²╧▐å↑┌C

Exécutez-le et déboguez-le en ligne

La représentation ascii correspondante du même programme est la suivante.

%r|TF2/{xsE@i^H/m:*+

Le programme souffre d'une erreur d'arrondi en virgule flottante. En particulier, il rapporte 33673.5000000011au lieu de 33673.5. Mais je pense que la précision est acceptable étant donné que ce programme fonctionne sur des valeurs à virgule flottante. C'est aussi très lent, prenant presque une minute pour les exemples d'entrées sur cette machine.

%                             get size of matrix
 r|T                          get all permutations of [0 ... size-1]
    F                         for each, execute the rest of the program
     2/                       get consecutive pairs
       {        m             map each pair... 
        xsE@                      the matrix element at that location
            i^H/                  divided by 2*(i+1) where i=iteration index
                 :*           product of array
                   +          add to running total
récursif
la source
1
Très impressionnant!
5

05AB1E , 21 octets

ā<œε2ô{}Ùεε`Isèsè;]PO

Essayez-le en ligne!


Ancienne version, 32 octets

āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/

Essayez-le en ligne!

Comment ça marche?

āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/ – Full program. Argument: A matrix M.
ā                                – The range [1 ... len(M)].
 œ                               – Permutations.
  v                    }         – Iterate over the above with a variable y.
   Ig;©                          – Push len(M) / 2 and also store it in register c.
       Lε            }           – For each integer in the range [1 ... ^]:
         ·U                      – Double it and store it in a variable X.
            yX<                  – Push the element of y at index X-1.
           I   è                 – And index with the result into M.
                yXè              – Push the element of y at index X.
                   è             – And index with the result into ^^.
                      P          – Take the product of the resulting list.
                        O        – Sum the result of the mapping.
                         θ       – And take the last element*.
                          ®!     – Take the factorial of the last item in register c.
                             ®o  – Raise 2 to the power of the last item in register c.
                            /  / – And divide the sum of the mapping accordingly.

* – Yeah, this is needed because I mess up the stack when pushing so many values in the loop and not popping correctly ;P
M. Xcoder
la source
1
Sans blague èsè, hah ... haha ​​... je suis ringard.
Magic Octopus Urn
@MagicOctopusUrn Fixed ... J'ai oublié que 05AB1E est indexé 0> _ <
M. Xcoder
3

Gelée , 19 octets

LŒ!s€2Ṣ€QḅL_LịFHP€S

Essayez-le en ligne!

Version alternative, 15 octets, défi postdates

LŒ!s€2Ṣ€QœịHP€S

Jelly a finalement obtenu l'indexation des tableaux à n dimensions.

Essayez-le en ligne!

Comment ça marche

LŒ!s€2Ṣ€QœiHP€S  Main link. Argument: M (matrix / 2D array)

L                Take the length, yielding 2n.
 Œ!              Generate all permutations of [1, ..., 2n].
   s€2           Split each permutation into pairs.
      Ṣ€         Sort the pair arrays.
        Q        Unique; deduplicate the array of pair arrays.
                 This avoids dividing by n! at the end.
           H     Halve; yield M, with all of its elements divided by 2.
                 This avoids dividing by 2**n at the end.
         œị      At-index (n-dimensional); take each pair of indices [i, j] and
                 yield M[i][j].
            P€   Take the product the results corresponding the same permutation.
              S  Take the sum of the products.

La version 19 octets fonctionne de manière similaire; il suffit de se mettre en œuvre œị.

...ḅL_LịFH...    Return value: Array of arrays of index pairs. Argument: M

    L            Length; yield 2n.
   ḅ             Convert each pair of indices [i, j] from base 2n to integer,
                 yielding ((2n)i + j).
     _L          Subtract 2n, yielding ((2n)(i - 1) + j).
                 This is necessary because indexing is 1-based in Jelly, so the
                 index pair [1, 1] must map to index 1.
        F        Yield M, flattened.
       ị         Take the indices to the left and get the element at these indices
                 from the array to the right.
         H       Halve; divide all retrieved elements by 2.
Dennis
la source
3

C (gcc) , 288 285 282 293 292 272 271 octets

  • Enregistré trois octets en jouant avec deux post-incréments et pour le placement de la boucle.
  • Enregistrement de trois octets en manipulant une autre post-incrémentation, en déplaçant les deux initialisations de variables avant la branche - joué if(...)...k=0...else...,j=0...au golf if(k=j=0,...)...else...- et en effectuant un décalage d'index.
  • Onze octets requis en prenant en charge les floatmatrices.
  • Enregistré un octet grâce à M. Xcoder ; jouer 2*j+++1au golf j-~j++.
  • Enregistrement de vingt octets en supprimant une intdéclaration de type de variable superflue et en n'utilisant pas de fonction factorielle, mais en calculant la valeur factorielle à l'aide d'une boucle for déjà existante.
  • Enregistré un octet en jouant S=S/F/(1<<n);au golf S/=F*(1<<n);.
float S,p,F;j,i;s(A,n,P,l,o,k)float*A;int*P;{if(k=j=0,o-l)for(;k<l;s(A,n,P,l,o+1))P[o]=k++;else{for(p=-l;j<l;j++)for(i=0;i<l;)p+=P[j]==P[i++];if(!p){for(F=p=1,j=0;j<n;F*=j)p*=A[P[2*j]*2*n+P[j-~j++]];S+=p;}}}float h(A,n)float*A;{int P[j=2*n];S=0;s(A,n,P,j,0);S/=F*(1<<n);}

Essayez-le en ligne!

Explication

float S,p,F;                    // global float variables: total sum, temporary, factorial
j,i;                            // global integer variables: indices
s(A,n,P,l,o,k)float*A;int*P;{   // recursively look at every permutation in S_n
 if(k=j=0,o-l)                  // initialize k and j, check if o != l (possible  permutation not yet fully generated)
  for(;k<l;s(A,n,P,l,o+1))      // loop through possible values for current possible  permuation position
   P[o]=k++;                    // set possible  permutation, recursively call (golfed into the for loop)
 else{for(p=-l;j<l;j++)         // there exists a possible permutation fully generated
  for(i=0;i<l;)                 // test if the possible permutation is a bijection
   p+=P[j]==P[i++];             // check for unique elements
  if(!p){                       // indeed, it is a permutation
   for(F=p=1,j=0;j<n;F*=j)      // Hafnian product loop and calculate the factorial (over and over to save bytes)
    p*=A[P[2*j]*2*n+P[j-~j++]]; // Hafnian product
   S+=p;}}}                     // add to sum
float h(A,n)float*A;{           // Hafnian function
 int P[j=2*n];S=0;              // allocate permutation memory, initialize sum
 s(A,n,P,j,0);                  // calculate Hafnian sum
 S/=F*(1<<n);}                  // calculate Hafnian

Essayez-le en ligne!

Au cœur du programme se trouve le générateur de permutation suivant qui boucle S_n. Tous les calculs hafniens sont simplement construits dessus - et plus loin encore.

j,i,p;Sn(A,l,o,k)int*A;{          // compute every element in S_n
 if(o-l)                          // o!=l, the permutation has not fully been generated
  for(k=0;k<l;k++)                // loop through the integers [0, n)
   A[o]=k,Sn(A,l,o+1);            // modify permutation, call recursively
 else{                            // possible permutation has been generated
  for(p=-l,j=0;j<l;j++)           // look at the entire possible permutation
   for(i=0;i<l;i++)p+=A[j]==A[i]; // check that all elements appear uniquely
  if(!p)                          // no duplicat elements, it is indeed a permutation
   for(printf("["),j=0;j<l        // print
   ||printf("]\n")*0;)            //  the
    printf("%d, ",A[j++]);}}      //   permutation
main(){int l=4,A[l];Sn(A,l,0);}   // all permutations in S_4

Essayez-le en ligne!

Jonathan Frech
la source
1
Super d'avoir une réponse C mais, comme vous le suggérez, elle n'est pas conforme actuellement.
@Lembik Fixed. Prend désormais en charge les floatmatrices.
Jonathan Frech
2*j+++1est équivalent à j+j+++1, ce qui est identique à j-(-j++-1), nous pouvons donc utiliser efficacement le complément au niveau du bit pour enregistrer un octet: j-~j++( Essayez-le en ligne )
M. Xcoder
3

R , 84 78 octets

h=function(m)"if"(n<-nrow(m),{for(j in 2:n)F=F+m[1,j]*h(m[v<--c(1,j),v]);F},1)

Essayez-le en ligne!

Edit: Merci à Vlo pour -6 octets.

Il semble que tout le monde ici implémente l'algorithme de référence standard avec permutations, mais j'ai essayé de tirer parti des connaissances de la communauté acquises dans le défi associé , qui est essentiellement la même tâche ciblée pour le code le plus rapide au lieu du golf.

Il s'avère que pour un langage qui est bon pour trancher des matrices (comme R), l'algorithme récursif: hafnian(m) = sum(m[i,j] * hafnian(m[-rows and columns at i,j])est non seulement plus rapide, mais aussi assez golfique. Voici le code non golfé:

hafnian<-function(m)
{
    n=nrow(m)
    #Exits one step earlier than golfed version
    if(n == 2) return(m[1,2])
    h = 0
    for(j in 2:n) {
        if(m[1,j] == 0) next
        h = h + m[1,j] * hafnian(m[c(-1,-j),c(-1,-j)])
    }
    h
}
Kirill L.
la source
Très belle réponse. -1 pour appeler Ifentre parenthèses, -4 pour utiliser Fcomme variable initialisée, -1 pour attribuer ndans le if. tio.run/##XU/LCsIwELz7FcFTVtOQl1pf1/…
Vlo
soigné! Je dirais de le poster dans le défi de la vitesse, mais il y a probablement d'autres optimisations (comme le filetage) qui peuvent être faites, et bien que R ne soit pas exactement connu pour sa vitesse, il serait bon de l'avoir là comme référence .
Giuseppe
Faites-le à des fins de référence!
Vlo
J'ai en fait essayé de tester cela pour la vitesse, mais j'ai été rapidement découragé par les résultats. La soumission Python la plus lente dans le défi de la vitesse en utilisant le même algorithme exact craque la matrice 24x24 en quelques secondes sur TIO, mais R expire. Sur ma machine locale, il n'a pas non plus répondu dans un délai raisonnable, même lorsqu'il a été aidé avec la mémorisation du paquet 'mémo' ...
Kirill L.
2

Gelée , 29 octets

LHµ2*×!
LŒ!s€2;@€€Wị@/€€P€S÷Ç

Essayez-le en ligne!

Je pense que la ;@€€Wị@/€€P€partie peut probablement être tournée vers le bas. Besoin de revenir plus tard pour vérifier et ajouter une explication.

dylnan
la source
Identique à ma solution (sauf la J) avant de jouer au golf . Les esprits gelés pensent de la même façon ? source
user202729
J'ai pu le réduire un peu plus en refactorisant la partie que vous avez mentionnée ainsi que la division par 2 et la factorielle. LḶŒ!s€2ḅL‘ịFZµPS÷JḤ$P$ TIO
miles
@ user202729 haha ​​nice
dylnan
@miles Wow, c'est beaucoup d'économies. Je vais l'éditer dans ma réponse mais c'est assez différent alors n'hésitez pas à soumettre votre propre réponse si vous le souhaitez
dylnan
2

Haskell , 136 octets

-14 octets grâce aux ovs.

import Data.List
f m|n<-length m`div`2=sum[product[m!!(s!!(2*j-2)-1)!!(s!!(2*j-1)-1)/realToFrac(2*j)|j<-[1..n]]|s<-permutations[1..n*2]]

Essayez-le en ligne!

Pouah...

totalement humain
la source
2

MATL , 29 24 22 octets

Zy:Y@!"G@2eZ{)tn:E/pvs

Essayez-le en ligne! Ou vérifiez tous les cas de test: 1 , 2 , 3 .

Comment ça marche

Zy       % Size of (implicit) input: pushes [2*n 2*n], where the
         % input is a 2*n × 2*n matrix. 
:        % Range: gives row vector [1 2 ... 2*n]
Y@       % All permutation of that vector as rows of a matrix
!"       % For each permutation 
  G      %   Push input matrix
  @      %   Push current permutation
  2e     %   Reshape as a 2-row array
  Z{     %   Split rows into a cell array of size 2
  )      %   Reference indexing. With a cell array as index this
         %   applies element-wise indexing (similar to sub2ind).
         %   Gives a row vector with the n matrix entries selected
         %   by the current permutation
  t      %   Duplicate
  n:     %   Number of elements, range: this gives [1 2 ... n]
  E      %   Double, element-wise: gives [2 4 ... 2*n]
  /      %   Divide, element-wise
  p      %   Product
  vs     %   Vertically concatenate and sum
         % End (implicit). Display (implicit)
Luis Mendo
la source
1

Noix de coco , 165 145 128 127 octets

-1 octet merci à M. Xcoder

m->sum(reduce((o,j)->o*m[p[2*j]][p[j-~j]]/-~j/2,range(len(m)//2),1)for p in permutations(range(len(m))))
from itertools import*

Essayez-le en ligne!

ovs
la source
1

Perl 6, 86 octets

{my \n=$^m/2;^$m .permutations.map({[*] .map(->\a,\b{$m[a][b]})}).sum/(2**n*[*] 1..n)}
bb94
la source