leading_boundary cleanup

Author: lkdvos

src/MPSKit.jl

@@ -147,6 +147,7 @@ include("algorithms/excitation/dmrgexcitation.jl")
 include("algorithms/excitation/chepigaansatz.jl")
 include("algorithms/excitation/exci_transfer_system.jl")
 
+include("algorithms/statmech/leading_boundary.jl")
 include("algorithms/statmech/vumps.jl")
 include("algorithms/statmech/vomps.jl")
 include("algorithms/statmech/gradient_grassmann.jl")

src/algorithms/statmech/gradient_grassmann.jl

@@ -1,14 +1,8 @@
-function leading_boundary(state::InfiniteMPS, H::DenseMPO, alg::GradientGrassmann,
-                          envs=environments(state, H))
-    (multi, envs, err) = leading_boundary(convert(MultilineMPS, state),
-                                          convert(MultilineMPO, H), alg, envs)
-    state = convert(InfiniteMPS, multi)
-    return (state, envs, err)
-end
-
-function leading_boundary(state::MultilineMPS, H, alg::GradientGrassmann,
-                          envs=environments(state, H))
-    fg(x) = GrassmannMPS.fg(x, H, envs)
+function leading_boundary(state::MultilineMPS,
+                          operator::MultilineMPO,
+                          alg::GradientGrassmann,
+                          envs::MultilineEnvironments=environments(state, H))
+    fg(x) = GrassmannMPS.fg(x, operator, envs)
     x, _, _, _, normgradhistory = optimize(fg, state,
                                            alg.method;
                                            (transport!)=GrassmannMPS.transport!,

src/algorithms/statmech/leading_boundary.jl

@@ -0,0 +1,33 @@
+
+@doc """
+    leading_boundary(ψ, O, [environments]; kwargs...)
+    leading_boundary(ψ, O, algorithm, environments)
+
+Compute the leading boundary MPS for operator `O` with initial guess `ψ`. If not specified, an
+optimization algorithm will be attempted based on the supplied keywords.
+
+## Arguments
+- `ψ::AbstractMPS`: initial guess
+- `O::AbstractMPO`: operator for which to find the leading_boundary
+- `[environments]`: MPS environment manager
+- `algorithm`: optimization algorithm
+
+## Keywords
+- `tol::Float64`: tolerance for convergence criterium
+- `maxiter::Int`: maximum amount of iterations
+- `verbosity::Int`: display progress information
+""" leading_boundary
+
+# TODO: alg selector
+
+# implementation always in terms of Multiline objects
+function leading_boundary(state::InfiniteMPS, operator::InfiniteMPO, alg,
+                          envs=environments(state, operator))
+    state_multi = convert(MultilineMPS, state)
+    operator_multi = convert(MultilineMPO, operator)
+    envs_multi = Multiline([envs])
+    state_multi′, envs_multi′, err = leading_boundary(state_multi, operator_multi, alg,
+                                                      envs_multi)
+    state′ = convert(InfiniteMPS, state_multi′)
+    return state′, envs, err
+end

src/environments/infinite_envs.jl

@@ -13,6 +13,8 @@ struct InfiniteEnvironments{V<:GenericMPSTensor} <: AbstractMPSEnvironments
     GRs::PeriodicVector{V}
 end
 
+Base.length(envs::InfiniteEnvironments) = length(envs.GLs)
+
 function environments(above::InfiniteMPS,
                       operator::Union{InfiniteMPO,InfiniteMPOHamiltonian},
                       below::InfiniteMPS=above;

src/utility/multiline.jl

@@ -28,6 +28,7 @@ function Base.axes(m::Multiline, i::Int)
     return i == 1 ? axes(parent(m), 1) :
            i == 2 ? axes(parent(m)[1], 1) : throw(ArgumentError("Invalid index $i"))
 end
+Base.eachindex(m::Multiline) = CartesianIndices(size(m))
 
 Base.getindex(m::Multiline, i::Int) = getindex(parent(m), i)
 Base.setindex!(m::Multiline, v, i::Int) = (setindex!(parent(m), v, i); m)