TFIM Performance Comparison

quickbench.jl

using DataFrames
using Gadfly
using PyCall
@pyimport timeit
include("tfim_1D_OBC.jl")

# Max number of sites, large values might run out of memory
maxN = 15

# Time python code
pytimes = collect(1.0:maxN)
for n in 1:maxN
  reps = Int(ceil(1E7/n^5))
  pyexpr = "tfim_1D_OBC.gen_TFIM_ham($(n),-1,-1)"
  pytimes[n]=timeit.timeit(pyexpr,setup="import tfim_1D_OBC",number=reps)/reps
end

# Function to time julia code
function repfunctime(n,reps)
  tic()
  for i in 1:reps
    gen_TFIM_ham(n,-1.0,-1.0)
  end
  toc()
end

# Time julia code
jltimes = collect(1.0:maxN)
for n in 1:maxN
  reps = Int(ceil(1E7/n^5))
  jltimes[n]=repfunctime(n,reps)/reps
end

# Gather everything in a dataframe
df1 = DataFrame(N=collect(1:maxN), t=jltimes, Language="Julia")
df2 = DataFrame(N=collect(1:maxN), t=pytimes, Language="Python")
df = vcat(df1, df2)

# Plot with nice colours
p = plot(df,x="N", y="t", color="Language",Geom.line,Geom.point,Scale.y_log10)
draw(PDF("tfimbench.pdf", 12cm, 9cm), p)

tfim_1D_OBC.jl

function gen_TFIM_ham(N,J,H)
# A simple diagonalization for the 1D TFIM with OBCs
# N is the number of sites in chain
# J and H are the Hamiltonian parameters
# Hilbert Space Dimension
dim = 2^N
# Hamiltonian Matrix
ham = zeros(dim,dim)
#Build the diagonal part of the Hamiltonian
#First loop through the hilbert space
for bra in 0:dim-1
  # Next, loop through all the Hamiltonian elements
  diag_term = 0.0
  for s in 0:N-2
    S0 = 2*((bra>>s)&1) - 1
    S1 = 2*((bra>>(s+1))&1) - 1
    diag_term = diag_term + J*S0*S1
  end
  ham[bra+1,bra+1] = diag_term
  #Now let's do the off-diagonal terms
  for s in 0:N-1
    y = 2^s  #This labels the bit to be flipped
    ket = bra $ y     #Flip that bit in the bra
    ham[bra+1,ket+1] = H
  end
end
return ham
end

tfim_1D_OBC.py

import numpy as np
# A simple diagonalization for the 1D TFIM with OBCs
# N is the number of sites in chain
# J and H are the Hamiltonian parameters
def gen_TFIM_ham(N,J,H):
    #Hilbert Space Dimension
    Dim = 2**N
    #Hamiltonian Matrix
    Ham = np.zeros((Dim,Dim)) ##creates a Dim by Dim matrix set to zero
    #Build the diagonal part of the Hamiltonian
    #First loop through the hilbert space
    for bra in range (0,Dim):
        #Next, loop through all the Hamiltonian elements
        DiagTerm = 0
        for s in range (0,N-1):    #OBC only
            S0 = 2*((bra>>s)&1) - 1
            S1 = 2*((bra>>(s+1))&1) - 1
            DiagTerm += J*S0*S1
        Ham[bra,bra] = DiagTerm   #diagonal matrix element
        #Now let's do the off-diagonal terms
        for s in range (0,N):
                y = 2**s  #This labels the bit to be flipped
                ket = bra^y     #Flip that bit in the bra
                Ham[bra,ket] = H
    return Ham

These files were used for a simple benchmark, see the relevant discussion here:
https://groups.google.com/forum/#!topic/julia-users/aW4rjUIFq6w

I updated the code to have the Ham[bra,bra] = DiagTerm line outside the loop, as in the Julia code.