#!/usr/bin/env python
# coding: utf-8

# # Comparing the ground state energies obtained by density matrix renormalization group, exact diagonalization, and an SDP hierarchy

# We would like to compare the ground state energy of the following spinless fermionic system [1]:
# 
# $H_{\mathrm{free}}=\sum_{<rs>}\left[c_{r}^{\dagger} c_{s}+c_{s}^{\dagger} c_{r}-\gamma(c_{r}^{\dagger} c_{s}^{\dagger}+c_{s}c_{r} )\right]-2\lambda\sum_{r}c_{r}^{\dagger}c_{r},$
# 
# where $<rs>$ goes through nearest neighbour pairs in a two-dimensional lattice. The fermionic operators are subject to the following constraints:
# 
# $\{c_{r}, c_{s}^{\dagger}\}=\delta_{rs}I_{r}$
# 
# $\{c_r^\dagger, c_s^\dagger\}=0,$
# 
# $\{c_{r}, c_{s}\}=0.$
# 

# Our primary goal is to benchmark the SDP hierarchy of Reference [2]. The baseline methods are density matrix renormalization group (DMRG) and exact diagonalization (ED), both of which are included in Algorithms and Libraries for Physics Simulations (ALPS, [3]). The range of predefined Hamiltonians is limited, so we simplify the equation by setting $\gamma=0$.

# # Prerequisites

# To run this notebook, [ALPS](http://alps.comp-phys.org/mediawiki/index.php/Main_Page), [Sympy](http://sympy.org/), [Scipy](http://scipy.org/), and [SDPA](http://sdpa.sourceforge.net/) must be installed. A recent version of [Ncpol2sdpa](https://pypi.python.org/pypi/ncpol2sdpa) is also necessary.

# # Calculating the ground state energy with DMRG and ED

# DMRG and ED are included in ALPS. To start the calculations, we need to import the Python interface:

# In[1]:


import pyalps


# For now, we are only interested in relatively small systems, we will try lattice sizes between $2\times 2$ and $5\times 5$. With this, we set the parameters for DMRG and ED:

# In[2]:


lattice_range = [2, 3, 4, 5]
parms = [{ 
  'LATTICE'        : "open square lattice",  # Set up the lattice
  'MODEL'          : "spinless fermions",    # Select the model 
  'L'              : L,                      # Lattice dimension
  't'              : -1 ,                    # This and the following
  'mu'             : 2,                      # are parameters to the
  'U'              : 0 ,                     # Hamiltonian.
  'V'              : 0,
  'Nmax'           : 2 ,                     # These parameters are
  'SWEEPS'         : 20,                      # specific to the DMRG
  'MAXSTATES'      : 300,                    # solver.
  'NUMBER_EIGENVALUES' : 1,          
  'MEASURE_ENERGY' : 1
} for L in lattice_range ]


# We will need a helper function to extract the ground state energy from the solutions:

# In[3]:


def extract_ground_state_energies(data):
    E0 = []
    for Lsets in data:
        allE = []
        for q in pyalps.flatten(Lsets):
            allE.append(q.y[0])
        E0.append(allE[0])
    return sorted(E0, reverse=True)


# We invoke the solvers and extract the ground state energies from the solutions. First we use exact diagonalization, which, unfortunately does not scale beyond a lattice size of $4\times 4$. 

# In[4]:


prefix_sparse = 'comparison_sparse'
input_file_sparse = pyalps.writeInputFiles(prefix_sparse, parms[:-1])

res = pyalps.runApplication('sparsediag', input_file_sparse)
sparsediag_data = pyalps.loadEigenstateMeasurements(
                     pyalps.getResultFiles(prefix=prefix_sparse)) 

sparsediag_ground_state_energy = extract_ground_state_energies(sparsediag_data)
sparsediag_ground_state_energy.append(0)


# DMRG scales to all the lattice sizes we want:

# In[5]:


prefix_dmrg = 'comparison_dmrg'
input_file_dmrg = pyalps.writeInputFiles(prefix_dmrg, parms)
res = pyalps.runApplication('dmrg',input_file_dmrg)
dmrg_data = pyalps.loadEigenstateMeasurements(
                  pyalps.getResultFiles(prefix=prefix_dmrg)) 
dmrg_ground_state_energy = extract_ground_state_energies(dmrg_data)


# # Calculating the ground state energy with SDP

# The ground state energy problem can be rephrased as a polynomial optimiziation problem of noncommuting variables. We use Ncpol2sdpa to translate this optimization problem to a sparse SDP relaxation [4]. The relaxation is solved with SDPA, a high-performance SDP solver that deals with sparse problems efficiently [5]. First we need to import a few more functions:

# In[6]:


from sympy.physics.quantum.dagger import Dagger
from ncpol2sdpa import SdpRelaxation, generate_operators, \
                       fermionic_constraints, get_neighbors


# We set the additional parameters for this formulation, including the order of the relaxation:

# In[7]:


level = 1
gam, lam = 0, 1


# Then we iterate over the lattice range, defining a new Hamiltonian and new constraints in each step:

# In[8]:


sdp_ground_state_energy = []
for lattice_dimension in lattice_range:
    n_vars = lattice_dimension * lattice_dimension
    C = generate_operators('C%s' % (lattice_dimension), n_vars)
    
    hamiltonian = 0
    for r in range(n_vars):
        hamiltonian -= 2*lam*Dagger(C[r])*C[r]
        for s in get_neighbors(r, lattice_dimension):
            hamiltonian += Dagger(C[r])*C[s] + Dagger(C[s])*C[r]
            hamiltonian -= gam*(Dagger(C[r])*Dagger(C[s]) + C[s]*C[r])
    
    substitutions = fermionic_constraints(C)
        
    sdpRelaxation = SdpRelaxation(C)
    sdpRelaxation.get_relaxation(level, objective=hamiltonian, substitutions=substitutions)
    sdpRelaxation.solve()
    sdp_ground_state_energy.append(sdpRelaxation.primal)


# # Comparison

# The level-one relaxation matches the ground state energy given by DMRG and ED.

# In[9]:


data = [dmrg_ground_state_energy,\
        sparsediag_ground_state_energy,\
        sdp_ground_state_energy]
labels = ["DMRG", "ED", "SDP"]
print ("{:>4} {:>9} {:>10} {:>10} {:>10}").format("", *lattice_range)
for label, row in zip(labels, data):
    print ("{:>4} {:>7.6f} {:>7.6f} {:>7.6f} {:>7.6f}").format(label, *row)


# # References

# [1] Corboz, P.; Evenbly, G.; Verstraete, F. & Vidal, G. [Simulation of interacting fermions with entanglement renormalization](http://arxiv.org/abs/0904.4151). _Physics Review A_, 2010, 81, pp. 010303.
# 
# [2] Pironio, S.; Navascués, M. & Acín, A. [Convergent relaxations of polynomial optimization problems with noncommuting variables](http://arxiv.org/abs/0903.4368). _SIAM Journal on Optimization_, 2010, 20, pp. 2157-2180.
# 
# [3] Bauer, B.; Carr, L.; Evertz, H.; Feiguin, A.; Freire, J.; Fuchs, S.; Gamper, L.; Gukelberger, J.; Gull, E.; Guertler, S. & others. [The ALPS project release 2.0: Open source software for strongly correlated systems](http://arxiv.org/abs/1101.2646). _Journal of Statistical Mechanics: Theory and Experiment_, IOP Publishing, 2011, 2011, P05001.
# 
# [4] Wittek, P. [Ncpol2sdpa -- Sparse Semidefinite Programming Relaxations for Polynomial Optimization Problems of Noncommuting Variables](http://arxiv.org/abs/1308.6029). _arXiv:1308.6029_, 2013.
# 
# [5] Yamashita, M.; Fujisawa, K. & Kojima, M. [Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0)](http://dx.doi.org/10.1080/1055678031000118482). _Optimization Methods and Software_, 2003, 18, 491-505.