Codegolf the Hafnian The challenge is to write codegolf for the Hafnian of a matrix.
The Hafnian of an 2n-by-2n symmetric matrix A is defined as:



Here S2n represents the set of all permutations of the integers from 1 to 2n, that is [1, 2n].

The wikipedia link talks about adjacency matrices but your code should work for any real valued symmetric input matrices.  

For those interested in applications of the Hafnian, the mathoverflow link discusses some more.

Your code can take input however it wishes and give output in any sensible format but please include in your answer a full worked example including clear instructions for how to supply input to your code. 

The input matrix is always square and will be at most 16 by 16. There is no need to be able to handle the empty matrix or matrices of odd dimension.

Reference implementation

Here is some example python code from Mr. 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




The wiki page has now (March 2 2018) been updated by  ShreevatsaR to include a different way of calculating the Hafnian.  It would be very interesting to see this golfed.
R, 150 142 127 119 bytes



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)


Try it online!

Uses the same trick I discovered golfing down this answer to index the matrix P, and @Vlo suggested an approach to entirely remove the for loop for -6 bytes!

To create a new test case, you can do matrix(c(values,separated,by,commas,going,across,rows),nrow=2n,ncol=2n,byrow=T).

Explanation: (the code is the same; it uses an apply rather than a for loop but the logic is otherwise identical).

function(A){
N &lt;- nrow(A)                   #N = 2*n
k &lt;- 1:(N/2) * 2               #k = c(2,4,...,N) -- aka 2*j in the formula
P &lt;- gtools::permutations(N,N) #all N-length permutations of 1:N
for(i in 1:nrow(P))
 F &lt;- 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
}