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 density_2s(self, n1, n2):
     """Returns a reduced density matrix for a pair of sites.
     
     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])
     
     for s2 in xrange(self.q[n2]):
         for t2 in xrange(self.q[n2]):
             r_n2 = m.mmul(self.A[n2][t2], self.r[n2], m.H(self.A[n2][s2]))
             
             r_n = r_n2
             for n in reversed(xrange(n1 + 1, n2)):
                 r_n = self.eps_r(n, r_n)        
                 
             for s1 in xrange(self.q[n1]):
                 for t1 in xrange(self.q[n1]):
                     r_n1 = m.mmul(self.A[n1][t1], r_n, m.H(self.A[n1][s1]))
                     tmp = m.mmul(self.l[n1 - 1], r_n1)
                     rho[s1 * self.q[n1] + s2, t1 * self.q[n1] + t2] = tmp.trace()
     return rho

    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 restore_ONR_n(self, n, G_n_i):
        """Transforms a single A[n] to obtain right orthonormalization.

        Implements the condition for right-orthonormalization from sub-section
        3.1, theorem 1 of arXiv:quant-ph/0608197v2.

        This function must be called for each n in turn, starting at N + 1,
        passing the gauge transformation matrix from the previous step
        as an argument.

        Finds a G[n-1] such that ON_R is fulfilled for n.

        Eigenvalues = 0 are a problem here... IOW rank-deficient matrices.
        Apparently, they can turn up during a run, but if they do we're screwed.

        The fact that M should be positive definite is used to optimize this.

        Parameters
        ----------
        n : int
            The site number.
        G_n_i : ndarray
            The inverse gauge transform matrix for site n obtained in the previous step (for n + 1).

        Returns
        -------
        G_n_m1_i : ndarray
            The inverse gauge transformation matrix for the site n - 1.
        """
        if G_n_i is None:
            GGh_n_i = self.r[n]
        else:
            GGh_n_i = mm.mmul(G_n_i, self.r[n], mm.H(G_n_i))

        M = self.eps_r(n, GGh_n_i)

        try:
            tu = la.cholesky(M) #Assumes M is pos. def.. It should raise LinAlgError if not.
            G_nm1 = mm.H(mm.invtr(tu)) #G is now lower-triangular
            G_nm1_i = mm.H(tu)
        except sp.linalg.LinAlgError:
            print "Restore_ON_R_%u: Falling back to eigh()!" % n
            e,Gh = la.eigh(M)
            G_nm1 = mm.H(mm.mmul(Gh, sp.diag(1/sp.sqrt(e) + 0.j)))
            G_nm1_i = la.inv(G_nm1)

        if G_n_i is None:
            G_n_i = G_nm1_i

        if self.sanity_checks:
            if not sp.allclose(sp.dot(G_nm1, G_nm1_i), sp.eye(G_nm1.shape[0]), atol=1E-13, rtol=1E-13):
                print "Sanity Fail in restore_ONR_%u!: Bad GT at n=%u" % (n, n)

        for s in xrange(self.q[n]):
            self.A[n][s] = mm.mmul(G_nm1, self.A[n][s], G_n_i)

        return G_nm1_i, G_nm1

    def restore_ONR_n(self, n, G_n_i):
        """Transforms a single A[n] to obtain right orthonormalization.
        
        Implements the condition for right-orthonormalization from sub-section
        3.1, theorem 1 of arXiv:quant-ph/0608197v2.
        
        This function must be called for each n in turn, starting at N + 1,
        passing the gauge transformation matrix from the previous step
        as an argument.
        
        Finds a G[n-1] such that ON_R is fulfilled for n.
        
        Eigenvalues = 0 are a problem here... IOW rank-deficient matrices. 
        Apparently, they can turn up during a run, but if they do we're screwed.    
        
        The fact that M should be positive definite is used to optimize this.
        
        Parameters
        ----------
        n : int
            The site number.
        G_n_i : ndarray
            The inverse gauge transform matrix for site n obtained in the previous step (for n + 1).
    
        Returns
        -------
        G_n_m1_i : ndarray
            The inverse gauge transformation matrix for the site n - 1.
        """
        GGh_n_i = m.mmul(
            G_n_i, m.H(G_n_i)
        )  #r[n] does not belong here. The condition is for sum(AA). r[n] = 1 is a consequence.

        M = self.eps_r(n, GGh_n_i)

        #The following should be more efficient than eigh():
        try:
            tu = la.cholesky(
                M
            )  #Assumes M is pos. def.. It should raise LinAlgError if not.
            G_nm1 = m.H(m.invtr(tu))  #G is now lower-triangular
            G_nm1_i = m.H(tu)
        except sp.linalg.LinAlgError:
            print "restore_ONR_n: Falling back to eigh()!"
            e, Gh = la.eigh(M)
            G_nm1 = m.H(m.mmul(Gh, sp.diag(1 / sp.sqrt(e) + 0.j)))
            G_nm1_i = la.inv(G_nm1)

        for s in xrange(self.q[n]):
            self.A[n][s] = m.mmul(G_nm1, self.A[n][s], G_n_i)
            #It's ok to use the same matrix as out and as an operand here
            #since there are > 2 matrices in the chain and it is not the last argument.

        return G_nm1_i

 def calc_x(self, n, Vsh, sqrt_l, sqrt_r, sqrt_l_inv, sqrt_r_inv):
     """Calculate the parameter matrix x* giving the desired B.
     
     This is equivalent to eqn. (49) of arXiv:1103.0936v2 [cond-mat.str-el] except 
     that, here, norm-preservation is not enforced, such that the optimal 
     parameter matrices x*_n (for the parametrization of B) are given by the 
     derivative w.r.t. x_n of <Phi[B, A]|Ĥ|Psi[A]>, rather than 
     <Phi[B, A]|Ĥ - H|Psi[A]> (with H = <Psi|Ĥ|Psi>).
     
     An additional sum was added for the single-site hamiltonian.
     
     Some multiplications have been pulled outside of the sums for efficiency.
     
     Direct dependencies: 
         - A[n - 1], A[n], A[n + 1]
         - r[n], r[n + 1], l[n - 2], l[n - 1]
         - C[n], C[n - 1]
         - K[n + 1]
         - V[n]
     """
     x = sp.zeros((self.D[n - 1], self.q[n] * self.D[n] - self.D[n - 1]), dtype=self.typ, order=self.odr)
     x_part = sp.empty_like(x)
     x_subpart = sp.empty_like(self.A[n][0])
     x_subsubpart = sp.empty_like(self.A[n][0])
     
     x_part.fill(0)
     for s in xrange(self.q[n]):
         x_subpart.fill(0)    
         
         if n < self.N:
             x_subsubpart.fill(0)
             for t in xrange(self.q[n + 1]):
                 x_subsubpart += m.mmul(self.C[n][s,t], self.r[n + 1], m.H(self.A[n + 1][t])) #~1st line
                 
             x_subsubpart += m.mmul(self.A[n][s], self.K[n + 1]) #~3rd line               
             
             x_subpart += m.mmul(x_subsubpart, sqrt_r_inv)
         
         if not self.h_ext is None:
             x_subsubpart.fill(0)
             for t in xrange(self.q[n]):                         #Extra term to take care of h_ext..
                 x_subsubpart += self.h_ext(n, s, t) * self.A[n][t] #it may be more effecient to squeeze this into the nn term...
             x_subpart += m.mmul(x_subsubpart, sqrt_r)
         
         x_part += m.mmul(x_subpart, Vsh[s])
             
     x += m.mmul(sqrt_l, x_part)
         
     if n > 1:
         x_part.fill(0)
         for s in xrange(self.q[n]):     #~2nd line
             x_subsubpart.fill(0)
             for t in xrange(self.q[n + 1]):
                 x_subsubpart += m.mmul(m.H(self.A[n - 1][t]), self.l[n - 2], self.C[n - 1][t, s])
             x_part += m.mmul(x_subsubpart, sqrt_r, Vsh[s])
         x += m.mmul(sqrt_l_inv, x_part)
             
     return x

    def restore_RCF_l(self):
        G_nm1 = None
        l_nm1 = self.l[0]
        for n in xrange(self.N + 1):
            if n == 0:
                x = l_nm1
            else:
                x = mm.mmul(mm.H(G_nm1), l_nm1, G_nm1)
            M = self.eps_l(n, x)
            ev, EV = la.eigh(M)

            self.l[n] = mm.simple_diag_matrix(ev, dtype=self.typ)
            G_n_i = EV

            if n == 0:
                G_nm1 = mm.H(EV) #for left uniform case
                l_nm1 = self.l[n] #for sanity check
                self.u_gnd_l.r = mm.mmul(G_nm1, self.u_gnd_l.r, G_n_i) #since r is not eye

            for s in xrange(self.q[n]):
                self.A[n][s] = mm.mmul(G_nm1, self.A[n][s], G_n_i)

            if self.sanity_checks:
                l = self.eps_l(n, l_nm1)
                if not sp.allclose(l, self.l[n], atol=1E-12, rtol=1E-12):
                    print "Sanity Fail in restore_RCF_l!: l_%u is bad" % n
                    print la.norm(l - self.l[n])

            G_nm1 = mm.H(EV)
            l_nm1 = self.l[n]

            if self.sanity_checks:
                if not sp.allclose(sp.dot(G_nm1, G_n_i), sp.eye(G_n_i.shape[0]),
                                   atol=1E-12, rtol=1E-12):
                    print "Sanity Fail in restore_RCF_l!: Bad GT for l_%u" % n

        #Now G_nm1 = G_N
        G_nm1_i = mm.H(G_nm1)
        for s in xrange(self.q[self.N + 1]):
            self.A[self.N + 1][s] = mm.mmul(G_nm1, self.A[self.N + 1][s], G_nm1_i)

        ##This should not be necessary if G_N is really unitary
        #self.r[self.N] = mm.mmul(G_nm1, self.r[self.N], mm.H(G_nm1))
        #self.r[self.N + 1] = self.r[self.N]
        self.u_gnd_r.l[:] = mm.mmul(mm.H(G_nm1_i), self.u_gnd_r.l, G_nm1_i)
        
        self.S_hc = sp.zeros((self.N), dtype=sp.complex128)
        for n in xrange(1, self.N + 1):
            self.S_hc[n-1] = -sp.sum(self.l[n].diag * sp.log2(self.l[n].diag))

    def restore_SCF(self):
        X = la.cholesky(self.r, lower=True)
        Y = la.cholesky(self.l, lower=False)
        
        U, sv, Vh = la.svd(Y.dot(X))
        
        #s contains the Schmidt coefficients,
        lam = sv**2
        self.S_hc = - np.sum(lam * sp.log2(lam))
        
        S = m.simple_diag_matrix(sv, dtype=self.typ)
        Srt = S.sqrt()
        
        g = m.mmul(Srt, Vh, m.invtr(X, lower=True))
        
        g_i = m.mmul(m.invtr(Y, lower=False), U, Srt)
        
        for s in xrange(self.q):
            self.A[s] = m.mmul(g, self.A[s], g_i)
                
        if self.sanity_checks:
            Sfull = np.asarray(S)
            
            if not np.allclose(g.dot(g_i), np.eye(self.D)):
                print "Sanity check failed! Restore_SCF, bad GT!"
            
            l = m.mmul(m.H(g_i), self.l, g_i)
            r = m.mmul(g, self.r, m.H(g))
            
            if not np.allclose(Sfull, l):
                print "Sanity check failed: Restorce_SCF, left failed!"
                
            if not np.allclose(Sfull, r):
                print "Sanity check failed: Restorce_SCF, right failed!"
                
            l = self.eps_l(Sfull)
            r = self.eps_r(Sfull)
            
            if not np.allclose(Sfull, l, rtol=self.itr_rtol*self.check_fac, 
                               atol=self.itr_atol*self.check_fac):
                print "Sanity check failed: Restorce_SCF, left bad!"
                
            if not np.allclose(Sfull, r, rtol=self.itr_rtol*self.check_fac, 
                               atol=self.itr_atol*self.check_fac):
                print "Sanity check failed: Restorce_SCF, right bad!"

        self.l = S
        self.r = S

    def calc_Vsh(self, n, sqrt_r):
        """Generates m.H(V[n][s]) for a given n, used for generating B[n][s]
        
        This is described on p. 14 of arXiv:1103.0936v2 [cond-mat.str-el] for left 
        gauge fixing. Here, we are using right gauge fixing.
        
        Array slicing and reshaping is used to manipulate the indices as necessary.
        
        Each V[n] directly depends only on A[n] and r[n].
        
        We return the conjugate m.H(V) because we use it in more places than V.
        """
        R = sp.zeros((self.D[n], self.q[n], self.D[n-1]), dtype=self.typ, order='C')
        
        for s in xrange(self.q[n]):
            R[:,s,:] = m.mmul(sqrt_r, m.H(self.A[n][s]))

        R = R.reshape((self.q[n] * self.D[n], self.D[n-1]))
        V = m.H(ns.nullspace_qr(m.H(R)))
        #print (q[n]*D[n] - D[n-1], q[n]*D[n])
        #print V.shape
        #print sp.allclose(mat(V) * mat(V).H, sp.eye(q[n]*D[n] - D[n-1]))
        #print sp.allclose(mat(V) * mat(Rh).H, 0)
        V = V.reshape((self.q[n] * self.D[n] - self.D[n - 1], self.D[n], self.q[n])) #this works with the above form for R
        
        #prepare for using V[s] and already take the adjoint, since we use it more often
        Vsh = sp.empty((self.q[n], self.D[n], self.q[n] * self.D[n] - self.D[n - 1]), dtype=self.typ, order=self.odr)
        for s in xrange(self.q[n]):
            Vsh[s] = m.H(V[:,:,s])
        
        return Vsh

    def eps_l(self, n, x):
        """Implements the left epsilon map

        Parameters
        ----------
        res : ndarray
            A matrix to hold the result (with the same dimensions as l[n]). May be None.
        n : int
            The site number.
        x : ndarray
            The argument matrix. For example, using l[n - 1] gives a result l[n]

        Returns
        -------
        res : ndarray
            The resulting matrix.
        """
        if n > self.N + 1:
            n = self.N + 1
        elif n < 0:
            n = 0

        res = sp.zeros_like(self.l[n])

        for s in xrange(self.q[n]):
            res += mm.mmul(mm.H(self.A[n][s]), x, self.A[n][s])
        return res

    def calc_C(self, n_low=-1, n_high=-1):
        """Generates the C matrices used to calculate the K's and ultimately the B's
        
        These are to be used on one side of the super-operator when applying the
        nearest-neighbour Hamiltonian, similarly to C in eqn. (44) of 
        arXiv:1103.0936v2 [cond-mat.str-el], except being for the non-norm-preserving case.

        Makes use only of the nearest-neighbour hamiltonian, and of the A's.
        
        C[n] depends on A[n] and A[n + 1].
        
        """
        if self.h_nn is None:
            return 0
        
        if n_low < 1:
            n_low = 1
        if n_high < 1:
            n_high = self.N
        
        for n in xrange(n_low, n_high):
            self.C[n].fill(0)
            for u in xrange(self.q[n]):
                for v in xrange(self.q[n + 1]):
                    AA = m.mmul(self.A[n][u], self.A[n + 1][v]) #only do this once for each 
                    for s in xrange(self.q[n]):
                        for t in xrange(self.q[n + 1]):                
                            h_nn_stuv = self.h_nn(n, s, t, u, v)
                            if h_nn_stuv != 0:
                                self.C[n][s, t] += h_nn_stuv * AA

    def expect_1s_cor(self, o1, o2, n1, n2):
        """Computes the correlation of two single site operators acting on two different sites.
        
        See expect_1s().
        
        n1 must be smaller than n2.
        
        Assumes that the state is normalized.
        
        Parameters
        ----------
        o1 : function
            The first operator, acting on the first site.
        o2 : function
            The second operator, acting on the second site.
        n1 : int
            The site number of the first site.
        n2 : int
            The site number of the second site (must be > n1).
        """        
        r_n = self.eps_r(n2, self.r[n2], o2)

        for n in reversed(xrange(n1 + 1, n2)):
            r_n = self.eps_r(n, r_n)

        r_n = self.eps_r(n1, r_n, o1)   
         
        res = m.mmul(self.l[n1 - 1], r_n)
        return res.trace()

    def eps_l(self, n, x, out=None):
        """Implements the left epsilon map
        
        FIXME: Ref.
        
        Parameters
        ----------
        n : int
            The site number.
        x : ndarray
            The argument matrix. For example, using l[n - 1] gives a result l[n]
        out : ndarray
            A matrix to hold the result (with the same dimensions as l[n]). May be None.
    
        Returns
        -------
        res : ndarray
            The resulting matrix.
        """
        if out is None:
            out = sp.zeros_like(self.l[n])
        else:
            out.fill(0.)

        for s in xrange(self.q[n]):
            out += m.mmul(m.H(self.A[n][s]), x, self.A[n][s])
        return out

    def calc_B(self, n, set_eta=True):
        """Generates the B[n] tangent vector corresponding to physical evolution of the state.

        In other words, this returns B[n][x*] (equiv. eqn. (47) of
        arXiv:1103.0936v2 [cond-mat.str-el])
        with x* the parameter matrices satisfying the Euler-Lagrange equations
        as closely as possible.
        
        In the case of Bc, use the general Bc generated in calc_B_centre().
        """
        if n == self.N_centre:
            B, eta_c = self.calc_B_centre()
            if set_eta:
                self.eta[self.N_centre] = eta_c
        else:
            l_sqrt, r_sqrt, l_sqrt_inv, r_sqrt_inv = self.calc_l_r_roots(n)
            
            if n > self.N_centre:
                Vsh = tm.calc_Vsh(self.A[n], r_sqrt, sanity_checks=self.sanity_checks)
                x = self.calc_x(n, Vsh, l_sqrt, r_sqrt, l_sqrt_inv, r_sqrt_inv, right=True)
                
                B = sp.empty_like(self.A[n])
                for s in xrange(self.q[n]):
                    B[s] = mm.mmul(l_sqrt_inv, x, mm.H(Vsh[s]), r_sqrt_inv)
                    
                if self.sanity_checks:
                    M = tm.eps_r_noop(self.r[n], B, self.A[n])
                    if not sp.allclose(M, 0):
                        print "Sanity Fail in calc_B!: B_%u does not satisfy GFC!" % n
            else:
                Vsh = tm.calc_Vsh_l(self.A[n], l_sqrt, sanity_checks=self.sanity_checks)
                x = self.calc_x(n, Vsh, l_sqrt, r_sqrt, l_sqrt_inv, r_sqrt_inv, right=False)
                
                B = sp.empty_like(self.A[n])
                for s in xrange(self.q[n]):
                    B[s] = mm.mmul(l_sqrt_inv, mm.H(Vsh[s]), x, r_sqrt_inv)
                    
                if self.sanity_checks:
                    M = tm.eps_l_noop(self.l[n - 1], B, self.A[n])
                    if not sp.allclose(M, 0):
                        print "Sanity Fail in calc_B!: B_%u does not satisfy GFC!" % n
            
            if set_eta:
                self.eta[n] = sp.sqrt(mm.adot(x, x))


        return B

 def get_B_from_x(self, x, Vsh, l_sqrt_i, r_sqrt_i, out=None):
     if out is None:
         out = np.zeros_like(self.A)
         
     for s in xrange(self.q):
         out[s] = m.mmul(l_sqrt_i, x, m.H(Vsh[s]), r_sqrt_i)
         
     return out

 def density_1s(self, n):
     """Returns a reduced density matrix for a single site.
     
     Parameters
     ----------
     n1 : int
         The site number.
     """
     rho = sp.empty((self.q[n], self.q[n]), dtype=sp.complex128)
                 
     r_n = self.r[n]
     r_nm1 = sp.empty_like(self.r[n - 1])
     for s in xrange(self.q[n]):
         for t in xrange(self.q[n]):
             r_nm1 = m.mmul(self.A[n][t], r_n, m.H(self.A[n][s]))                
             rho[s, t] = m.mmul(self.l[n - 1], r_nm1).trace()
     return rho

 def eps_l(self, x, out=None):
     if out is None:
         out = np.zeros_like(self.A[0])
     else:
         out.fill(0.)
         
     for s in xrange(self.q):
         out += m.mmul(m.H(self.A[s]), x, self.A[s])        
         
     return out

    def calc_C(self, n_low=-1, n_high=-1):
        """Generates the C matrices used to calculate the K's and ultimately the B's

        These are to be used on one side of the super-operator when applying the
        nearest-neighbour Hamiltonian, similarly to C in eqn. (44) of
        arXiv:1103.0936v2 [cond-mat.str-el], except being for the non-norm-preserving case.

        Makes use only of the nearest-neighbour hamiltonian, and of the A's.

        C[n] depends on A[n] and A[n + 1].
        
        This calculation can be significantly faster if a matrix form for h_nn
        is available. See gen_h_matrix().

        """
        if self.h_nn is None:
            return 0

        if n_low < 1:
            n_low = 0
        if n_high < 1:
            n_high = self.N + 1
        
        if self.h_nn_mat is None:
            for n in xrange(n_low, n_high):
                self.C[n].fill(0)
                for u in xrange(self.q[n]):
                    for v in xrange(self.q[n + 1]):
                        AA = mm.mmul(self.A[n][u], self.A[n + 1][v]) #only do this once for each
                        for s in xrange(self.q[n]):
                            for t in xrange(self.q[n + 1]):
                                h_nn_stuv = self.h_nn(n, s, t, u, v)
                                if h_nn_stuv != 0:
                                    self.C[n][s, t] += h_nn_stuv * AA
        else:
            dot = sp.dot
            for n in xrange(n_low, n_high):
                An = self.A[n]
                Anp1 = self.A[n + 1]
                
                AA = sp.empty_like(self.C[n])
                for u in xrange(self.q[n]):
                    for v in xrange(self.q[n + 1]):
                        AA[u, v] = dot(An[u], Anp1[v])
                        
                if n == 0: #FIXME: Temp. hack
                    self.AA0 = AA
                elif n == 1:
                    self.AA1 = AA
                
                res = sp.tensordot(AA, self.h_nn_mat[n], ((0, 1), (2, 3)))
                res = sp.rollaxis(res, 3)
                res = sp.rollaxis(res, 3)
                
                self.C[n][:] = res

    def calc_K(self, n_low=-1, n_high=-1):
        """Generates the 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-gauge-preserving case, and includes a single-site
        Hamiltonian term.
        
        K[1] is, assuming a normalized state, the expectation value H of Ĥ.
        
        Instead of an explicit single-site term here, one could also include the 
        single-site Hamiltonian in the nearest-neighbour term, which may be more 
        efficient.
        """
        if n_low < 1:
            n_low = 1
        if n_high < 1:
            n_high = self.N + 1
            
        for n in reversed(xrange(n_low, n_high)):
            self.K[n].fill(0)

            if n < self.N:
                for s in xrange(self.q[n]): 
                    for t in xrange(self.q[n+1]):
                        self.K[n] += m.mmul(self.C[n][s, t],
                                              self.r[n + 1], m.H(self.A[n+1][t]), 
                                              m.H(self.A[n][s]))
                    self.K[n] += m.mmul(self.A[n][s], self.K[n + 1], 
                                          m.H(self.A[n][s]))
            
            if not self.h_ext is None:
                for s in xrange(self.q[n]):
                    for t in xrange(self.q[n]):
                        h_ext_st = self.h_ext(n, s, t)
                        if h_ext_st != 0:
                            self.K[n] += h_ext_st * m.mmul(self.A[n][t], 
                                                    self.r[n], m.H(self.A[n][s]))

    def calc_l_r_roots(self, n):
        """Returns the matrix square roots (and inverses) needed to calculate B.

        Hermiticity of l[n] and r[n] is used to speed this up.
        If an exception occurs here, it is probably because these matrices
        are no longer Hermitian (enough).
        
        If l[n] or r[n] are diagonal or the identity, further optimizations are
        used.
        """
        try:
            l_sqrt = self.l[n - 1].sqrt()
        except AttributeError:
            l_sqrt, evd = mm.sqrtmh(self.l[n - 1], ret_evd=True)

        try:
            l_sqrt_inv = l_sqrt.inv()
        except AttributeError:
            l_sqrt_inv = mm.invmh(l_sqrt, evd=evd)

        try:
            r_sqrt = self.r[n].sqrt()
        except AttributeError:
            r_sqrt, evd =  mm.sqrtmh(self.r[n], ret_evd=True)

        try:
            r_sqrt_inv = r_sqrt.inv()
        except AttributeError:
            r_sqrt_inv = mm.invmh(r_sqrt, evd=evd)

        if self.sanity_checks:
            if not sp.allclose(mm.mmul(l_sqrt, l_sqrt), self.l[n - 1]):
                print "Sanity Fail in calc_l_r_roots: Bad l_sqrt_%u" % (n - 1)
            if not sp.allclose(mm.mmul(r_sqrt, r_sqrt), self.r[n]):
                print "Sanity Fail in calc_l_r_roots: Bad r_sqrt_%u" % (n)
            if not sp.allclose(mm.mmul(l_sqrt, l_sqrt_inv), sp.eye(l_sqrt.shape[0])):
                print "Sanity Fail in calc_l_r_roots: Bad l_sqrt_inv_%u" % (n - 1)
            if not sp.allclose(mm.mmul(r_sqrt, r_sqrt_inv), sp.eye(r_sqrt.shape[0])):
                print "Sanity Fail in calc_l_r_roots: Bad r_sqrt_inv_%u" % (n)

        return l_sqrt, r_sqrt, l_sqrt_inv, r_sqrt_inv