def restore_CF_dbg(self):
        for n in xrange(self.N + 2):
            print (n, self.l[n].trace().real, self.r[n].trace().real,
                   mm.adot(self.l[n], self.r[n]).real)

        norm_r = mm.adot(self.uni_r.l[-1], self.r[self.N])
        norm_l = mm.adot(self.l[0], self.uni_l.r[-1])
        print "norm of uni_r: ", norm_r
        print "norm of uni_l: ", norm_l

        #r_m1 = self.eps_r(0, self.r[0])
        #print mm.adot(self.l[0], r_m1).real

        norm = mm.adot(self.l[0], self.r[0]).real
        
        try:
            h = sp.empty((self.N + 1), dtype=self.typ)
            for n in xrange(self.N + 1):
                h[n] = self.expect_2s(self.h_nn[n], n)
            h *= 1/norm
    
            #self.uni_l.A = self.A[0] #FIXME: Not sure how to handle this yet...
            self.uni_l.l[-1] = self.l[0]
            self.uni_l.calc_AA()
            h_left = self.uni_l.expect_2s(self.uni_l.ham.copy()) / norm_l
    
            #self.uni_r.A = self.A[self.N + 1]
            self.uni_r.r[-1] = self.r[self.N]
            self.uni_r.calc_AA()
            h_right = self.uni_r.expect_2s(self.uni_r.ham.copy()) / norm_r
    
            return h, h_left, h_right
        except AttributeError:
            return sp.array([0]), 0, 0

    def restore_CF_dbg(self):
        for n in xrange(self.N + 2):
            print (n, self.l[n].trace().real, self.r[n].trace().real,
                   mm.adot(self.l[n], self.r[n]).real)

        norm_r = mm.adot(self.uni_r.l, self.r[self.N])
        norm_l = mm.adot(self.l[0], self.uni_l.r)
        print "norm of uni_r: ", norm_r
        print "norm of uni_l: ", norm_l

        #r_m1 = self.eps_r(0, self.r[0])
        #print mm.adot(self.l[0], r_m1).real

        norm = mm.adot(self.l[0], self.r[0]).real

        h = sp.empty((self.N + 1), dtype=self.typ)
        for n in xrange(self.N + 1):
            h[n] = self.expect_2s(self.h_nn[n], n)
        h *= 1/norm

        self.uni_l.A = self.A[0]
        self.uni_l.l = self.l[0]
        self.uni_l.calc_AA()
        h_left = self.uni_l.expect_2s(self.uni_l.ham.copy()) / norm_l

        self.uni_r.A = self.A[self.N + 1]
        self.uni_r.r = self.r[self.N]
        self.uni_r.calc_AA()
        h_right = self.uni_l.expect_2s(self.uni_r.ham.copy()) / norm_r

        return h, h_left, h_right

 def calc_B1(self):
     """Calculate the optimal B1 given right gauge-fixing on B2..N and
     no gauge-fixing on B1.
     
     We use the non-norm-preserving K's, since the norm-preservation
     is not needed elsewhere. It is cleaner to subtract the relevant
     norm-changing terms from the K's here than to generate all K's
     with norm-preservation.
     """
     B1 = sp.empty_like(self.A[1])
     
     try:
         r1_i = self.r[1].inv()
     except AttributeError:
         r1_i = mm.invmh(self.r[1])
         
     try:
         l0_i = self.l[0].inv()
     except AttributeError:
         l0_i = mm.invmh(self.l[0])
     
     A0 = self.A[0]
     A1 = self.A[1]
     A2 = self.A[2]
     r1 = self.r[1]
     r2 = self.r[2]
     l0 = self.l[0]
     
     KLh = mm.H(self.u_gnd_l.K_left - l0 * mm.adot(self.u_gnd_l.K_left, self.r[0]))
     K2 = self.K[2] - r1 * mm.adot(self.l[1], self.K[2])
     
     C1 = self.C[1] - self.h_expect[1] * self.AA1
     C0 = self.C[0] - self.h_expect[0] * self.AA0
     
     for s in xrange(self.q[1]):
         try:
             B1[s] = A1[s].dot(r1_i.dot_left(K2))
         except AttributeError:
             B1[s] = A1[s].dot(K2.dot(r1_i))
         
         for t in xrange(self.q[2]):
             try:
                 B1[s] += C1[s, t].dot(r2.dot(r1_i.dot_left(mm.H(A2[t]))))
             except AttributeError:
                 B1[s] += C1[s, t].dot(r2.dot(mm.H(A2[t]).dot(r1_i)))                    
             
         B1sbit = KLh.dot(A1[s])
                         
         for t in xrange(self.q[0]):
             B1sbit += mm.H(A0[t]).dot(l0.dot(C0[t,s]))
             
         B1[s] += l0_i.dot(B1sbit)
        
     rb = sp.zeros_like(self.r[0])
     for s in xrange(self.q[1]):
         rb += B1[s].dot(r1.dot(mm.H(B1[s])))
     eta = sp.sqrt(mm.adot(l0, rb))
             
     return B1, eta

    def calc_K(self):
        """Generates the right K matrices used to calculate the B's

        K[n] is recursively defined. It depends on C[m] and A[m] for all m >= n.

        It directly depends on A[n], A[n + 1], r[n], r[n + 1], C[n] and K[n + 1].

        This is equivalent to K on p. 14 of arXiv:1103.0936v2 [cond-mat.str-el], except
        that it is for the non-norm-preserving case.

        K[1] is, assuming a normalized state, the expectation value H of Ĥ.
        
        Return the excess energy.
        """
        n_low = 0
        n_high = self.N + 1
            
        H = mm.H
        
        self.h_expect = sp.zeros((self.N + 1), dtype=self.typ)
        
        self.u_gnd_r.calc_AA()
        self.u_gnd_r.calc_C()
        self.u_gnd_r.calc_K()
        self.K[self.N + 1][:] = self.u_gnd_r.K

        for n in reversed(xrange(n_low, n_high)):
            self.K[n].fill(0)
            
            K = self.K[n]
            Kp1 = self.K[n + 1]
            C = self.C[n]
            rp1 = self.r[n + 1]
            A = self.A[n]
            Ap1 = self.A[n + 1]
            
            Hr = sp.zeros_like(K)

            for s in xrange(self.q[n]):
                Ash = H(A[s])
                for t in xrange(self.q[n+1]):
                    Hr += C[s, t].dot(rp1.dot(H(Ap1[t]).dot(Ash)))

                K += A[s].dot(Kp1.dot(Ash))
                
            self.h_expect[n] = mm.adot(self.get_l(n), Hr)
                
            K += Hr
            
        self.u_gnd_l.calc_AA()
        self.u_gnd_l.calc_C()
        K_left, h_left_uni = self.u_gnd_l.calc_K_l()

        h = (mm.adot(K_left, self.r[0]) + mm.adot(self.l[0], self.K[0]) 
             - (self.N + 1) * self.u_gnd_r.h)
             
        return h

    def expect_3s(self, op, n, AAA=None):
        """Computes the expectation value of a nearest-neighbour three-site operator.

        The operator should be a q[n] x q[n + 1] x q[n + 2] x q[n] x
        q[n + 1] x q[n + 2] array such that op[s, t, u, v, w, x] =
        <stu|op|vwx> or a function of the form op(s, t, u, v, w, x) =
        <stu|op|vwx>.

        The state must be up-to-date -- see self.update()!

        Parameters
        ----------
        op : ndarray or callable
            The operator array or function.
        n : int
            The leftmost site number (operator acts on n, n + 1, n + 2).

        Returns
        -------
        expval : floating point number
            The expectation value (data type may be complex)
        """
        if AAA is None:
            if not self.AAA[n] is None:
                AAA = self.AAA[n]
            else:
                AAA = tm.calc_AAA_AA(self.AA[n], self.A[n + 2])
                
        if op is self.ham[n] and self.ham_sites == 3:
            res = tm.eps_r_op_3s_C123_AAA456(self.r[n + 2], self.C[n], AAA)
            return m.adot(self.l[n - 1], res)
        else:
            return super(EvoMPS_MPS_Generic).expect_3s(op, n, AAA=AAA)

 def expect_1s(self, op):
     """Computes the expectation value of a single-site operator.
     
     The operator should be a self.q x self.q matrix or generating function 
     such that op[s, t] or op(s, t) equals <s|op|t>.
     
     The state must be up-to-date -- see self.update()!
     
     Parameters
     ----------
     op : ndarray or callable
         The operator.
         
     Returns
     -------
     expval : floating point number
         The expectation value (data type may be complex)
     """        
     if callable(op):
         op = np.vectorize(op, otypes=[np.complex128])
         op = np.fromfunction(op, (self.q, self.q))
         
     Or = tm.eps_r_op_1s(self.r, self.A, self.A, op)
     
     return m.adot(self.l, Or)

 def density_2s(self, d):
     """Returns a reduced density matrix for a pair of (seperated) sites.
     
     The site number basis is used: rho[s * q + u, t * q + v]
     with 0 <= s, t < q and 0 <= u, v < q.
     
     The state must be up-to-date -- see self.update()!
     
     Parameters
     ----------
     d : int
         The distance between the first and the second sites considered (d = n2 - n1).
         
     Returns
     -------
     rho : ndarray
         Reduced density matrix in the number basis.
     """
     rho = sp.empty((self.q * self.q, self.q * self.q), dtype=sp.complex128)
     
     for s2 in xrange(self.q):
         for t2 in xrange(self.q):
             r_n2 = m.mmul(self.A[t2], self.r, m.H(self.A[s2]))
             
             r_n = r_n2
             for n in xrange(d - 1):
                 r_n = tm.eps_r_noop(r_n, self.A, self.A)
                 
             for s1 in xrange(self.q):
                 for t1 in xrange(self.q):
                     r_n1 = m.mmul(self.A[t1], r_n, m.H(self.A[s1]))
                     tmp = m.adot(self.l, r_n1)
                     rho[s1 * self.q + s2, t1 * self.q + t2] = tmp
     return rho        

 def expect_2s(self, op):
     """Computes the expectation value of a nearest-neighbour two-site operator.
     
     The operator should be a q x q x q x q array 
     such that op[s, t, u, v] = <st|op|uv> or a function of the form 
     op(s, t, u, v) = <st|op|uv>.
     
     The state must be up-to-date -- see self.update()!
     
     Parameters
     ----------
     op : ndarray or callable
         The operator array or function.
         
     Returns
     -------
     expval : floating point number
         The expectation value (data type may be complex)
     """
     if callable(op):
         op = np.vectorize(op, otypes=[np.complex128])
         op = np.fromfunction(op, (self.q, self.q, self.q, self.q))        
     
     C = tm.calc_C_mat_op_AA(op, self.AA)
     res = tm.eps_r_op_2s_C12_AA34(self.r, C, self.AA)
     
     return m.adot(self.l, res)

    def calc_K(self):
        """Generates the right K matrices used to calculate the B's
        
        K[n] contains 'column-vectors' such that <l[n]|K[n]> = trace(l[n].dot(K[n])).
        K_l[n] contains 'bra-vectors' such that <K_l[n]|r[n]> = trace(K_l[n].dot(r[n])).
        """
   
        self.h_expect = sp.zeros((self.N + 1), dtype=self.typ)
        
        self.uni_r.calc_AA()
        self.uni_r.calc_C()
        self.uni_r.calc_K()
        self.K[self.N + 1][:] = self.uni_r.K[0]

        self.uni_l.calc_AA()
        self.uni_l.calc_C()
        K_left, h_left_uni = self.uni_l.calc_K_l()
        
        self.K_l[0][:] = K_left[-1]

        for n in xrange(self.N, self.N_centre - 1, -1):
            self.K[n], he = tm.calc_K(self.K[n + 1], self.C[n], self.get_l(n - 1),
                                      self.r[n + 1], self.A[n], self.get_AA(n))
                
            self.h_expect[n] = he

        for n in xrange(1, self.N_centre + 1):
            self.K_l[n], he = tm.calc_K_l(self.K_l[n - 1], self.C[n - 1], self.get_l(n - 2),
                                          self.r[n], self.A[n], self.get_AA(n - 1))
                
            self.h_expect[n - 1] = he

        self.dH_expect = (mm.adot_noconj(self.K_l[self.N_centre], self.r[self.N_centre]) 
                          + mm.adot(self.l[self.N_centre - 1], self.K[self.N_centre]) 
                          - (self.N + 1) * self.uni_r.h_expect)

    def expect_3s(self, op, n):
        """Computes the expectation value of a nearest-neighbour three-site operator.

        The operator should be a q[n] x q[n + 1] x q[n + 2] x q[n] x
        q[n + 1] x q[n + 2] array such that op[s, t, u, v, w, x] =
        <stu|op|vwx> or a function of the form op(s, t, u, v, w, x) =
        <stu|op|vwx>.

        The state must be up-to-date -- see self.update()!

        Parameters
        ----------
        op : ndarray or callable
            The operator array or function.
        n : int
            The leftmost site number (operator acts on n, n + 1, n + 2).

        Returns
        -------
        expval : floating point number
            The expectation value (data type may be complex)
        """
        A = self.A[n]
        Ap1 = self.A[n + 1]
        Ap2 = self.A[n + 2]
        AAA = tm.calc_AAA(A, Ap1, Ap2)

        if callable(op):
            op = sp.vectorize(op, otypes=[sp.complex128])
            op = sp.fromfunction(op, (A.shape[0], Ap1.shape[0], Ap2.shape[0],
                                      A.shape[0], Ap1.shape[0], Ap2.shape[0]))

        C = tm.calc_C_3s_mat_op_AAA(op, AAA)
        res = tm.eps_r_op_3s_C123_AAA456(self.r[n + 2], C, AAA)
        return m.adot(self.l[n - 1], res)

    def expect_2s(self, op, n):
        """Computes the expectation value of a nearest-neighbour two-site operator.

        The operator should be a q[n] x q[n + 1] x q[n] x q[n + 1] array
        such that op[s, t, u, v] = <st|op|uv> or a function of the form
        op(s, t, u, v) = <st|op|uv>.

        Parameters
        ----------
        o : ndarray or callable
            The operator array or function.
        n : int
            The leftmost site number (operator acts on n, n + 1).
        """
        A = self.get_A(n)
        Ap1 = self.get_A(n + 1)
        AA = tm.calc_AA(A, Ap1)

        if callable(op):
            op = sp.vectorize(op, otypes=[sp.complex128])
            op = sp.fromfunction(op, (A.shape[0], Ap1.shape[0], A.shape[0], Ap1.shape[0]))

        C = tm.calc_C_mat_op_AA(op, AA)
        res = tm.eps_r_op_2s_C12_AA34(self.get_r(n + 1), C, AA)
        return mm.adot(self.get_l(n - 1), res)

    def density_2s(self, n1, n2):
        """Returns a reduced density matrix for a pair of sites.
        
        Currently only supports sites in the nonuniform window.

        Parameters
        ----------
        n1 : int
            The site number of the first site.
        n2 : int
            The site number of the second site (must be > n1).
        """
        rho = sp.empty((self.q[n1] * self.q[n2], self.q[n1] * self.q[n2]), dtype=sp.complex128)
        r_n2 = sp.empty_like(self.r[n2 - 1])
        r_n1 = sp.empty_like(self.r[n1 - 1])
        ln1m1 = self.get_l(n1 - 1)

        for s2 in xrange(self.q[n2]):
            for t2 in xrange(self.q[n2]):
                r_n2 = mm.mmul(self.A[n2][t2], self.r[n2], mm.H(self.A[n2][s2]))

                r_n = r_n2
                for n in reversed(xrange(n1 + 1, n2)):
                    r_n = tm.eps_r_noop(r_n, self.A[n], self.A[n])

                for s1 in xrange(self.q[n1]):
                    for t1 in xrange(self.q[n1]):
                        r_n1 = mm.mmul(self.A[n1][t1], r_n, mm.H(self.A[n1][s1]))
                        tmp = mm.adot(ln1m1, r_n1)
                        rho[s1 * self.q[n1] + s2, t1 * self.q[n1] + t2] = tmp
        return rho

 def matvec(self, v):
     x = v.reshape((self.D, self.D))
     
     if self.left:
         xE = self.tdvp._eps_l_noop_dense_A(x, self.out)
         QEQ = xE - m.H(self.l) * m.adot(self.r, x)
     else:
         Ex = self.tdvp._eps_r_noop_dense_A(x, self.out)
         QEQ = Ex - self.r * m.adot(self.l, x)        
     
     
     if not self.p == 0:
         QEQ *= np.exp(1.j * self.p)
     
     res = x - QEQ
     
     return res.ravel()

    def matvec(self, v):
        x = v.reshape((self.D, self.D))

        res = self.tmp
        out = self.out
        res[:] = x

        if self.left:  # Multiplying from the left, but x is a col. vector, so use mat_dagger
            for k in xrange(len(self.A1)):
                out = tm.eps_l_noop_inplace(res, self.A1[k], self.A2[k], out)
                tmp = res
                res = out
                out = tmp

            if self.pseudo:
                out[:] = self.lL
                out *= m.adot(self.rL, x)
                res -= out
                res *= -sp.exp(-1.0j * self.p)
                res += x
            else:
                res *= -sp.exp(-1.0j * self.p)
                res += x
        else:
            for k in xrange(len(self.A1) - 1, -1, -1):
                out = tm.eps_r_noop_inplace(res, self.A1[k], self.A2[k], out)
                tmp = res
                res = out
                out = tmp

            if self.pseudo:
                out[:] = self.rL
                out *= m.adot(self.lL, x)
                res -= out
                res *= -sp.exp(1.0j * self.p)
                res += x
            else:
                res *= -sp.exp(1.0j * self.p)
                res += x

        return (
            res.copy().ravel()
        )  # apparently we can't reuse the memory we return (seems to blow up lgmres), so we make a copy

 def matvec(self, v):
     x = v.reshape((self.D, self.D))
     
     if self.left: #Multiplying from the left, but x is a col. vector, so use mat_dagger
         Ehx = tm.eps_l_noop_inplace(x, self.A1, self.A2, self.out)
         if self.pseudo:
             QEQhx = Ehx - self.l * m.adot(self.r, x)
             res = x - sp.exp(-1.j * self.p) * QEQhx
         else:
             res = x - sp.exp(-1.j * self.p) * Ehx
     else:
         Ex = tm.eps_r_noop_inplace(x, self.A1, self.A2, self.out)
         if self.pseudo:
             QEQx = Ex - self.r * m.adot(self.l, x)
             res = x - sp.exp(1.j * self.p) * QEQx
         else:
             res = x - sp.exp(1.j * self.p) * Ex
     
     return res.ravel()

 def calc_K(self):
     Hr = np.zeros_like(self.A[0])
     
     for s in xrange(self.q):
         for t in xrange(self.q):
             Hr += m.mmul(self.C[s, t], self.r, m.H(self.AA[s, t]))
     
     self.h = m.adot(self.l, Hr)
     
     QHr = Hr - self.r * self.h
     
     self.calc_PPinv(QHr, out=self.K)
     
     if self.sanity_checks:
         Ex = self.eps_r(self.K)
         QEQ = Ex - self.r * m.adot(self.l, self.K)
         res = self.K - QEQ
         if not np.allclose(res, QHr):
             print "Sanity check failed: Bad K!"
             print "Off by: " + str(la.norm(res - QHr))

    def simple_renorm(self, update_lr=True):
        """Renormalizes the state in by multiplying one of the parameter
        tensors by a factor.
        """
        norm = mm.adot(self.l[self.N_centre - 1], self.r[self.N_centre - 1])

        if abs(1 - norm) > 1E-15:
            self.A[self.N_centre] *= 1 / sp.sqrt(norm)

        if update_lr:
            self.calc_l(n_low=self.N_centre)
            self.calc_r(n_high=self.N_centre)

 def f(tau, *args):
     if tau < 0:
         if use_tangvec_overlap:
             res = tau**2 + self.eta_sq.sum().real
         else:
             res = tau**2 + h_expect_0.real
         log.debug((tau, res, "punishing negative tau!"))
         taus.append(tau)
         ress.append(res)
         hs.append(h_expect_0.real)
         return res
     try:
         i = taus.index(tau)
         log.debug((tau, ress[i], "from stored"))
         return ress[i]
     except ValueError:
         for n in xrange(1, self.N + 1):
             if not Bs[n] is None:
                 self.A[n] = As0[n] - tau * Bs[n]
             
         if use_tangvec_overlap:
             self.update(restore_CF=False)
             Bsg = self.calc_B(set_eta=False)
             res = 0
             for n in xrange(1, self.N + 1):
                 if not Bs[n] is None:
                     res += abs(m.adot(self.l[n - 1], tm.eps_r_noop(self.r[n], Bsg[n], Bs[n])))
             h_exp = self.H_expect.real
         else:
             self.calc_l()
             self.calc_r()
             self.simple_renorm()
             self.calc_C()
             
             h_exp = 0
             if self.ham_sites == 2:
                 for n in xrange(1, self.N):
                     h_exp += self.expect_2s(self.ham[n], n).real
             else:
                 for n in xrange(1, self.N - 1):
                     h_exp += self.expect_3s(self.ham[n], n).real
             res = h_exp
         
         log.debug((tau, res, h_exp, h_exp - h_expect_0.real))
         
         taus.append(tau)
         ress.append(res)
         hs.append(h_exp)
         
         return res

 def matvec(self, v):
     x = v.reshape((self.D, self.D))
     
     if self.left: #Multiplying from the left, but x is a col. vector, so use mat_dagger
         Ehx = x
         for k in xrange(len(self.A1)):
             Ehx = tm.eps_l_noop(Ehx, self.A1[k], self.A2[k])
         if self.pseudo:
             QEQhx = Ehx - self.lL * m.adot(self.rL, x)
             res = x - sp.exp(-1.j * self.p) * QEQhx
         else:
             res = x - sp.exp(-1.j * self.p) * Ehx
     else:
         Ex = x
         for k in xrange(len(self.A1) - 1, -1, -1):
             Ex = tm.eps_r_noop(Ex, self.A1[k], self.A2[k])
         if self.pseudo:
             QEQx = Ex - self.rL * m.adot(self.lL, x)
             res = x - sp.exp(1.j * self.p) * QEQx
         else:
             res = x - sp.exp(1.j * self.p) * Ex
     
     return res.ravel()

 def calc_K_l(self):
     lH = np.zeros_like(self.A[0])
     
     for s in xrange(self.q):
         for t in xrange(self.q):
             lH += m.mmul(m.H(self.AA[s, t]), self.l, self.C[s, t])
     
     h = m.adot(lH, self.r)
     
     lHQ = lH - self.l * h
     
     self.K_left = self.calc_PPinv(lHQ, left=True, out=self.K_left)
     
     if self.sanity_checks:
         xE = self.eps_l(self.K_left)
         QEQ = xE - self.l * m.adot(self.K_left, self.r)
         res = self.K_left - QEQ
         if not np.allclose(res, lHQ):
             print "Sanity check failed: Bad K_left!"
             print "Off by: " + str(la.norm(res - lHQ))
     
     
     return self.K_left, h