#include "dmrg.jl"
function tdvptwosite(H::MPO, psi0::MPS, dt::Number, order::Int; kwargs...)::Tuple{Number,MPS}
    check_hascommoninds(siteinds, H, psi)
    check_hascommoninds(siteinds, H, psi')
    # Permute the indices to have a better memory layout
    # and minimize permutations
    H = permute(H, (linkind, siteinds, linkind))
    PH = ProjMPO(H)
    psi = copy(psi0)
    return tdvptwosite(PH, psi, dt, order; kwargs...)
  end


function tdvptwosite!(H::MPO, psi::MPS, dt::Number, order::Int; kwargs...)::Tuple{Number,MPS}
    check_hascommoninds(siteinds, H, psi)
    check_hascommoninds(siteinds, H, psi')
    # Permute the indices to have a better memory layout
    # and minimize permutations
    H = permute(H, (linkind, siteinds, linkind))
    PH = ProjMPO(H)
    return tdvptwosite(PH, psi, dt, order; kwargs...)
  end



function forward!(PH::ProjMPO, psi::MPS,dt::Number,kwargs)
    which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
    #svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
    svd_alg::String = get(kwargs, :svd_alg, "recursive")
    
    obs = get(kwargs, :observer, NoObserver())
    outputlevel::Int = get(kwargs, :outputlevel, 1)
  
    write_when_maxdim_exceeds::Union{Int,Nothing} = get(
      kwargs, :write_when_maxdim_exceeds, nothing
    )
  
    # eigsolve kwargs
    eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
    eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
    eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)
  
    # SVD kwargseigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    maxdim::Int = get(kwargs, :maxdim,100)
    cutoff::Float64 = get(kwargs, :cutoff,0.0)

    


    # TODO: add support for non-Hermitian DMRG
    ishermitian::Bool = get(kwargs, :ishermitian, true)
  

  
    # TODO: use this as preferred syntax for passing arguments
    # to eigsolve
    #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
    #                         verbosity = 0, ishermitian = true,
    #                         which_eigenvalue = :SR)
    #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)
  
    # Keyword argument deprecations
    if haskey(kwargs, :maxiter)
      error("""maxiter keyword has been replaced by eigsolve_krylovdim.
               Note: compared to the C++ version of ITensor,
               setting eigsolve_krylovdim 3 is the same as setting
               a maxiter of 2.""")
    end
  
    if haskey(kwargs, :errgoal)
      error("errgoal keyword has been replaced by eigsolve_tol.")
    end
  
    if haskey(kwargs, :quiet)
      error("quiet keyword has been replaced by outputlevel")
    end
  
    N = length(psi)    
    for (b, ha) in sweepnext(N)
        if ha == 2 
            return PH, psi
        end
        ##Taken care of in two steps below, assuming that before the sweep position was set
        #@timeit_debug timer "dmrg: position!" begin
        #  position!(PH, psi, b)
        #end

        @debug_check begin
            checkflux(psi)
            checkflux(PH)
        end
        
        @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
            phi = psi[b] * psi[b + 1]
        end

        @timeit_debug timer "tdvp2: expintegrator-2site" begin
            phi, propinfo = exponentiate(
            PH,
            dt,
            phi,
            ishermitian=ishermitian,
            tol=eigsolve_tol,
            krylovdim=eigsolve_krylovdim,
            #maxiter=eigsolve_maxiter,
            )
        end
        
        ortho = ha == 1 ? "left" : "right"
        
        @debug_check begin
            checkflux(phi)
        end
        
        @timeit_debug timer "dmrg: replacebond!" begin
            spec = replacebond!(
            psi,
            b,
            phi;
            maxdim=maxdim,
            cutoff=cutoff,
            #eigen_perturbation=drho,
            ortho=ortho,
            normalize=true,
            which_decomp=which_decomp,
            svd_alg=svd_alg,
            )
        end
        
        maxtruncerr = max(maxtruncerr, spec.truncerr)
        #@show inds(psi[b+1])
        
        
        if b+1<N 
            position_leftonly!(PH, psi, b+1)
            phi=psi[b + 1]
            #@show inds(psi[b+1])
            #@show inds(phi)
            @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                psi[b+1], propinfo = exponentiate(
                PH,
                -dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                #verbosity=1
                #maxiter=eigsolve_maxiter,
                )
            end
            #@show inds(vecs)
            position_rightonly!(PH, psi, b+1)
        else
            position!(PH,psi,b)
        end 

        @debug_check begin
            checkflux(psi)
            checkflux(PH)
        end
    end
  end

function backward!(PH::ProjMPO, psi::MPS,dt::Number,kwargs)
    which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
    #svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
    svd_alg::String = get(kwargs, :svd_alg, "recursive")
    
    obs = get(kwargs, :observer, NoObserver())
    outputlevel::Int = get(kwargs, :outputlevel, 1)
  
    write_when_maxdim_exceeds::Union{Int,Nothing} = get(
      kwargs, :write_when_maxdim_exceeds, nothing
    )
  
    # eigsolve kwargs
    eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
    eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
    eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)
  
    # SVD kwargseigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    maxdim::Int = get(kwargs, :maxdim,100)
    cutoff::Float64 = get(kwargs, :cutoff,0.0)

    


    # TODO: add support for non-Hermitian DMRG
    ishermitian::Bool = get(kwargs, :ishermitian, true)
  

  
    # TODO: use this as preferred syntax for passing arguments
    # to eigsolve
    #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
    #                         verbosity = 0, ishermitian = true,
    #                         which_eigenvalue = :SR)
    #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)
  
    # Keyword argument deprecations
    if haskey(kwargs, :maxiter)
      error("""maxiter keyword has been replaced by eigsolve_krylovdim.
               Note: compared to the C++ version of ITensor,
               setting eigsolve_krylovdim 3 is the same as setting
               a maxiter of 2.""")
    end
  
    if haskey(kwargs, :errgoal)
      error("errgoal keyword has been replaced by eigsolve_tol.")
    end
  
    if haskey(kwargs, :quiet)
      error("quiet keyword has been replaced by outputlevel")
    end
  
    N = length(psi)
    for (b, ha) in sweepnext(N)
        if ha == 2
            ortho = ha == 1 ? "left" : "right"
            if b+1<N
                phi=psi[b + 1]
                
            

                @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                    vecs, propinfo = exponentiate(
                    PH,
                    -dt,
                    phi,
                    ishermitian=ishermitian,
                    tol=eigsolve_tol,
                    krylovdim=eigsolve_krylovdim,
                    #verbosity=1
                    #maxiter=eigsolve_maxiter,
                    )
                end
                #@show inds(vecs)
                psi[b+1] = vecs
                position_leftonly!(PH, psi, b)
            end    
            #position_leftonly!(PH, psi, b)
            #position_rightonly!(PH, psi, b)
            
            @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
                phi = psi[b] * psi[b + 1]
            end
            

            @timeit_debug timer "tdvp2: expintegrator-2site" begin
                vecs, propinfo = exponentiate(
                PH,
                dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                #verbosity=1
                #maxiter=eigsolve_maxiter,
                )
            end
            #energy::Number = vals[1]
            phi = vecs
            @timeit_debug timer "dmrg: replacebond!" begin
                spec = replacebond!(
                psi,
                b,
                phi;
                maxdim=maxdim,
                cutoff=cutoff,
                #eigen_perturbation=drho,
                ortho=ortho,
                normalize=true,
                which_decomp=which_decomp,
                svd_alg=svd_alg,
                )
            end
            maxtruncerr = max(maxtruncerr, spec.truncerr)
            
            if b > 1    #after completed backward-forward sweep we want to be left with 2-site PH
                position_rightonly!(PH, psi, b-1)
                #position_leftonly!(PH, psi, b-1)
            else
                return PH,psi   
            end

        end
    end
  end
            
    
function tdvptwosite!(PH::ProjMPO, psi::MPS, dt::Number, order::Int; kwargs...)::Tuple{Number,MPS}
    if length(psi) == 1
      error(
        "`tdvp2` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.eigsolve`, etc.",
      )
    end
  
    @debug_check begin
      # Debug level checks
      # Enable with ITensors.enable_debug_checks()
      checkflux(psi)
      checkflux(PH)
    end
  
    which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
    #svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
    svd_alg::String = get(kwargs, :svd_alg, "recursive")
    
    obs = get(kwargs, :observer, NoObserver())
    outputlevel::Int = get(kwargs, :outputlevel, 1)
  
    write_when_maxdim_exceeds::Union{Int,Nothing} = get(
      kwargs, :write_when_maxdim_exceeds, nothing
    )
  
    # eigsolve kwargs
    eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
    eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
    eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)
  
    # SVD kwargseigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    maxdim::Int = get(kwargs, :maxdim,100)
    cutoff::Float64 = get(kwargs, :cutoff,0.0)

    


    # TODO: add support for non-Hermitian DMRG
    ishermitian::Bool = get(kwargs, :ishermitian, true)
  

  
    # TODO: use this as preferred syntax for passing arguments
    # to eigsolve
    #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
    #                         verbosity = 0, ishermitian = true,
    #                         which_eigenvalue = :SR)
    #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)
  
    # Keyword argument deprecations
    if haskey(kwargs, :maxiter)
      error("""maxiter keyword has been replaced by eigsolve_krylovdim.
               Note: compared to the C++ version of ITensor,
               setting eigsolve_krylovdim 3 is the same as setting
               a maxiter of 2.""")
    end
  
    if haskey(kwargs, :errgoal)
      error("errgoal keyword has been replaced by eigsolve_tol.")
    end
  
    if haskey(kwargs, :quiet)
      error("quiet keyword has been replaced by outputlevel")
    end
  
    N = length(psi)
  
    if !isortho(psi) || orthocenter(psi) != 1
      orthogonalize!(psi, 1)
    end
    @assert isortho(psi) && orthocenter(psi) == 1
  
    position!(PH, psi, 1)
    energy = 0.0
    if order==2
        sweeps=1
    end
    for sw in 1:sweeps
      sw_time = @elapsed begin
        maxtruncerr = 0.0
        #=
        if !isnothing(write_when_maxdim_exceeds) &&
          maxdim(sweeps, sw) > write_when_maxdim_exceeds
          if outputlevel >= 2
            println(
              "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sw)), writing environment tensors to disk",
            )
          end
          PH = disk(PH)
        end
        =#
        print("forward\n")
        PH,psi=forward!(PH,psi,dt/2.,kwargs)

        print("backward\n")
        PH,psi=backward!(PH,psi,dt/2.,kwargs)
        #=
        for (b, ha) in sweepnext(N)
          @debug_check begin
            checkflux(psi)
            checkflux(PH)
          end
          
          if ha == 1
            ##Taken care of in two steps below, assuming that before the sweep position was set
            #@timeit_debug timer "dmrg: position!" begin
            #  position!(PH, psi, b)
            #end
    
            @debug_check begin
                checkflux(psi)
                checkflux(PH)
            end
            
            @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
                phi = psi[b] * psi[b + 1]
            end
            #@show inds(phi)
            
            #@show phi
            @timeit_debug timer "tdvp2: expintegrator-2site" begin
                vecs, propinfo = exponentiate(
                PH,
                dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                
                #maxiter=eigsolve_maxiter,
                )
            end
            #print(size(PH),size(vecs))
            #energy::Number = vals[1]
            #print(vecs)
            phi = vecs
            #@show inds(phi);
            ortho = ha == 1 ? "left" : "right"
    
            #drho = nothing
            #if noise(sweeps, sw) > 0.0
            #  @timeit_debug timer "dmrg: noiseterm" begin
            #    # Use noise term when determining new MPS basis
            #    drho = noise(sweeps, sw) * noiseterm(PH, phi, ortho)
            #  end
            #end
    
            @debug_check begin
                checkflux(phi)
            end
            
            #for j in 1:N
            #    @show inds(psi[j])
            #end
            @timeit_debug timer "dmrg: replacebond!" begin
                spec = replacebond!(
                psi,
                b,
                phi;
                maxdim=maxdim,
                cutoff=cutoff,
                #eigen_perturbation=drho,
                ortho=ortho,
                normalize=true,
                which_decomp=which_decomp,
                svd_alg=svd_alg,
                )
            end
            #for j in 1:N
            #    @show inds(psi[j])
            #end
            maxtruncerr = max(maxtruncerr, spec.truncerr)
            #@show inds(psi[b+1])
            
            
            if b+1<N 
                position_leftonly!(PH, psi, b+1)
                phi=psi[b + 1]
                #@show inds(psi[b+1])
                #@show inds(phi)
                @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                    vecs, propinfo = exponentiate(
                    PH,
                    -dt,
                    phi,
                    ishermitian=ishermitian,
                    tol=eigsolve_tol,
                    krylovdim=eigsolve_krylovdim,
                    #verbosity=1
                    #maxiter=eigsolve_maxiter,
                    )
                end
                #@show inds(vecs)
                psi[b+1]=vecs        ##FIXME: necessary? was the prior assignment occurence of this assignment to new var?
                position_rightonly!(PH, psi, b+1)
            else
                position!(PH,psi,b)
            end
                
            
            @debug_check begin
                checkflux(psi)
                checkflux(PH)
            end

            
          elseif ha == 2
            #return (energy, psi)
            ortho = ha == 1 ? "left" : "right"
            if b+1<N
                phi=psi[b + 1]
                #@show inds(phi)
                #println(b,' ', ha)
                #println(indices(phi))
                #println(indices(PH))
            

                @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                    vecs, propinfo = exponentiate(
                    PH,
                    -dt,
                    phi,
                    ishermitian=ishermitian,
                    tol=eigsolve_tol,
                    krylovdim=eigsolve_krylovdim,
                    #verbosity=1
                    #maxiter=eigsolve_maxiter,
                    )
                end
                #@show inds(vecs)
                psi[b+1] = vecs
                position_leftonly!(PH, psi, b)
            end    
            #position_leftonly!(PH, psi, b)
            #position_rightonly!(PH, psi, b)
            
            @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
                phi = psi[b] * psi[b + 1]
            end
            
        
            @timeit_debug timer "tdvp2: expintegrator-2site" begin
                vecs, propinfo = exponentiate(
                PH,
                dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                #verbosity=1
                #maxiter=eigsolve_maxiter,
                )
            end
            #energy::Number = vals[1]
            phi = vecs
            @timeit_debug timer "dmrg: replacebond!" begin
                spec = replacebond!(
                psi,
                b,
                phi;
                maxdim=maxdim,
                cutoff=cutoff,
                #eigen_perturbation=drho,
                ortho=ortho,
                normalize=true,
                which_decomp=which_decomp,
                svd_alg=svd_alg,
                )
            end
            maxtruncerr = max(maxtruncerr, spec.truncerr)
            
            if b > 1    #after completed backward-forward sweep we want to be left with 2-site PH
                position_rightonly!(PH, psi, b-1)
                #position_leftonly!(PH, psi, b-1)
                
            end

            end
          if outputlevel >= 2
            @printf(
              "Sweep %d, half %d, bond (%d,%d) energy=%.12f\n", sw, ha, b, b + 1, energy
            )
            @printf(
              "  Truncated using cutoff=%.1E maxdim=%d mindim=%d\n",
              cutoff(sweeps, sw),
              maxdim(sweeps, sw),
              mindim(sweeps, sw)
            )
            @printf(
              "  Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b))
            )
            flush(stdout)
          end
          
          sweep_is_done = (b == 1 && ha == 2)
          =#
          #=
          measure!(
            obs;
            energy=energy,
            psi=psi,
            bond=b,
            sweep=sw,
            half_sweep=ha,
            spec=spec,
            outputlevel=outputlevel,
            sweep_is_done=sweep_is_done,
          )
          =#
        end
      end
      if outputlevel >= 1
        #=
        @printf(
          "After sweep %d energy=%.12f maxlinkdim=%d maxerr=%.2E time=%.3f\n",
          sw,
          energy,
          maxlinkdim(psi),
          maxtruncerr,
          sw_time
        )
        =#
        flush(stdout)
      end

    
    return (energy, psi)
  end