From eef7d0f1d83555883f9ee10c396b72ef9d50e20c Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Apr 2025 13:08:36 -0400 Subject: [PATCH 01/28] Delete and move some files to make space for new code --- .../alternating_update/alternating_update.jl | 143 ------------------ .../alternating_update/region_update.jl | 77 ---------- src/solvers/extract/extract.jl | 28 ---- src/solvers/insert/insert.jl | 46 ------ src/solvers/local_solvers/eigsolve.jl | 34 ----- src/solvers/local_solvers/exponentiate.jl | 34 ----- .../{ => previous_interfaces}/contract.jl | 0 src/solvers/{ => previous_interfaces}/dmrg.jl | 0 .../{ => previous_interfaces}/dmrg_x.jl | 0 .../{ => previous_interfaces}/linsolve.jl | 0 src/solvers/{ => previous_interfaces}/tdvp.jl | 0 11 files changed, 362 deletions(-) delete mode 100644 src/solvers/alternating_update/alternating_update.jl delete mode 100644 src/solvers/alternating_update/region_update.jl delete mode 100644 src/solvers/extract/extract.jl delete mode 100644 src/solvers/insert/insert.jl delete mode 100644 src/solvers/local_solvers/eigsolve.jl delete mode 100644 src/solvers/local_solvers/exponentiate.jl rename src/solvers/{ => previous_interfaces}/contract.jl (100%) rename src/solvers/{ => previous_interfaces}/dmrg.jl (100%) rename src/solvers/{ => previous_interfaces}/dmrg_x.jl (100%) rename src/solvers/{ => previous_interfaces}/linsolve.jl (100%) rename src/solvers/{ => previous_interfaces}/tdvp.jl (100%) diff --git a/src/solvers/alternating_update/alternating_update.jl b/src/solvers/alternating_update/alternating_update.jl deleted file mode 100644 index 69f965f0..00000000 --- a/src/solvers/alternating_update/alternating_update.jl +++ /dev/null @@ -1,143 +0,0 @@ -using ITensors: state -using NamedGraphs.GraphsExtensions: GraphsExtensions - -function alternating_update( - operator, - init_state::AbstractTTN; - nsweeps, # define default for each solver implementation - nsites, # define default for each level of solver implementation - updater, # this specifies the update performed locally - outputlevel=default_outputlevel(), - region_printer=default_region_printer, - sweep_printer=default_sweep_printer, - (sweep_observer!)=nothing, - (region_observer!)=nothing, - root_vertex=GraphsExtensions.default_root_vertex(init_state), - extracter_kwargs=(;), - extracter=default_extracter(), - updater_kwargs=(;), - inserter_kwargs=(;), - inserter=default_inserter(), - transform_operator_kwargs=(;), - transform_operator=default_transform_operator(), - kwargs..., -) - inserter_kwargs = (; inserter_kwargs..., kwargs...) - sweep_plans = default_sweep_plans( - nsweeps, - init_state; - root_vertex, - extracter, - extracter_kwargs, - updater, - updater_kwargs, - inserter, - inserter_kwargs, - transform_operator, - transform_operator_kwargs, - nsites, - ) - return alternating_update( - operator, - init_state, - sweep_plans; - outputlevel, - sweep_observer!, - region_observer!, - sweep_printer, - region_printer, - ) -end - -function alternating_update( - projected_operator, - init_state::AbstractTTN, - sweep_plans; - outputlevel=default_outputlevel(), - checkdone=default_checkdone(), # - (sweep_observer!)=nothing, - sweep_printer=default_sweep_printer,#? - (region_observer!)=nothing, - region_printer=default_region_printer, -) - state = copy(init_state) - @assert !isnothing(sweep_plans) - for which_sweep in eachindex(sweep_plans) - sweep_plan = sweep_plans[which_sweep] - sweep_time = @elapsed begin - for which_region_update in eachindex(sweep_plan) - state, projected_operator = region_update( - projected_operator, - state; - which_sweep, - sweep_plan, - region_printer, - (region_observer!), - which_region_update, - outputlevel, - ) - end - end - update_observer!( - sweep_observer!; state, which_sweep, sweep_time, outputlevel, sweep_plans - ) - !isnothing(sweep_printer) && - sweep_printer(; state, which_sweep, sweep_time, outputlevel, sweep_plans) - checkdone(; - state, - which_sweep, - outputlevel, - sweep_plan, - sweep_plans, - sweep_observer!, - region_observer!, - ) && break - end - return state -end - -function alternating_update(operator::AbstractTTN, init_state::AbstractTTN; kwargs...) - projected_operator = ProjTTN(operator) - return alternating_update(projected_operator, init_state; kwargs...) -end - -function alternating_update( - operator::AbstractTTN, init_state::AbstractTTN, sweep_plans; kwargs... -) - projected_operator = ProjTTN(operator) - return alternating_update(projected_operator, init_state, sweep_plans; kwargs...) -end - -#ToDo: Fix docstring. -""" - tdvp(Hs::Vector{MPO},init_state::MPS,t::Number; kwargs...) - tdvp(Hs::Vector{MPO},init_state::MPS,t::Number, sweeps::Sweeps; kwargs...) - -Use the time dependent variational principle (TDVP) algorithm -to compute `exp(t*H)*init_state` using an efficient algorithm based -on alternating optimization of the MPS tensors and local Krylov -exponentiation of H. - -This version of `tdvp` accepts a representation of H as a -Vector of MPOs, Hs = [H1,H2,H3,...] such that H is defined -as H = H1+H2+H3+... -Note that this sum of MPOs is not actually computed; rather -the set of MPOs [H1,H2,H3,..] is efficiently looped over at -each step of the algorithm when optimizing the MPS. - -Returns: -* `state::MPS` - time-evolved MPS -""" -function alternating_update( - operators::Vector{<:AbstractTTN}, init_state::AbstractTTN; kwargs... -) - projected_operators = ProjTTNSum(operators) - return alternating_update(projected_operators, init_state; kwargs...) -end - -function alternating_update( - operators::Vector{<:AbstractTTN}, init_state::AbstractTTN, sweep_plans; kwargs... -) - projected_operators = ProjTTNSum(operators) - return alternating_update(projected_operators, init_state, sweep_plans; kwargs...) -end diff --git a/src/solvers/alternating_update/region_update.jl b/src/solvers/alternating_update/region_update.jl deleted file mode 100644 index c741c82a..00000000 --- a/src/solvers/alternating_update/region_update.jl +++ /dev/null @@ -1,77 +0,0 @@ -function region_update( - projected_operator, - state; - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - region_printer, - (region_observer!), -) - (region, region_kwargs) = sweep_plan[which_region_update] - (; - extracter, - extracter_kwargs, - updater, - updater_kwargs, - inserter, - inserter_kwargs, - transform_operator, - transform_operator_kwargs, - internal_kwargs, - ) = region_kwargs - - # ToDo: remove orthogonality center on vertex for generality - # region carries same information - if !isnothing(transform_operator) - projected_operator = transform_operator( - state, projected_operator; outputlevel, transform_operator_kwargs... - ) - end - state, projected_operator, phi = extracter( - state, projected_operator, region; extracter_kwargs..., internal_kwargs - ) - # create references, in case solver does (out-of-place) modify PH or state - state! = Ref(state) - projected_operator! = Ref(projected_operator) - # args passed by reference are supposed to be modified out of place - phi, info = updater( - phi; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - updater_kwargs..., - internal_kwargs, - ) - state = state![] - projected_operator = projected_operator![] - # ToDo: implement noise term as updater - #drho = nothing - #ortho = "left" #i guess with respect to ordered vertices that's valid but may be cleaner to use next_region logic - #if noise > 0.0 && isforward(direction) - # drho = noise * noiseterm(PH, phi, ortho) # TODO: actually implement this for trees... - # so noiseterm is a solver - #end - #if isa(region, AbstractEdge) && - state, spec = inserter(state, phi, region; inserter_kwargs..., internal_kwargs) - all_kwargs = (; - which_region_update, - sweep_plan, - total_sweep_steps=length(sweep_plan), - end_of_sweep=(which_region_update == length(sweep_plan)), - state, - region, - which_sweep, - spec, - outputlevel, - info..., - region_kwargs..., - internal_kwargs..., - ) - update_observer!(region_observer!; all_kwargs...) - !(isnothing(region_printer)) && region_printer(; all_kwargs...) - return state, projected_operator -end diff --git a/src/solvers/extract/extract.jl b/src/solvers/extract/extract.jl deleted file mode 100644 index 1013d1bd..00000000 --- a/src/solvers/extract/extract.jl +++ /dev/null @@ -1,28 +0,0 @@ -# Here extract_local_tensor and insert_local_tensor -# are essentially inverse operations, adapted for different kinds of -# algorithms and networks. -# -# In the simplest case, exact_local_tensor contracts together a few -# tensors of the network and returns the result, while -# insert_local_tensors takes that tensor and factorizes it back -# apart and puts it back into the network. -# - -function default_extracter(state, projected_operator, region; internal_kwargs) - if isa(region, AbstractEdge) - # TODO: add functionality for orthogonalizing onto a bond so that can be called instead - vsrc, vdst = src(region), dst(region) - state = orthogonalize(state, vsrc) - left_inds = uniqueinds(state[vsrc], state[vdst]) - U, S, V = svd( - state[vsrc], left_inds; lefttags=tags(state, region), righttags=tags(state, region) - ) - state[vsrc] = U - local_tensor = S * V - else - state = orthogonalize(state, region) - local_tensor = prod(state[v] for v in region) - end - projected_operator = position(projected_operator, state, region) - return state, projected_operator, local_tensor -end diff --git a/src/solvers/insert/insert.jl b/src/solvers/insert/insert.jl deleted file mode 100644 index 01fb35bd..00000000 --- a/src/solvers/insert/insert.jl +++ /dev/null @@ -1,46 +0,0 @@ -# Here extract_local_tensor and insert_local_tensor -# are essentially inverse operations, adapted for different kinds of -# algorithms and networks. - -# TODO: use dense TTN constructor to make this more general. -function default_inserter( - state::AbstractTTN, - phi::ITensor, - region; - normalize=false, - maxdim=nothing, - mindim=nothing, - cutoff=nothing, - internal_kwargs, -) - state = copy(state) - spec = nothing - if length(region) == 2 - v = last(region) - e = edgetype(state)(first(region), last(region)) - indsTe = inds(state[first(region)]) - L, phi, spec = factorize(phi, indsTe; tags=tags(state, e), maxdim, mindim, cutoff) - state[first(region)] = L - else - v = only(region) - end - state[v] = phi - state = set_ortho_region(state, [v]) - normalize && (state[v] /= norm(state[v])) - return state, spec -end - -function default_inserter( - state::AbstractTTN, - phi::ITensor, - region::NamedEdge; - cutoff=nothing, - maxdim=nothing, - mindim=nothing, - normalize=false, - internal_kwargs, -) - state[dst(region)] *= phi - state = set_ortho_region(state, [dst(region)]) - return state, nothing -end diff --git a/src/solvers/local_solvers/eigsolve.jl b/src/solvers/local_solvers/eigsolve.jl deleted file mode 100644 index ed993d80..00000000 --- a/src/solvers/local_solvers/eigsolve.jl +++ /dev/null @@ -1,34 +0,0 @@ -using KrylovKit: eigsolve - -function eigsolve_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - internal_kwargs, - which_eigval=:SR, - ishermitian=true, - tol=1e-14, - krylovdim=3, - maxiter=1, - verbosity=0, - eager=false, -) - howmany = 1 - vals, vecs, info = eigsolve( - projected_operator![], - init, - howmany, - which_eigval; - ishermitian, - tol, - krylovdim, - maxiter, - verbosity, - eager, - ) - return vecs[1], (; info, eigvals=vals) -end diff --git a/src/solvers/local_solvers/exponentiate.jl b/src/solvers/local_solvers/exponentiate.jl deleted file mode 100644 index c70a91c5..00000000 --- a/src/solvers/local_solvers/exponentiate.jl +++ /dev/null @@ -1,34 +0,0 @@ -using KrylovKit: exponentiate - -function exponentiate_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - internal_kwargs, - krylovdim=30, - maxiter=100, - verbosity=0, - tol=1E-12, - ishermitian=true, - issymmetric=true, - eager=true, -) - (; time_step) = internal_kwargs - result, exp_info = exponentiate( - projected_operator![], - time_step, - init; - eager, - krylovdim, - maxiter, - verbosity, - tol, - ishermitian, - issymmetric, - ) - return result, (; info=exp_info) -end diff --git a/src/solvers/contract.jl b/src/solvers/previous_interfaces/contract.jl similarity index 100% rename from src/solvers/contract.jl rename to src/solvers/previous_interfaces/contract.jl diff --git a/src/solvers/dmrg.jl b/src/solvers/previous_interfaces/dmrg.jl similarity index 100% rename from src/solvers/dmrg.jl rename to src/solvers/previous_interfaces/dmrg.jl diff --git a/src/solvers/dmrg_x.jl b/src/solvers/previous_interfaces/dmrg_x.jl similarity index 100% rename from src/solvers/dmrg_x.jl rename to src/solvers/previous_interfaces/dmrg_x.jl diff --git a/src/solvers/linsolve.jl b/src/solvers/previous_interfaces/linsolve.jl similarity index 100% rename from src/solvers/linsolve.jl rename to src/solvers/previous_interfaces/linsolve.jl diff --git a/src/solvers/tdvp.jl b/src/solvers/previous_interfaces/tdvp.jl similarity index 100% rename from src/solvers/tdvp.jl rename to src/solvers/previous_interfaces/tdvp.jl From b15d04559e2f43518eefc7bc93a7f70e2495a1d9 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Apr 2025 22:07:55 -0400 Subject: [PATCH 02/28] Add redesigned solver codes --- src/solvers/adapters.jl | 76 +++++++++++++++++++++++++++++++ src/solvers/alternating_update.jl | 23 ++++++++++ src/solvers/eigsolve.jl | 40 ++++++++++++++++ src/solvers/iterators.jl | 70 ++++++++++++++++++++++++++++ src/solvers/region_plans.jl | 71 +++++++++++++++++++++++++++++ src/solvers/subspace.jl | 52 +++++++++++++++++++++ src/solvers/tdvp.jl | 36 +++++++++++++++ 7 files changed, 368 insertions(+) create mode 100644 src/solvers/adapters.jl create mode 100644 src/solvers/alternating_update.jl create mode 100644 src/solvers/eigsolve.jl create mode 100644 src/solvers/iterators.jl create mode 100644 src/solvers/region_plans.jl create mode 100644 src/solvers/subspace.jl create mode 100644 src/solvers/tdvp.jl diff --git a/src/solvers/adapters.jl b/src/solvers/adapters.jl new file mode 100644 index 00000000..7278d245 --- /dev/null +++ b/src/solvers/adapters.jl @@ -0,0 +1,76 @@ + +# +# RepeatIterator +# +# An "iterator of iterators", plugging in initialization +# arguments into an `init` function each repetition +# + +mutable struct RepeatIterator{Iterator} + iterator::Iterator + keyword_args::Vector{<:NamedTuple} +end + +function RepeatIterator(iter, n::Integer, kwargs::NamedTuple) + return RepeatIterator(iter, fill(kwargs, n)) +end + +# +# Version 1: output the initialized iterator each time +# +function Base.iterate(R::RepeatIterator) + inner_next = iterate(R.iterator) + isnothing(inner_next) && return nothing + (item, inner_state) = inner_next + return item, (1, inner_state) +end +function Base.iterate(R::RepeatIterator, state) + state = Base.iterate(keyword_args) + isnothing(state) && return nothing + (kwargs, which) = state + R.iterator = init(R.iterator; kwargs...) + return R, which +end + +## +## Version 2: output the initialized iterator each time +## +#function Base.iterate(R::RepeatIterator, state=1) +# state = Base.iterate(keyword_args) +# isnothing(state) && return nothing +# (kwargs, which) = state +# R.iterator = init(R.iterator; kwargs...) +# return R, which +#end + +# +# TupleRegionIterator +# +# Adapts outputs to be (region, region_kwargs) tuples +# +# More generic design? maybe just assuming RegionIterator +# or its outputs implement some interface function that +# generates each tuple? +# + +mutable struct TupleRegionIterator{RegionIter} + region_iterator::RegionIter +end + +region_iterator(T::TupleRegionIterator) = T.region_iterator + +function Base.iterate(T::TupleRegionIterator, which=1) + state = iterate(region_iterator(T), which) + isnothing(state) && return nothing + (current_region, region_kwargs) = current_region_plan(region_iterator(T)) + return (current_region, region_kwargs), last(state) +end + +""" + region_tuples(R::RegionIterator) + +The `region_tuples` adapter converts a RegionIterator into an +iterator which outputs a tuple of the form (current_region, current_region_kwargs) +at each step. +""" +region_tuples(R::RegionIterator) = TupleRegionIterator(R) diff --git a/src/solvers/alternating_update.jl b/src/solvers/alternating_update.jl new file mode 100644 index 00000000..1e2b1013 --- /dev/null +++ b/src/solvers/alternating_update.jl @@ -0,0 +1,23 @@ + +function sweep_callback(problem; outputlevel, sweep, kws...) + if outputlevel >= 1 + println("Done with sweep $sweep") + #print(" cpu_time=", round(sweep_time; digits=3)) + #println() + flush(stdout) + end + return problem +end + +function alternating_update(sweep_iterator; outputlevel=0, kwargs...) + for (sweep, region_iter) in enumerate(sweep_iterator) + prob = problem(region_iter) + for (region, region_kwargs) in region_tuples(region_iter) + end + prob = sweep_callback( + prob; nsweeps=length(sweep_iterator), outputlevel, sweep, kwargs... + ) + #TODO need to update `problem` reference held within region iterator(s) + end + return problem(last(sweep_iterator)) +end diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl new file mode 100644 index 00000000..28e6903d --- /dev/null +++ b/src/solvers/eigsolve.jl @@ -0,0 +1,40 @@ +import ITensorNetworks as itn +import ITensors as it +using Printf + +@kwdef mutable struct EigsolveProblem{State,Operator} + state::State + operator::Operator + eigenvalue::Number = Inf +end + +eigenvalue(E::EigsolveProblem) = E.eigenvalue +state(E::EigsolveProblem) = E.state +operator(E::EigsolveProblem) = E.operator + +function set( + E::EigsolveProblem; state=state(E), operator=operator(E), eigenvalue=eigenvalue(E) +) + return EigsolveProblem(; state, operator, eigenvalue) +end + +function updater!(E::EigsolveProblem, local_tensor, region; outputlevel, kws...) + E.eigenvalue, local_tensor = eigsolve_updater(operator(E), local_tensor; kws...) + if outputlevel >= 2 + @printf(" Region %s: energy = %.12f\n", region, eigenvalue(E)) + end + return local_tensor +end + +function eigsolve( + H, init_state; nsweeps, nsites=2, outputlevel=0, updater_kwargs=(;), inserter_kwargs=(;), kws... +) + init_prob = EigsolveProblem(; state=copy(init_state), operator=itn.ProjTTN(H)) + common_sweep_kwargs = (; nsites, outputlevel, updater_kwargs, inserter_kwargs) + kwargs_array = [(; common_sweep_kwargs..., sweep=s) for s in 1:nsweeps] + sweep_iter = sweep_iterator(init_prob, kwargs_array) + converged_prob = alternating_update(sweep_iter; outputlevel, kws...) + return eigenvalue(converged_prob), state(converged_prob) +end + +dmrg(args...; kws...) = eigsolve(args...; kws...) diff --git a/src/solvers/iterators.jl b/src/solvers/iterators.jl new file mode 100644 index 00000000..746de67d --- /dev/null +++ b/src/solvers/iterators.jl @@ -0,0 +1,70 @@ + +# +# sweep_iterator +# + +function sweep_iterator(problem, sweep_kwargs_array) + return [region_iterator(problem; sweep_kwargs...) for sweep_kwargs in sweep_kwargs_array] +end + +function sweep_iterator(problem, nsweeps::Integer) + return sweep_iterator(problem, Iterators.repeated((;), nsweeps)) +end + +# +# step_iterator +# + +step_iterator(args...; kws...) = Iterators.flatten(sweep_iterator(args...; kws...)) + +# +# RegionIterator +# + +@kwdef mutable struct RegionIterator{Problem,RegionPlan} + problem::Problem + region_plan::RegionPlan + which_region_plan::Int = 1 + prev_region = nothing + #extra_kwargs::NamedTuple = (;) +end + +problem(R::RegionIterator) = R.problem +current_region_plan(R::RegionIterator) = R.region_plan[R.which_region_plan] +current_region(R::RegionIterator) = R.region_plan[R.which_region_plan][1] +region_kwargs(R::RegionIterator) = R.region_plan[R.which_region_plan][2] + +function Base.iterate(R::RegionIterator, which=1) + R.which_region_plan = which + region_plan_state = iterate(R.region_plan, which) + isnothing(region_plan_state) && return nothing + (current_region, region_kwargs), next = region_plan_state + + region_iterator_action!(problem(R); region=current_region, prev_region=R.prev_region, region_kwargs...) + R.prev_region = current_region + return R, next +end + +# +# Functions associated with RegionIterator +# + +function region_iterator(problem; nsites=1, sweep_kwargs...) + return RegionIterator(; + problem, region_plan=region_plan(problem; nsites, sweep_kwargs...) + ) +end + +function region_iterator_action!( + problem; region, prev_region=nothing, extracter_kwargs=(;), updater_kwargs=(;), inserter_kwargs=(;), kwargs... +) + local_tensor = extracter!(problem, region; extracter_kwargs..., kwargs...) + local_tensor = prepare_subspace!(problem, local_tensor, region; prev_region, extracter_kwargs..., kwargs...) + local_tensor = updater!(problem, local_tensor, region; updater_kwargs..., kwargs...) + inserter!(problem, local_tensor, region; inserter_kwargs..., kwargs...) + return +end + +function region_plan(problem; nsites, sweep_kwargs...) + return basic_region_plan(state(problem); nsites, sweep_kwargs...) +end diff --git a/src/solvers/region_plans.jl b/src/solvers/region_plans.jl new file mode 100644 index 00000000..ac4ab733 --- /dev/null +++ b/src/solvers/region_plans.jl @@ -0,0 +1,71 @@ +import Graphs: AbstractGraph, AbstractEdge, edges, dst, src, vertices +import NamedGraphs: GraphsExtensions + +#function basic_path_regions(g::AbstractGraph; sweep_kwargs...) +# fwd_sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges(g)] +# return [fwd_sweep..., reverse(fwd_sweep)...] +#end + +function tdvp_regions(g::AbstractGraph, time_step; nsites=1, updater_kwargs, sweep_kwargs...) + @assert nsites==1 + fwd_up_args = (; time=(time_step / 2), updater_kwargs...) + rev_up_args = (; time=(-time_step / 2), updater_kwargs...) + + fwd_sweep = [] + for e in edges(g) + push!(fwd_sweep, ([src(e)], (; updater_kwargs=fwd_up_args, sweep_kwargs...))) + push!(fwd_sweep, (e, (; updater_kwargs=rev_up_args, sweep_kwargs...))) + end + push!(fwd_sweep, ([dst(last(edges(g)))], (; updater_kwargs=fwd_up_args, sweep_kwargs...))) + + # Reverse regions as well as ordering of regions: + rev_sweep = [(reverse(rk[1]), rk[2]) for rk in reverse(fwd_sweep)] + + return [fwd_sweep..., rev_sweep...] +end + +function overlap(ea::AbstractEdge, eb::AbstractEdge) + return intersect([src(ea),dst(ea)], [src(eb),dst(eb)]) +end + +function forward_region(edges, which_edge; nsites=1, region_kwargs=(;)) + current_edge = edges[which_edge] + if nsites == 1 + #handle edge case + if current_edge == last(edges) + overlapping_vertex = only( + union([overlap(e, current_edge) for e in edges[1:(which_edge - 1)]]...) + ) + nonoverlapping_vertex = only( + setdiff([src(current_edge), dst(current_edge)], [overlapping_vertex]) + ) + return [ + ([overlapping_vertex], region_kwargs), ([nonoverlapping_vertex], region_kwargs) + ] + else + future_edges = edges[(which_edge + 1):end] + future_edges = isa(future_edges, AbstractEdge) ? [future_edges] : future_edges + overlapping_vertex = only(union([overlap(e, current_edge) for e in future_edges]...)) + nonoverlapping_vertex = only( + setdiff([src(current_edge), dst(current_edge)], [overlapping_vertex]) + ) + return [([nonoverlapping_vertex], region_kwargs)] + end + elseif nsites == 2 + return [([src(current_edge), dst(current_edge)], region_kwargs)] + end +end + +function basic_region_plan( + graph::AbstractGraph; + nsites, + root_vertex=GraphsExtensions.default_root_vertex(graph), + sweep_kwargs..., +) + edges = GraphsExtensions.post_order_dfs_edges(graph, root_vertex) + fwd_sweep = [ + forward_region(edges, i; nsites, region_kwargs=sweep_kwargs) for i in 1:length(edges) + ] + fwd_sweep = collect(Iterators.flatten(fwd_sweep)) + return [fwd_sweep..., reverse(fwd_sweep)...] +end diff --git a/src/solvers/subspace.jl b/src/solvers/subspace.jl new file mode 100644 index 00000000..43eb80bb --- /dev/null +++ b/src/solvers/subspace.jl @@ -0,0 +1,52 @@ +using ITensors: commonind, dag, dim, directsum, dot, hascommoninds, Index, norm, onehot, uniqueinds, random_itensor + +# TODO: hoist num_expand default value out to a function or similar +function subspace_expand!(problem::EigsolveProblem, local_tensor, region; prev_region, num_expand=4, kws...) + + if isnothing(prev_region) || isa(region, AbstractEdge) + return local_tensor + end + + prev_vertex_set = setdiff(prev_region, region) + (length(prev_vertex_set) != 1) && return local_tensor + prev_vertex = only(prev_vertex_set) + + psi = state(problem) + A = psi[prev_vertex] + + next_vertex = only(filter(v -> (it.hascommoninds(psi[v], A)), region)) + C = psi[next_vertex] + + # Analyze indices of A + # TODO: if "a" is missing, could supply a 1-dim index and put on both A and C? + a = commonind(A, C) + isnothing(a) && return local_tensor + basis_inds = uniqueinds(A, C) + + # Determine maximum value of num_expand + dim_basis = prod(dim.(basis_inds)) + num_expand = min(num_expand, dim_basis - dim(a)) + (num_expand <= 0) && return local_tensor + + # Build new subspace + function linear_map(w) + return w = w - A * (dag(A) * w) + end + random_vector() = random_itensor(basis_inds...) + Q = range_finder(linear_map, random_vector; max_rank=num_expand, oversample=0) + + # Direct sum new space with A to make Ax + qinds = [Index(1, "q$j") for j in 1:num_expand] + Q = [Q[j] * onehot(qinds[j] => 1) => qinds[j] for j in 1:num_expand] + Ax, sa = directsum(A => a, Q...) + + expander = dag(Ax) * A + psi[prev_vertex] = Ax + psi[next_vertex] = expander * C + + # TODO: avoid computing local tensor twice + # while also handling AbstractEdge region case + local_tensor = prod(psi[v] for v in region) + + return local_tensor +end diff --git a/src/solvers/tdvp.jl b/src/solvers/tdvp.jl new file mode 100644 index 00000000..db5a36e8 --- /dev/null +++ b/src/solvers/tdvp.jl @@ -0,0 +1,36 @@ +import ITensorNetworks as itn +using Printf + +@kwdef mutable struct TDVPProblem{State} + state::State + operator + #current_time::Number = 0.0 +end + +state(tdvp::TDVPProblem) = tdvp.state +operator(tdvp::TDVPProblem) = tdvp.operator + +function set(tdvp::TDVPProblem; state=state(tdvp), operator=operator(tdvp)) + return TDVPProblem(; state, operator) +end + +function region_plan(tdvp::TDVPProblem; nsites, time_step, sweep_kwargs...) + return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) +end + +function updater!(tdvp::TDVPProblem, local_tensor, region; outputlevel, kws...) + local_tensor, info = exponentiate_updater(operator(tdvp), local_tensor; kws...) + return local_tensor +end + +function applyexp( + H, init_state, time_points; updater_kwargs=(;), inserter_kwargs=(;), outputlevel=0, kws... +) + init_prob = TDVPProblem(; state=copy(init_state), operator=itn.ProjTTN(H)) + time_steps = diff([0.0, time_points...]) + common_sweep_kwargs = (; outputlevel, updater_kwargs, inserter_kwargs) + kwargs_array = [(; common_sweep_kwargs..., time_step=t) for t in time_steps] + sweep_iter = sweep_iterator(init_prob, kwargs_array) + converged_prob = alternating_update(sweep_iter; outputlevel, kws...) + return state(converged_prob) +end From 3edcc3d7aef409b7b7c3d79df8815288c6b86e30 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Jul 2025 16:16:41 -0400 Subject: [PATCH 03/28] Update solvers codes to latest versions --- src/solvers/adapters.jl | 44 ---- src/solvers/alternating_update.jl | 23 -- src/solvers/applyexp.jl | 103 ++++++++ src/solvers/defaults.jl | 66 ------ src/solvers/eigsolve.jl | 68 ++++-- src/solvers/extracter.jl | 17 ++ src/solvers/fitting.jl | 112 +++++++++ src/solvers/inserter.jl | 33 +++ src/solvers/iterators.jl | 100 +++++--- src/solvers/local_solvers/contract.jl | 13 -- src/solvers/local_solvers/dmrg_x.jl | 22 -- src/solvers/local_solvers/eigsolve.jl | 29 +++ src/solvers/local_solvers/exponentiate.jl | 29 +++ src/solvers/local_solvers/linsolve.jl | 26 --- src/solvers/local_solvers/runge_kutta.jl | 25 ++ src/solvers/operator_map.jl | 42 ++++ src/solvers/permute_indices.jl | 18 ++ src/solvers/region_plans/dfs_plans.jl | 22 ++ src/solvers/region_plans/euler_plans.jl | 13 ++ src/solvers/region_plans/euler_tour.jl | 48 ++++ src/solvers/region_plans/tdvp_region_plans.jl | 44 ++++ src/solvers/subspace/densitymatrix.jl | 74 ++++++ src/solvers/subspace/ortho_subspace.jl | 77 ++++++ src/solvers/subspace/subspace.jl | 48 ++++ src/solvers/sweep_plans/sweep_plans.jl | 221 ------------------ src/solvers/sweep_solve.jl | 40 ++++ src/solvers/truncation_parameters.jl | 14 ++ 27 files changed, 907 insertions(+), 464 deletions(-) delete mode 100644 src/solvers/alternating_update.jl create mode 100644 src/solvers/applyexp.jl delete mode 100644 src/solvers/defaults.jl create mode 100644 src/solvers/extracter.jl create mode 100644 src/solvers/fitting.jl create mode 100644 src/solvers/inserter.jl delete mode 100644 src/solvers/local_solvers/contract.jl delete mode 100644 src/solvers/local_solvers/dmrg_x.jl create mode 100644 src/solvers/local_solvers/eigsolve.jl create mode 100644 src/solvers/local_solvers/exponentiate.jl delete mode 100644 src/solvers/local_solvers/linsolve.jl create mode 100644 src/solvers/local_solvers/runge_kutta.jl create mode 100644 src/solvers/operator_map.jl create mode 100644 src/solvers/permute_indices.jl create mode 100644 src/solvers/region_plans/dfs_plans.jl create mode 100644 src/solvers/region_plans/euler_plans.jl create mode 100644 src/solvers/region_plans/euler_tour.jl create mode 100644 src/solvers/region_plans/tdvp_region_plans.jl create mode 100644 src/solvers/subspace/densitymatrix.jl create mode 100644 src/solvers/subspace/ortho_subspace.jl create mode 100644 src/solvers/subspace/subspace.jl delete mode 100644 src/solvers/sweep_plans/sweep_plans.jl create mode 100644 src/solvers/sweep_solve.jl create mode 100644 src/solvers/truncation_parameters.jl diff --git a/src/solvers/adapters.jl b/src/solvers/adapters.jl index 7278d245..7c033d8e 100644 --- a/src/solvers/adapters.jl +++ b/src/solvers/adapters.jl @@ -1,48 +1,4 @@ -# -# RepeatIterator -# -# An "iterator of iterators", plugging in initialization -# arguments into an `init` function each repetition -# - -mutable struct RepeatIterator{Iterator} - iterator::Iterator - keyword_args::Vector{<:NamedTuple} -end - -function RepeatIterator(iter, n::Integer, kwargs::NamedTuple) - return RepeatIterator(iter, fill(kwargs, n)) -end - -# -# Version 1: output the initialized iterator each time -# -function Base.iterate(R::RepeatIterator) - inner_next = iterate(R.iterator) - isnothing(inner_next) && return nothing - (item, inner_state) = inner_next - return item, (1, inner_state) -end -function Base.iterate(R::RepeatIterator, state) - state = Base.iterate(keyword_args) - isnothing(state) && return nothing - (kwargs, which) = state - R.iterator = init(R.iterator; kwargs...) - return R, which -end - -## -## Version 2: output the initialized iterator each time -## -#function Base.iterate(R::RepeatIterator, state=1) -# state = Base.iterate(keyword_args) -# isnothing(state) && return nothing -# (kwargs, which) = state -# R.iterator = init(R.iterator; kwargs...) -# return R, which -#end - # # TupleRegionIterator # diff --git a/src/solvers/alternating_update.jl b/src/solvers/alternating_update.jl deleted file mode 100644 index 1e2b1013..00000000 --- a/src/solvers/alternating_update.jl +++ /dev/null @@ -1,23 +0,0 @@ - -function sweep_callback(problem; outputlevel, sweep, kws...) - if outputlevel >= 1 - println("Done with sweep $sweep") - #print(" cpu_time=", round(sweep_time; digits=3)) - #println() - flush(stdout) - end - return problem -end - -function alternating_update(sweep_iterator; outputlevel=0, kwargs...) - for (sweep, region_iter) in enumerate(sweep_iterator) - prob = problem(region_iter) - for (region, region_kwargs) in region_tuples(region_iter) - end - prob = sweep_callback( - prob; nsweeps=length(sweep_iterator), outputlevel, sweep, kwargs... - ) - #TODO need to update `problem` reference held within region iterator(s) - end - return problem(last(sweep_iterator)) -end diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl new file mode 100644 index 00000000..8832cfcb --- /dev/null +++ b/src/solvers/applyexp.jl @@ -0,0 +1,103 @@ +import ITensorNetworks as itn +using Printf + +@kwdef mutable struct ApplyExpProblem{State} + state::State + operator + current_time::Number = 0.0 +end + +state(tdvp::ApplyExpProblem) = tdvp.state +operator(tdvp::ApplyExpProblem) = tdvp.operator +current_time(tdvp::ApplyExpProblem) = tdvp.current_time + +function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...) + return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) +end + +function updater( + T::ApplyExpProblem, + local_state, + region_iterator; + nsites, + time_step, + solver=runge_kutta_solver, + outputlevel, + kws..., +) + local_state, info = solver(x->optimal_map(operator(T), x), time_step, local_state; kws...) + + if nsites==1 + curr_reg = current_region(region_iterator) + next_reg = next_region(region_iterator) + if !isnothing(next_reg) && next_reg != curr_reg + next_edge = first(itn.edge_sequence_between_regions(state(T), curr_reg, next_reg)) + v1, v2 = itn.src(next_edge), itn.dst(next_edge) + psi = copy(state(T)) + psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2])) + shifted_operator = itn.position(operator(T), psi, itn.NamedEdge(v1=>v2)) + R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...) + local_state = psi[v1]*R_t + end + end + + curr_time = current_time(T) + time_step + T = setproperties(T; current_time=curr_time) + + return T, local_state +end + +function applyexp_sweep_printer( + region_iterator; outputlevel, sweep, nsweeps, process_time=identity, kws... +) + if outputlevel >= 1 + T = problem(region_iterator) + @printf(" Current time = %s, ", process_time(current_time(T))) + @printf("maxlinkdim=%d", itn.maxlinkdim(state(T))) + println() + flush(stdout) + end +end + +function applyexp( + init_prob, + exponents; + extracter_kwargs=(;), + updater_kwargs=(;), + inserter_kwargs=(;), + outputlevel=0, + nsites=1, + tdvp_order=4, + sweep_printer=applyexp_sweep_printer, + kws..., +) + time_steps = diff([0.0, exponents...])[2:end] + sweep_kws = (; + outputlevel, extracter_kwargs, inserter_kwargs, nsites, tdvp_order, updater_kwargs + ) + kws_array = [(; sweep_kws..., time_step=t) for t in time_steps] + sweep_iter = sweep_iterator(init_prob, kws_array) + converged_prob = sweep_solve(sweep_iter; outputlevel, sweep_printer, kws...) + return state(converged_prob) +end + +function applyexp(H, init_state, exponents; kws...) + init_prob = ApplyExpProblem(; + state=permute_indices(init_state), operator=itn.ProjTTN(permute_indices(H)) + ) + return applyexp(init_prob, exponents; kws...) +end + +process_real_times(z) = round(-imag(z); digits=10) + +function tdvp( + H, + init_state, + time_points; + process_time=process_real_times, + sweep_printer=(a...; k...)->applyexp_sweep_printer(a...; process_time, k...), + kws..., +) + exponents = [-im*t for t in time_points] + return applyexp(H, init_state, exponents; sweep_printer, kws...) +end diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl deleted file mode 100644 index 5f25e087..00000000 --- a/src/solvers/defaults.jl +++ /dev/null @@ -1,66 +0,0 @@ -using Printf: @printf, @sprintf -default_outputlevel() = 0 -default_nsites() = 2 -default_nsweeps() = 1 #? or nothing? -default_extracter() = default_extracter -default_inserter() = default_inserter -default_checkdone() = (; kws...) -> false -default_transform_operator() = nothing - -format(x) = @sprintf("%s", x) -format(x::AbstractFloat) = @sprintf("%.1E", x) - -function default_region_printer(; - inserter_kwargs, - outputlevel, - state, - sweep_plan, - spec, - which_region_update, - which_sweep, - kwargs..., -) - if outputlevel >= 2 - region = first(sweep_plan[which_region_update]) - @printf("Sweep %d, region=%s \n", which_sweep, region) - print(" Truncated using") - for key in [:cutoff, :maxdim, :mindim] - if haskey(inserter_kwargs, key) - print(" ", key, "=", format(inserter_kwargs[key])) - end - end - println() - if spec != nothing - @printf( - " Trunc. err=%.2E, bond dimension %d\n", - spec.truncerr, - linkdim(state, edgetype(state)(region...)) - ) - end - flush(stdout) - end -end - -#ToDo: Implement sweep_time_printer more generally -#ToDo: Implement more printers -#ToDo: Move to another file? -function default_sweep_time_printer(; outputlevel, which_sweep, kwargs...) - if outputlevel >= 1 - sweeps_per_step = order ÷ 2 - if which_sweep % sweeps_per_step == 0 - current_time = (which_sweep / sweeps_per_step) * time_step - println("Current time (sweep $which_sweep) = ", round(current_time; digits=3)) - end - end - return nothing -end - -function default_sweep_printer(; outputlevel, state, which_sweep, sweep_time, kwargs...) - if outputlevel >= 1 - print("After sweep ", which_sweep, ":") - print(" maxlinkdim=", maxlinkdim(state)) - print(" cpu_time=", round(sweep_time; digits=3)) - println() - flush(stdout) - end -end diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 28e6903d..d2b0b042 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,5 +1,4 @@ import ITensorNetworks as itn -import ITensors as it using Printf @kwdef mutable struct EigsolveProblem{State,Operator} @@ -12,29 +11,66 @@ eigenvalue(E::EigsolveProblem) = E.eigenvalue state(E::EigsolveProblem) = E.state operator(E::EigsolveProblem) = E.operator -function set( - E::EigsolveProblem; state=state(E), operator=operator(E), eigenvalue=eigenvalue(E) +function updater( + E::EigsolveProblem, + local_state, + region_iterator; + outputlevel, + solver=eigsolve_solver, + kws..., ) - return EigsolveProblem(; state, operator, eigenvalue) + eigval, local_state = solver(ψ->optimal_map(operator(E), ψ), local_state; kws...) + E = setproperties(E; eigenvalue=eigval) + if outputlevel >= 2 + @printf(" Region %s: energy = %.12f\n", current_region(region_iterator), eigenvalue(E)) + end + return E, local_state end -function updater!(E::EigsolveProblem, local_tensor, region; outputlevel, kws...) - E.eigenvalue, local_tensor = eigsolve_updater(operator(E), local_tensor; kws...) - if outputlevel >= 2 - @printf(" Region %s: energy = %.12f\n", region, eigenvalue(E)) +function eigsolve_sweep_printer(region_iterator; outputlevel, sweep, nsweeps, kws...) + if outputlevel >= 1 + if nsweeps >= 10 + @printf("After sweep %02d/%d ", sweep, nsweeps) + else + @printf("After sweep %d/%d ", sweep, nsweeps) + end + E = problem(region_iterator) + @printf("eigenvalue=%.12f ", eigenvalue(E)) + @printf("maxlinkdim=%d", itn.maxlinkdim(state(E))) + println() + flush(stdout) end - return local_tensor end function eigsolve( - H, init_state; nsweeps, nsites=2, outputlevel=0, updater_kwargs=(;), inserter_kwargs=(;), kws... + init_prob; + nsweeps, + nsites=1, + outputlevel=0, + extracter_kwargs=(;), + updater_kwargs=(;), + inserter_kwargs=(;), + sweep_printer=eigsolve_sweep_printer, + kws..., ) - init_prob = EigsolveProblem(; state=copy(init_state), operator=itn.ProjTTN(H)) - common_sweep_kwargs = (; nsites, outputlevel, updater_kwargs, inserter_kwargs) - kwargs_array = [(; common_sweep_kwargs..., sweep=s) for s in 1:nsweeps] - sweep_iter = sweep_iterator(init_prob, kwargs_array) - converged_prob = alternating_update(sweep_iter; outputlevel, kws...) - return eigenvalue(converged_prob), state(converged_prob) + sweep_iter = sweep_iterator( + init_prob, + nsweeps; + nsites, + outputlevel, + extracter_kwargs, + updater_kwargs, + inserter_kwargs, + ) + prob = sweep_solve(sweep_iter; outputlevel, sweep_printer, kws...) + return eigenvalue(prob), state(prob) +end + +function eigsolve(H, init_state; kws...) + init_prob = EigsolveProblem(; + state=permute_indices(init_state), operator=itn.ProjTTN(permute_indices(H)) + ) + return eigsolve(init_prob; kws...) end dmrg(args...; kws...) = eigsolve(args...; kws...) diff --git a/src/solvers/extracter.jl b/src/solvers/extracter.jl new file mode 100644 index 00000000..50a396f6 --- /dev/null +++ b/src/solvers/extracter.jl @@ -0,0 +1,17 @@ +import ConstructionBase: setproperties + +function extracter(problem, region_iterator; sweep, trunc=(;), kws...) + trunc = truncation_parameters(sweep; trunc...) + region = current_region(region_iterator) + psi = itn.orthogonalize(state(problem), region) + local_state = prod(psi[v] for v in region) + problem = setproperties(problem; state=psi) + + problem, local_state = subspace_expand( + problem, local_state, region_iterator; sweep, trunc, kws... + ) + + shifted_operator = itn.position(operator(problem), state(problem), region) + + return setproperties(problem; operator=shifted_operator), local_state +end diff --git a/src/solvers/fitting.jl b/src/solvers/fitting.jl new file mode 100644 index 00000000..0995124b --- /dev/null +++ b/src/solvers/fitting.jl @@ -0,0 +1,112 @@ +import ITensorNetworks as itn +import ITensors as it +import Graphs: vertices +import NamedGraphs.PartitionedGraphs as npg +using NamedGraphs: AbstractNamedGraph, NamedEdge +using Printf + +@kwdef mutable struct FittingProblem{State<:itn.AbstractBeliefPropagationCache} + state::State + ket_graph::AbstractNamedGraph + overlap::Number = 0 + gauge_region +end + +overlap(F::FittingProblem) = F.overlap +state(F::FittingProblem) = F.state +ket_graph(F::FittingProblem) = F.ket_graph +gauge_region(F::FittingProblem) = F.gauge_region + +function ket(F::FittingProblem) + ket_vertices = vertices(ket_graph(F)) + return first(itn.induced_subgraph(itn.tensornetwork(state(F)), ket_vertices)) +end + +function extracter(problem::FittingProblem, region_iterator; sweep, kws...) + region = current_region(region_iterator) + prev_region = gauge_region(problem) + tn = state(problem) + path = itn.edge_sequence_between_regions(ket_graph(problem), prev_region, region) + tn = itn.gauge_walk(itn.Algorithm("orthogonalize"), tn, path) + pe_path = npg.partitionedges(itn.partitioned_tensornetwork(tn), path) + tn = itn.update( + itn.Algorithm("bp"), tn, pe_path; message_update_function_kwargs=(; normalize=false) + ) + local_tensor = itn.environment(tn, region) + sequence = itn.contraction_sequence(local_tensor; alg="optimal") + local_tensor = dag(it.contract(local_tensor; sequence)) + #problem, local_tensor = subspace_expand(problem, local_tensor, region; sweep, kws...) + return setproperties(problem; state=tn, gauge_region=region), local_tensor +end + +function updater(F::FittingProblem, local_tensor, region; outputlevel, kws...) + n = (local_tensor * dag(local_tensor))[] + F = setproperties(F; overlap=n / sqrt(n)) + if outputlevel >= 2 + @printf(" Region %s: squared overlap = %.12f\n", region, overlap(F)) + end + return F, local_tensor +end + +function region_plan(F::FittingProblem; nsites, sweep_kwargs...) + return euler_sweep(ket_graph(F); nsites, sweep_kwargs...) +end + +function fit_tensornetwork( + overlap_network, + args...; + nsweeps=25, + nsites=1, + outputlevel=0, + extracter_kwargs=(;), + updater_kwargs=(;), + inserter_kwargs=(;), + normalize=true, + kws..., +) + bpc = itn.BeliefPropagationCache(overlap_network, args...) + ket_vertices = itn.ket_vertices(overlap_network) + ket_graph = first( + itn.induced_subgraph( + itn.underlying_graph(overlap_network), itn.ket_vertices(overlap_network) + ), + ) + init_prob = FittingProblem(; + ket_graph=ket_graph, state=bpc, gauge_region=collect(vertices(ket_graph)) + ) + + inserter_kwargs = (; inserter_kwargs..., normalize, set_orthogonal_region=false) + common_sweep_kwargs = (; nsites, outputlevel, updater_kwargs, inserter_kwargs) + kwargs_array = [(; common_sweep_kwargs..., sweep=s) for s in 1:nsweeps] + sweep_iter = sweep_iterator(init_prob, kwargs_array) + converged_prob = sweep_solve(sweep_iter; outputlevel, kws...) + return itn.rename_vertices(itn.inv_vertex_map(overlap_network), ket(converged_prob)) +end + +function fit_tensornetwork(tn, init_state, args...; kwargs...) + return fit_tensornetwork(itn.inner_network(tn, init_state), args; kwargs...) +end + +function itn.truncate(tn; maxdim=default_maxdim(), cutoff=default_cutoff(), kwargs...) + init_state = itn.ITensorNetwork( + v -> inds -> it.delta(inds), itn.siteinds(tn); link_space=maxdim + ) + overlap_network = itn.inner_network(tn, init_state) + inserter_kwargs = (; trunc=(; cutoff, maxdim)) + return fit_tensornetwork(overlap_network; inserter_kwargs, kwargs...) +end + +function itn.apply( + A::itn.ITensorNetwork, + x::itn.ITensorNetwork; + maxdim=default_maxdim(), + cutoff=default_cutoff(), + kwargs..., +) + init_state = itn.ITensorNetwork( + v -> inds -> it.delta(inds), itn.siteinds(x); link_space=maxdim + ) + overlap_network = itn.inner_network(x, A, init_state) + inserter_kwargs = (; trunc=(; cutoff, maxdim)) + return fit_tensornetwork(overlap_network; inserter_kwargs, kwargs...) +end diff --git a/src/solvers/inserter.jl b/src/solvers/inserter.jl new file mode 100644 index 00000000..64a03921 --- /dev/null +++ b/src/solvers/inserter.jl @@ -0,0 +1,33 @@ +import ConstructionBase: setproperties + +function inserter( + problem, + local_tensor, + region_iterator; + normalize=false, + set_orthogonal_region=true, + sweep, + trunc=(;), + kws..., +) + trunc = truncation_parameters(sweep; trunc...) + + region = current_region(region_iterator) + psi = copy(state(problem)) + if length(region) == 1 + C = local_tensor + elseif length(region) == 2 + e = ng.edgetype(psi)(first(region), last(region)) + indsTe = it.inds(psi[first(region)]) + tags = it.tags(psi, e) + U, C, _ = it.factorize(local_tensor, indsTe; tags, trunc...) + itn.@preserve_graph psi[first(region)] = U + else + error("Region of length $(length(region)) not currently supported") + end + v = last(region) + itn.@preserve_graph psi[v] = C + psi = set_orthogonal_region ? itn.set_ortho_region(psi, [v]) : psi + normalize && itn.@preserve_graph psi[v] = psi[v] / norm(psi[v]) + return setproperties(problem; state=psi) +end diff --git a/src/solvers/iterators.jl b/src/solvers/iterators.jl index 746de67d..aaf530cf 100644 --- a/src/solvers/iterators.jl +++ b/src/solvers/iterators.jl @@ -1,21 +1,43 @@ - # -# sweep_iterator +# SweepIterator # -function sweep_iterator(problem, sweep_kwargs_array) - return [region_iterator(problem; sweep_kwargs...) for sweep_kwargs in sweep_kwargs_array] +mutable struct SweepIterator + sweep_kws + region_iter + which_sweep::Int end -function sweep_iterator(problem, nsweeps::Integer) - return sweep_iterator(problem, Iterators.repeated((;), nsweeps)) +problem(S::SweepIterator) = problem(S.region_iter) + +Base.length(S::SweepIterator) = length(S.sweep_kws) + +function Base.iterate(S::SweepIterator, which=nothing) + if isnothing(which) + sweep_kws_state = iterate(S.sweep_kws) + else + sweep_kws_state = iterate(S.sweep_kws, which) + end + isnothing(sweep_kws_state) && return nothing + current_sweep_kws, next = sweep_kws_state + + if !isnothing(which) + S.region_iter = region_iterator( + problem(S.region_iter); sweep=S.which_sweep, current_sweep_kws... + ) + end + S.which_sweep += 1 + return S.region_iter, next end -# -# step_iterator -# +function sweep_iterator(problem, sweep_kws) + region_iter = region_iterator(problem; sweep=1, first(sweep_kws)...) + return SweepIterator(sweep_kws, region_iter, 1) +end -step_iterator(args...; kws...) = Iterators.flatten(sweep_iterator(args...; kws...)) +function sweep_iterator(problem, nsweeps::Integer; sweep_kws...) + return sweep_iterator(problem, Iterators.repeated(sweep_kws, nsweeps)) +end # # RegionIterator @@ -24,24 +46,27 @@ step_iterator(args...; kws...) = Iterators.flatten(sweep_iterator(args...; kws.. @kwdef mutable struct RegionIterator{Problem,RegionPlan} problem::Problem region_plan::RegionPlan - which_region_plan::Int = 1 - prev_region = nothing - #extra_kwargs::NamedTuple = (;) + which_region::Int = 1 end problem(R::RegionIterator) = R.problem -current_region_plan(R::RegionIterator) = R.region_plan[R.which_region_plan] -current_region(R::RegionIterator) = R.region_plan[R.which_region_plan][1] -region_kwargs(R::RegionIterator) = R.region_plan[R.which_region_plan][2] +current_region_plan(R::RegionIterator) = R.region_plan[R.which_region] +current_region(R::RegionIterator) = current_region_plan(R)[1] +region_kwargs(R::RegionIterator) = current_region_plan(R)[2] +function previous_region(R::RegionIterator) + R.which_region==1 ? nothing : R.region_plan[R.which_region - 1][1] +end +function next_region(R::RegionIterator) + R.which_region==length(R.region_plan) ? nothing : R.region_plan[R.which_region + 1][1] +end +is_last_region(R::RegionIterator) = isnothing(next_region(R)) function Base.iterate(R::RegionIterator, which=1) - R.which_region_plan = which + R.which_region = which region_plan_state = iterate(R.region_plan, which) isnothing(region_plan_state) && return nothing (current_region, region_kwargs), next = region_plan_state - - region_iterator_action!(problem(R); region=current_region, prev_region=R.prev_region, region_kwargs...) - R.prev_region = current_region + R.problem = region_iterator_action(problem(R), R; region_kwargs...) return R, next end @@ -49,22 +74,31 @@ end # Functions associated with RegionIterator # -function region_iterator(problem; nsites=1, sweep_kwargs...) - return RegionIterator(; - problem, region_plan=region_plan(problem; nsites, sweep_kwargs...) - ) +function region_iterator(problem; sweep_kwargs...) + return RegionIterator(; problem, region_plan=region_plan(problem; sweep_kwargs...)) end -function region_iterator_action!( - problem; region, prev_region=nothing, extracter_kwargs=(;), updater_kwargs=(;), inserter_kwargs=(;), kwargs... +function region_iterator_action( + problem, + region_iterator; + extracter_kwargs=(;), + updater_kwargs=(;), + inserter_kwargs=(;), + sweep, + kws..., ) - local_tensor = extracter!(problem, region; extracter_kwargs..., kwargs...) - local_tensor = prepare_subspace!(problem, local_tensor, region; prev_region, extracter_kwargs..., kwargs...) - local_tensor = updater!(problem, local_tensor, region; updater_kwargs..., kwargs...) - inserter!(problem, local_tensor, region; inserter_kwargs..., kwargs...) - return + problem, local_state = extracter( + problem, region_iterator; extracter_kwargs..., sweep, kws... + ) + problem, local_state = updater( + problem, local_state, region_iterator; updater_kwargs..., kws... + ) + problem = inserter( + problem, local_state, region_iterator; sweep, inserter_kwargs..., kws... + ) + return problem end -function region_plan(problem; nsites, sweep_kwargs...) - return basic_region_plan(state(problem); nsites, sweep_kwargs...) +function region_plan(problem; kws...) + return euler_sweep(state(problem); kws...) end diff --git a/src/solvers/local_solvers/contract.jl b/src/solvers/local_solvers/contract.jl deleted file mode 100644 index bffefdef..00000000 --- a/src/solvers/local_solvers/contract.jl +++ /dev/null @@ -1,13 +0,0 @@ -function contract_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - internal_kwargs, -) - P = projected_operator![] - return contract_ket(P, ITensor(one(Bool))), (;) -end diff --git a/src/solvers/local_solvers/dmrg_x.jl b/src/solvers/local_solvers/dmrg_x.jl deleted file mode 100644 index 3c3ae429..00000000 --- a/src/solvers/local_solvers/dmrg_x.jl +++ /dev/null @@ -1,22 +0,0 @@ -using ITensors: ITensor, contract, dag, onehot, uniqueind -using ITensors.NDTensors: array -using LinearAlgebra: eigen - -function dmrg_x_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - internal_kwargs, -) - H = contract(projected_operator![], ITensor(true)) - D, U = eigen(H; ishermitian=true) - u = uniqueind(U, H) - max_overlap, max_ind = findmax(abs, array(dag(init) * U)) - U_max = U * dag(onehot(u => max_ind)) - eigvals = [((onehot(u => max_ind)' * D) * dag(onehot(u => max_ind)))[]] - return U_max, (; eigvals) -end diff --git a/src/solvers/local_solvers/eigsolve.jl b/src/solvers/local_solvers/eigsolve.jl new file mode 100644 index 00000000..96f1c9fb --- /dev/null +++ b/src/solvers/local_solvers/eigsolve.jl @@ -0,0 +1,29 @@ +using KrylovKit: KrylovKit + +function eigsolve_solver( + operator, + init, + howmany=1; + which_eigval=:SR, + ishermitian=true, + tol=1E-14, + krylovdim=3, + maxiter=1, + verbosity=0, + eager=false, + kws..., +) + vals, vecs, info = KrylovKit.eigsolve( + operator, + init, + howmany, + which_eigval; + ishermitian, + tol, + krylovdim, + maxiter, + verbosity, + eager, + ) + return vals[1], vecs[1] +end diff --git a/src/solvers/local_solvers/exponentiate.jl b/src/solvers/local_solvers/exponentiate.jl new file mode 100644 index 00000000..b0929916 --- /dev/null +++ b/src/solvers/local_solvers/exponentiate.jl @@ -0,0 +1,29 @@ +using KrylovKit: KrylovKit + +function exponentiate_solver( + operator, + time, + init; + krylovdim=30, + maxiter=100, + verbosity=0, + tol=1E-12, + ishermitian=true, + issymmetric=true, + eager=true, + kws..., +) + result, exp_info = KrylovKit.exponentiate( + operator, + time, + init; + eager, + krylovdim, + maxiter, + verbosity, + tol, + ishermitian, + issymmetric, + ) + return result, exp_info +end diff --git a/src/solvers/local_solvers/linsolve.jl b/src/solvers/local_solvers/linsolve.jl deleted file mode 100644 index c5b8c4c6..00000000 --- a/src/solvers/local_solvers/linsolve.jl +++ /dev/null @@ -1,26 +0,0 @@ -using KrylovKit: linsolve - -function linsolve_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - region_kwargs, - ishermitian=false, - tol=1E-14, - krylovdim=30, - maxiter=100, - verbosity=0, - a₀, - a₁, -) - P = projected_operator![] - b = dag(only(proj_mps(P))) - x, info = linsolve( - P, b, init, a₀, a₁; ishermitian=false, tol, krylovdim, maxiter, verbosity - ) - return x, (;) -end diff --git a/src/solvers/local_solvers/runge_kutta.jl b/src/solvers/local_solvers/runge_kutta.jl new file mode 100644 index 00000000..03294f0d --- /dev/null +++ b/src/solvers/local_solvers/runge_kutta.jl @@ -0,0 +1,25 @@ + +function runge_kutta_2(H, t, ψ0) + Hψ = H(ψ0) + H2ψ = H(Hψ) + return (ψ0 + t * Hψ + (t^2 / 2) * H2ψ) +end + +function runge_kutta_4(H, t, ψ0) + k1 = H(ψ0) + k2 = k1 + (t / 2) * H(k1) + k3 = k1 + (t / 2) * H(k2) + k4 = k1 + t * H(k3) + return ψ0 + (t / 6) * (k1 + 2 * k2 + 2 * k3 + k4) +end + +function runge_kutta_solver(H, time, ψ; order=4, kws...) + if order == 4 + Hψ = runge_kutta_4(H, time, ψ) + elseif order == 2 + Hψ = runge_kutta_2(H, time, ψ) + else + error("For runge_kutta_solver, must specify `order` keyword") + end + return Hψ, (;) +end diff --git a/src/solvers/operator_map.jl b/src/solvers/operator_map.jl new file mode 100644 index 00000000..309d85ad --- /dev/null +++ b/src/solvers/operator_map.jl @@ -0,0 +1,42 @@ +import ITensorNetworks as itn + +function optimal_map(P::itn.ProjTTN, ψ) + envs = [itn.environment(P, e) for e in itn.incident_edges(P)] + site_ops = [itn.operator(P)[s] for s in itn.sites(P)] + contract_list = [envs..., site_ops..., ψ] + sequence = itn.contraction_sequence(contract_list; alg="optimal") + Pψ = itn.contract(contract_list; sequence) + return noprime(Pψ) +end + +# This function is a workaround for the slow contraction order +# heuristic in ITensorNetworks/src/treetensornetworks/projttns/projttn.jl +# in the projected_operator_tensors(P::ProjTTN) function (line 97 or so) +function operator_map(P::itn.ProjTTN, ψ) + ψ = copy(ψ) + if itn.on_edge(P) + for edge in itn.incident_edges(P) + ψ *= itn.environment(P, edge) + end + else + region = itn.sites(P) + ie = itn.incident_edges(P) + # TODO: improvement ideas + # - check which vertex (first(region) vs. last(region) + # has more incident edges and contract those environments first + for edge in ie + if itn.dst(edge) == first(region) + ψ *= itn.environment(P, edge) + end + end + for s in itn.sites(P) + ψ *= itn.operator(P)[s] + end + for edge in ie + if itn.dst(edge) != first(region) + ψ *= itn.environment(P, edge) + end + end + end + return noprime(ψ) +end diff --git a/src/solvers/permute_indices.jl b/src/solvers/permute_indices.jl new file mode 100644 index 00000000..927f6b87 --- /dev/null +++ b/src/solvers/permute_indices.jl @@ -0,0 +1,18 @@ +import ITensorNetworks as itn +using ITensors + +function permute_indices(tn) + si = itn.siteinds(tn) + ptn = copy(tn) + for v in itn.vertices(tn) + is = inds(tn[v]) + ls = setdiff(is, si[v]) + isempty(ls) && continue + new_is = [first(ls), si[v]...] + if length(ls) >= 2 + new_is = vcat(new_is, ls[2:end]) + end + ptn[v] = permute(tn[v], new_is) + end + return ptn +end diff --git a/src/solvers/region_plans/dfs_plans.jl b/src/solvers/region_plans/dfs_plans.jl new file mode 100644 index 00000000..facbaca6 --- /dev/null +++ b/src/solvers/region_plans/dfs_plans.jl @@ -0,0 +1,22 @@ +import Graphs: dst, src +import NamedGraphs.GraphsExtensions: + default_root_vertex, post_order_dfs_edges, post_order_dfs_vertices + +function post_order_dfs_plan( + graph; nsites, root_vertex=default_root_vertex(graph), sweep_kwargs... +) + if nsites == 1 + vertices = post_order_dfs_vertices(graph, root_vertex) + fwd_sweep = [([v], sweep_kwargs) for v in vertices] + elseif nsites == 2 + edges = post_order_dfs_edges(graph, root_vertex) + fwd_sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges] + end + return fwd_sweep +end + +function post_order_dfs_sweep(args...; kws...) + fwd_sweep = post_order_dfs_plan(args...; kws...) + rev_sweep = [(reverse(reg_kws[1]), reg_kws[2]) for reg_kws in reverse(fwd_sweep)] + return [fwd_sweep..., rev_sweep...] +end diff --git a/src/solvers/region_plans/euler_plans.jl b/src/solvers/region_plans/euler_plans.jl new file mode 100644 index 00000000..f3b9a96e --- /dev/null +++ b/src/solvers/region_plans/euler_plans.jl @@ -0,0 +1,13 @@ +import Graphs: dst, src +import NamedGraphs.GraphsExtensions: default_root_vertex + +function euler_sweep(graph; nsites, root_vertex=default_root_vertex(graph), sweep_kwargs...) + if nsites == 1 + vertices = euler_tour_vertices(graph, root_vertex) + sweep = [([v], sweep_kwargs) for v in vertices] + elseif nsites == 2 + edges = euler_tour_edges(graph, root_vertex) + sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges] + end + return sweep +end diff --git a/src/solvers/region_plans/euler_tour.jl b/src/solvers/region_plans/euler_tour.jl new file mode 100644 index 00000000..6aeb0029 --- /dev/null +++ b/src/solvers/region_plans/euler_tour.jl @@ -0,0 +1,48 @@ +import Graphs: dst, edges, src, vertices +import NamedGraphs as ng + +function compute_adjacencies(G) + adj = Dict(v => Vector{ng.vertextype(G)}() for v in vertices(G)) + for e in edges(G) + push!(adj[src(e)], dst(e)) + push!(adj[dst(e)], src(e)) + end + return adj +end + +function euler_tour_edges(G, start_vertex) + adj = compute_adjacencies(G) + etype = ng.edgetype(G) + vtype = ng.vertextype(G) + visited = Set{Tuple{vtype,vtype}}() + tour = Vector{etype}() + stack = [start_vertex] + while !isempty(stack) + u = stack[end] + pushed = false + for v in adj[u] + if (u, v) ∉ visited + push!(visited, (u, v)) + push!(visited, (v, u)) + push!(tour, etype(u => v)) + push!(stack, v) + pushed = true + break # handle one neighbor at a time + end + end + if !pushed + pop!(stack) + if !isempty(stack) + v = stack[end] + push!(tour, etype(u => v)) # Backtracking step + end + end + end + return tour +end + +function euler_tour_vertices(G, start_vertex) + edges = euler_tour_edges(G, start_vertex) + isempty(edges) && return Vector{eltype(vertices(G))}[] + return [src(edges[1]), dst.(edges)...] +end diff --git a/src/solvers/region_plans/tdvp_region_plans.jl b/src/solvers/region_plans/tdvp_region_plans.jl new file mode 100644 index 00000000..0624d666 --- /dev/null +++ b/src/solvers/region_plans/tdvp_region_plans.jl @@ -0,0 +1,44 @@ +function tdvp_sub_time_steps(tdvp_order) + if tdvp_order == 1 + return [1.0] + elseif tdvp_order == 2 + return [1 / 2, 1 / 2] + elseif tdvp_order == 4 + s = (2 - 2^(1 / 3))^(-1) + return [s/2, s/2, 1/2 - s, 1/2 - s, s/2, s/2] + else + error("TDVP order of $tdvp_order not supported") + end +end + +function first_order_sweep( + graph, time_step, dir=Base.Forward; updater_kwargs, nsites, kws... +) + basic_fwd_sweep = post_order_dfs_plan(graph; nsites, kws...) + updater_kwargs = (; nsites, time_step, updater_kwargs...) + sweep = [] + for (j, (region, region_kws)) in enumerate(basic_fwd_sweep) + push!(sweep, (region, (; nsites, updater_kwargs, region_kws...))) + if length(region) == 2 && j < length(basic_fwd_sweep) + rev_kwargs = (; updater_kwargs..., time_step=(-updater_kwargs.time_step)) + push!(sweep, ([last(region)], (; updater_kwargs=rev_kwargs, region_kws...))) + end + end + if dir==Base.Reverse + # Reverse regions as well as ordering of regions + sweep = [(reverse(reg_kws[1]), reg_kws[2]) for reg_kws in reverse(sweep)] + end + return sweep +end + +function tdvp_regions(graph, time_step; updater_kwargs, tdvp_order, nsites, kws...) + sweep_plan = [] + for (step, weight) in enumerate(tdvp_sub_time_steps(tdvp_order)) + dir = isodd(step) ? Base.Forward : Base.Reverse + append!( + sweep_plan, + first_order_sweep(graph, weight*time_step, dir; updater_kwargs, nsites, kws...), + ) + end + return sweep_plan +end diff --git a/src/solvers/subspace/densitymatrix.jl b/src/solvers/subspace/densitymatrix.jl new file mode 100644 index 00000000..1944cf4d --- /dev/null +++ b/src/solvers/subspace/densitymatrix.jl @@ -0,0 +1,74 @@ +import ITensors as it +import ITensorNetworks as itn +using NamedGraphs.GraphsExtensions: incident_edges + +function subspace_expand( + ::Backend"densitymatrix", + problem, + local_state::ITensor, + region_iterator; + north_pass=1, + expansion_factor=default_expansion_factor(), + max_expand=default_max_expand(), + trunc, + kws..., +) + region = current_region(region_iterator) + psi = copy(state(problem)) + + prev_vertex_set = setdiff(itn.pos(operator(problem)), region) + (length(prev_vertex_set) != 1) && return problem, local_state + prev_vertex = only(prev_vertex_set) + A = psi[prev_vertex] + + next_vertices = filter(v -> (it.hascommoninds(psi[v], A)), region) + isempty(next_vertices) && return problem, local_state + next_vertex = only(next_vertices) + C = psi[next_vertex] + + a = commonind(A, C) + isnothing(a) && return problem, local_state + basis_size = prod(dim.(uniqueinds(A, C))) + + expanded_maxdim = compute_expansion( + dim(a), basis_size; expansion_factor, max_expand, trunc.maxdim + ) + expanded_maxdim <= 0 && return problem, local_state + trunc = (; trunc..., maxdim=expanded_maxdim) + + envs = itn.environments(operator(problem)) + H = itn.operator(operator(problem)) + sqrt_rho = A + for e in itn.incident_edges(operator(problem)) + (src(e) ∈ region || dst(e) ∈ region) && continue + sqrt_rho *= envs[e] + end + sqrt_rho *= H[prev_vertex] + + conj_proj_A(T) = (T - prime(A)*(dag(prime(A))*T)) + for pass in 1:north_pass + sqrt_rho = conj_proj_A(sqrt_rho) + end + rho = sqrt_rho * dag(noprime(sqrt_rho)) + D, U = eigen(rho; trunc..., ishermitian=true) + + Uproj(T) = (T - prime(A, a)*(dag(prime(A, a))*T)) + for pass in 1:north_pass + U = Uproj(U) + end + if norm(dag(U)*A) > 1E-10 + @printf("Warning: |U*A| = %.3E in subspace expansion\n", norm(dag(U)*A)) + return problem, local_state + end + + Ax, ax = directsum(A=>a, U=>commonind(U, D)) + #println("Old space: ", space(a)) + #println("New space: ", space(ax)) + #ITensors.pause() + expander = dag(Ax) * A + psi[prev_vertex] = Ax + psi[next_vertex] = expander * C + local_state = expander*local_state + + return setproperties(problem; state=psi), local_state +end diff --git a/src/solvers/subspace/ortho_subspace.jl b/src/solvers/subspace/ortho_subspace.jl new file mode 100644 index 00000000..7d7ca6c2 --- /dev/null +++ b/src/solvers/subspace/ortho_subspace.jl @@ -0,0 +1,77 @@ +using ITensors +using Graphs: AbstractEdge + +expand_space(χ::Integer, expansion_factor) = max(χ + 1, floor(Int, expansion_factor * χ)) + +function expand_space(χs::Vector{<:Pair}, expansion_factor) + return [q => expand_space(d, expansion_factor) for (q, d) in χs] +end + +# +# Alternative idea for "ortho" method: +# - Just make a random tensor with `basis_inds` on both sides +# (Kind of like a random density matrix.) +# - Symmetrize to make it Hermitian PSD. +# - Then do eigenvalue decomposition to get U at desired size. +# - Finally, project out space of A from U. +# + +function subspace_expand!( + ::Backend"ortho", + problem::EigsolveProblem, + local_tensor, + region_iterator; + cutoff=default_cutoff(), + maxdim=default_maxdim(), + mindim=default_mindim(), + expansion_factor=default_expansion_factor(), + max_expand=default_max_expand(), + kws..., +) + prev_region = previous_region(region_iterator) + region = current_region(region_iterator) + if isnothing(prev_region) || isa(region, AbstractEdge) + return local_tensor + end + + prev_vertex_set = setdiff(prev_region, region) + (length(prev_vertex_set) != 1) && return local_tensor + prev_vertex = only(prev_vertex_set) + + psi = state(problem) + A = psi[prev_vertex] + + next_vertices = filter(v -> (it.hascommoninds(psi[v], A)), region) + isempty(next_vertices) && return local_tensor + next_vertex = only(next_vertices) + C = psi[next_vertex] + + # Analyze indices of A + # TODO: if "a" is missing, could supply a 1-dim index and put on both A and C? + a = commonind(A, C) + isnothing(a) && return local_tensor + basis_inds = uniqueinds(A, C) + basis_size = prod(dim.(basis_inds)) + + ci = combinedind(combiner(basis_inds...)) + ax_space = expand_space(space(ci), expansion_factor) + ax = Index(ax_space, "ax") + + linear_map(w) = (w - A * (dag(A) * w)) + Y = linear_map(random_itensor(basis_inds, dag(ax))) + expand_maxdim = compute_expansion( + dim(a), basis_size; expansion_factor, max_expand, maxdim + ) + (norm(Y) <= 1E-15 || expand_maxdim <= 0) && return local_tensor + Ux, S, V = svd(Y, basis_inds; cutoff=1E-14, maxdim=expand_maxdim, lefttags="ux,Link") + + Ux = linear_map(Ux) + ux = commonind(Ux, S) + Ax, sa = directsum(A => a, Ux => ux) + expander = dag(Ax) * A + psi[prev_vertex] = Ax + psi[next_vertex] = expander * C + local_tensor = expander*local_tensor + + return local_tensor +end diff --git a/src/solvers/subspace/subspace.jl b/src/solvers/subspace/subspace.jl new file mode 100644 index 00000000..7a332c3b --- /dev/null +++ b/src/solvers/subspace/subspace.jl @@ -0,0 +1,48 @@ +using NDTensors: NDTensors +using NDTensors.BackendSelection: Backend, @Backend_str +import ConstructionBase: setproperties + +default_expansion_factor() = 1.5 +default_max_expand() = typemax(Int) + +function subspace_expand( + problem, local_state, region_iterator; subspace_algorithm=nothing, sweep, trunc, kws... +) + return subspace_expand( + Backend(subspace_algorithm), problem, local_state, region_iterator; trunc, kws... + ) +end + +function subspace_expand(backend, problem, local_state, region_iterator; kws...) + error( + "Subspace expansion (subspace_expand!) not defined for requested combination of subspace_algorithm and problem types", + ) +end + +function subspace_expand( + backend::Backend{:nothing}, problem, local_state, region_iterator; kws... +) + problem, local_state +end + +function compute_expansion( + current_dim, + basis_size; + expansion_factor=default_expansion_factor(), + max_expand=default_max_expand(), + maxdim=default_maxdim(), +) + # Note: expand_maxdim will be *added* to current bond dimension + # Obtain expand_maxdim from expansion_factor + expand_maxdim = ceil(Int, expansion_factor * current_dim) + # Enforce max_expand keyword + expand_maxdim = min(max_expand, expand_maxdim) + + # Restrict expand_maxdim below theoretical upper limit + expand_maxdim = min(basis_size-current_dim, expand_maxdim) + # Enforce total maxdim setting (e.g. used in inserter step) + expand_maxdim = min(maxdim-current_dim, expand_maxdim) + # Ensure expand_maxdim is non-negative + expand_maxdim = max(0, expand_maxdim) + return expand_maxdim +end diff --git a/src/solvers/sweep_plans/sweep_plans.jl b/src/solvers/sweep_plans/sweep_plans.jl deleted file mode 100644 index dda6dd96..00000000 --- a/src/solvers/sweep_plans/sweep_plans.jl +++ /dev/null @@ -1,221 +0,0 @@ -using Graphs: AbstractEdge, dst, src -using NamedGraphs.GraphsExtensions: GraphsExtensions - -direction(step_number) = isodd(step_number) ? Base.Forward : Base.Reverse - -function overlap(edge_a::AbstractEdge, edge_b::AbstractEdge) - return intersect(support(edge_a), support(edge_b)) -end - -function support(edge::AbstractEdge) - return [src(edge), dst(edge)] -end - -support(r) = r - -function reverse_region(edges, which_edge; reverse_edge=false, nsites=1, region_kwargs=(;)) - current_edge = edges[which_edge] - if nsites == 1 - !reverse_edge && return [(current_edge, region_kwargs)] - reverse_edge && return [(reverse(current_edge), region_kwargs)] - elseif nsites == 2 - if last(edges) == current_edge - return () - end - future_edges = edges[(which_edge + 1):end] - future_edges = isa(future_edges, AbstractEdge) ? [future_edges] : future_edges - #error if more than single vertex overlap - overlapping_vertex = only(union([overlap(e, current_edge) for e in future_edges]...)) - return [([overlapping_vertex], region_kwargs)] - end -end - -function forward_region(edges, which_edge; nsites=1, region_kwargs=(;)) - if nsites == 1 - current_edge = edges[which_edge] - #handle edge case - if current_edge == last(edges) - overlapping_vertex = only( - union([overlap(e, current_edge) for e in edges[1:(which_edge - 1)]]...) - ) - nonoverlapping_vertex = only( - setdiff([src(current_edge), dst(current_edge)], [overlapping_vertex]) - ) - return [ - ([overlapping_vertex], region_kwargs), ([nonoverlapping_vertex], region_kwargs) - ] - else - future_edges = edges[(which_edge + 1):end] - future_edges = isa(future_edges, AbstractEdge) ? [future_edges] : future_edges - overlapping_vertex = only(union([overlap(e, current_edge) for e in future_edges]...)) - nonoverlapping_vertex = only( - setdiff([src(current_edge), dst(current_edge)], [overlapping_vertex]) - ) - return [([nonoverlapping_vertex], region_kwargs)] - end - elseif nsites == 2 - current_edge = edges[which_edge] - return [([src(current_edge), dst(current_edge)], region_kwargs)] - end -end - -function forward_sweep( - dir::Base.ForwardOrdering, - graph::AbstractGraph; - root_vertex=GraphsExtensions.default_root_vertex(graph), - reverse_edges=false, - region_kwargs, - reverse_kwargs=region_kwargs, - reverse_step=false, - kwargs..., -) - edges = post_order_dfs_edges(graph, root_vertex) - regions = map(eachindex(edges)) do i - forward_region(edges, i; region_kwargs, kwargs...) - end - regions = collect(flatten(regions)) - if reverse_step - reverse_regions = map(eachindex(edges)) do i - reverse_region( - edges, i; reverse_edge=reverse_edges, region_kwargs=reverse_kwargs, kwargs... - ) - end - reverse_regions = collect(flatten(reverse_regions)) - _check_reverse_sweeps(regions, reverse_regions, graph; kwargs...) - regions = interleave(regions, reverse_regions) - end - - return regions -end - -#ToDo: is there a better name for this? unidirectional_sweep? traversal? -function forward_sweep(dir::Base.ReverseOrdering, args...; kwargs...) - return reverse(forward_sweep(Base.Forward, args...; reverse_edges=true, kwargs...)) -end - -function default_sweep_plans( - nsweeps, - init_state; - sweep_plan_func=default_sweep_plan, - root_vertex, - extracter, - extracter_kwargs, - updater, - updater_kwargs, - inserter, - inserter_kwargs, - transform_operator, - transform_operator_kwargs, - kwargs..., -) - extracter, updater, inserter, transform_operator = extend_or_truncate.( - (extracter, updater, inserter, transform_operator), nsweeps - ) - inserter_kwargs, updater_kwargs, extracter_kwargs, transform_operator_kwargs, kwargs = expand.( - ( - inserter_kwargs, - updater_kwargs, - extracter_kwargs, - transform_operator_kwargs, - NamedTuple(kwargs), - ), - nsweeps, - ) - sweep_plans = [] - for i in 1:nsweeps - sweep_plan = sweep_plan_func( - init_state; - root_vertex, - region_kwargs=(; - inserter=inserter[i], - inserter_kwargs=inserter_kwargs[i], - updater=updater[i], - updater_kwargs=updater_kwargs[i], - extracter=extracter[i], - extracter_kwargs=extracter_kwargs[i], - transform_operator=transform_operator[i], - transform_operator_kwargs=transform_operator_kwargs[i], - ), - kwargs[i]..., - ) - push!(sweep_plans, sweep_plan) - end - return sweep_plans -end - -function default_sweep_plan( - graph::AbstractGraph; - root_vertex=GraphsExtensions.default_root_vertex(graph), - region_kwargs, - nsites::Int, -) - return vcat( - [ - forward_sweep( - direction(half), - graph; - root_vertex, - nsites, - region_kwargs=(; internal_kwargs=(; half), region_kwargs...), - ) for half in 1:2 - ]..., - ) -end - -function tdvp_sweep_plan( - graph::AbstractGraph; - root_vertex=GraphsExtensions.default_root_vertex(graph), - region_kwargs, - reverse_step=true, - order::Int, - nsites::Int, - time_step::Number, - t_evolved::Number, -) - sweep_plan = [] - for (substep, fac) in enumerate(sub_time_steps(order)) - sub_time_step = time_step * fac - append!( - sweep_plan, - forward_sweep( - direction(substep), - graph; - root_vertex, - nsites, - region_kwargs=(; - internal_kwargs=(; substep, time_step=sub_time_step, t=t_evolved), - region_kwargs..., - ), - reverse_kwargs=(; - internal_kwargs=(; substep, time_step=(-sub_time_step), t=t_evolved), - region_kwargs..., - ), - reverse_step, - ), - ) - end - return sweep_plan -end - -#ToDo: Move to test. -function _check_reverse_sweeps(forward_sweep, reverse_sweep, graph; nsites, kwargs...) - fw_regions = first.(forward_sweep) - bw_regions = first.(reverse_sweep) - if nsites == 2 - fw_verts = flatten(fw_regions) - bw_verts = flatten(bw_regions) - for v in vertices(graph) - @assert isone(count(isequal(v), fw_verts) - count(isequal(v), bw_verts)) - end - elseif nsites == 1 - fw_verts = flatten(fw_regions) - bw_edges = bw_regions - for v in vertices(graph) - @assert isone(count(isequal(v), fw_verts)) - end - for e in edges(graph) - @assert isone(count(x -> (isequal(x, e) || isequal(x, reverse(e))), bw_edges)) - end - end - return true -end diff --git a/src/solvers/sweep_solve.jl b/src/solvers/sweep_solve.jl new file mode 100644 index 00000000..f5c5862d --- /dev/null +++ b/src/solvers/sweep_solve.jl @@ -0,0 +1,40 @@ + +default_region_callback(problem; kws...) = nothing + +default_sweep_callback(problem; kws...) = nothing + +function default_sweep_printer(problem; outputlevel, sweep, nsweeps, kws...) + if outputlevel >= 1 + println("Done with sweep $sweep/$nsweeps") + end +end + +function sweep_solve( + sweep_iterator; + outputlevel=0, + region_callback=default_region_callback, + sweep_callback=default_sweep_callback, + sweep_printer=default_sweep_printer, + kwargs..., +) + for (sweep, region_iter) in enumerate(sweep_iterator) + for (region, region_kwargs) in region_tuples(region_iter) + region_callback( + problem(region_iter); + nsweeps=length(sweep_iterator), + outputlevel, + region, + region_kwargs, + sweep, + kwargs..., + ) + end + sweep_callback( + region_iter; nsweeps=length(sweep_iterator), outputlevel, sweep, kwargs... + ) + sweep_printer( + region_iter; nsweeps=length(sweep_iterator), outputlevel, sweep, kwargs... + ) + end + return problem(sweep_iterator) +end diff --git a/src/solvers/truncation_parameters.jl b/src/solvers/truncation_parameters.jl new file mode 100644 index 00000000..bcd7d940 --- /dev/null +++ b/src/solvers/truncation_parameters.jl @@ -0,0 +1,14 @@ +default_maxdim() = typemax(Int) +default_mindim() = 1 +default_cutoff() = 0.0 + +get_or_last(x, i::Integer) = (i >= length(x)) ? last(x) : x[i] + +function truncation_parameters( + sweep; cutoff=default_cutoff(), maxdim=default_maxdim(), mindim=default_mindim() +) + cutoff = get_or_last(cutoff, sweep) + mindim = get_or_last(mindim, sweep) + maxdim = get_or_last(maxdim, sweep) + return (; cutoff, mindim, maxdim) +end From 8930d8042364a7a75ad576f6f1a62514abb03c46 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Jul 2025 16:18:28 -0400 Subject: [PATCH 04/28] Remove previous solver includes --- src/ITensorNetworks.jl | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 0ca52a12..90463df3 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -34,11 +34,7 @@ include("formnetworks/quadraticformnetwork.jl") include("gauging.jl") include("utils.jl") include("update_observer.jl") -include("solvers/local_solvers/eigsolve.jl") -include("solvers/local_solvers/exponentiate.jl") -include("solvers/local_solvers/dmrg_x.jl") -include("solvers/local_solvers/contract.jl") -include("solvers/local_solvers/linsolve.jl") + include("treetensornetworks/abstracttreetensornetwork.jl") include("treetensornetworks/treetensornetwork.jl") include("treetensornetworks/opsum_to_ttn/matelem.jl") @@ -48,18 +44,7 @@ include("treetensornetworks/projttns/abstractprojttn.jl") include("treetensornetworks/projttns/projttn.jl") include("treetensornetworks/projttns/projttnsum.jl") include("treetensornetworks/projttns/projouterprodttn.jl") -include("solvers/solver_utils.jl") -include("solvers/defaults.jl") -include("solvers/insert/insert.jl") -include("solvers/extract/extract.jl") -include("solvers/alternating_update/alternating_update.jl") -include("solvers/alternating_update/region_update.jl") -include("solvers/tdvp.jl") -include("solvers/dmrg.jl") -include("solvers/dmrg_x.jl") -include("solvers/contract.jl") -include("solvers/linsolve.jl") -include("solvers/sweep_plans/sweep_plans.jl") + include("apply.jl") include("inner.jl") include("normalize.jl") From 9235349a0814e7474723348ab8649dbf3943acb2 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Jul 2025 16:28:04 -0400 Subject: [PATCH 05/28] Add ConstructionBase dependency --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index c3d97358..12c8898f 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.13.12" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" @@ -53,6 +54,7 @@ AbstractTrees = "0.4.4" Adapt = "4" Combinatorics = "1" Compat = "3, 4" +ConstructionBase = "1.6.0" DataGraphs = "0.2.3" DataStructures = "0.18" Dictionaries = "0.4" From d70bf03ef4ce2b31e364f57dd100d7fcce782bda Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Jul 2025 17:08:55 -0400 Subject: [PATCH 06/28] Adapt new solver codes into ITensorNetworks module --- src/ITensorNetworks.jl | 17 +++ src/solvers/applyexp.jl | 16 +-- src/solvers/eigsolve.jl | 10 +- src/solvers/extracter.jl | 4 +- src/solvers/fitting.jl | 112 -------------------- src/solvers/inserter.jl | 8 +- src/solvers/operator_map.jl | 43 +------- src/solvers/permute_indices.jl | 6 +- src/solvers/region_plans.jl | 6 +- src/solvers/subspace.jl | 22 +++- src/solvers/subspace/densitymatrix.jl | 13 ++- src/solvers/tdvp.jl | 36 ------- test/solvers/test_tree_dmrg.jl | 68 ++++++++++++ test/solvers/utilities/simple_ed_methods.jl | 42 ++++++++ test/solvers/utilities/tree_graphs.jl | 23 ++++ 15 files changed, 203 insertions(+), 223 deletions(-) delete mode 100644 src/solvers/fitting.jl delete mode 100644 src/solvers/tdvp.jl create mode 100644 test/solvers/test_tree_dmrg.jl create mode 100644 test/solvers/utilities/simple_ed_methods.jl create mode 100644 test/solvers/utilities/tree_graphs.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 90463df3..02d4c25b 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -35,6 +35,23 @@ include("gauging.jl") include("utils.jl") include("update_observer.jl") +include("solvers/local_solvers/eigsolve.jl") +include("solvers/local_solvers/exponentiate.jl") +include("solvers/local_solvers/runge_kutta.jl") +include("solvers/truncation_parameters.jl") +include("solvers/sweep_solve.jl") +include("solvers/iterators.jl") +include("solvers/region_plans/dfs_plans.jl") +include("solvers/region_plans/euler_tour.jl") +include("solvers/region_plans/euler_plans.jl") +include("solvers/region_plans/tdvp_region_plans.jl") +include("solvers/extracter.jl") +include("solvers/inserter.jl") +include("solvers/subspace/subspace.jl") +include("solvers/subspace/densitymatrix.jl") +include("solvers/eigsolve.jl") +include("solvers/applyexp.jl") + include("treetensornetworks/abstracttreetensornetwork.jl") include("treetensornetworks/treetensornetwork.jl") include("treetensornetworks/opsum_to_ttn/matelem.jl") diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 8832cfcb..2cb76024 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -1,5 +1,5 @@ -import ITensorNetworks as itn -using Printf +using Printf: @printf +import ConstructionBase: setproperties @kwdef mutable struct ApplyExpProblem{State} state::State @@ -7,7 +7,7 @@ using Printf current_time::Number = 0.0 end -state(tdvp::ApplyExpProblem) = tdvp.state +ITensorNetworks.state(tdvp::ApplyExpProblem) = tdvp.state operator(tdvp::ApplyExpProblem) = tdvp.operator current_time(tdvp::ApplyExpProblem) = tdvp.current_time @@ -31,11 +31,11 @@ function updater( curr_reg = current_region(region_iterator) next_reg = next_region(region_iterator) if !isnothing(next_reg) && next_reg != curr_reg - next_edge = first(itn.edge_sequence_between_regions(state(T), curr_reg, next_reg)) - v1, v2 = itn.src(next_edge), itn.dst(next_edge) + next_edge = first(edge_sequence_between_regions(state(T), curr_reg, next_reg)) + v1, v2 = src(next_edge), dst(next_edge) psi = copy(state(T)) psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2])) - shifted_operator = itn.position(operator(T), psi, itn.NamedEdge(v1=>v2)) + shifted_operator = position(operator(T), psi, NamedEdge(v1=>v2)) R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...) local_state = psi[v1]*R_t end @@ -53,7 +53,7 @@ function applyexp_sweep_printer( if outputlevel >= 1 T = problem(region_iterator) @printf(" Current time = %s, ", process_time(current_time(T))) - @printf("maxlinkdim=%d", itn.maxlinkdim(state(T))) + @printf("maxlinkdim=%d", maxlinkdim(state(T))) println() flush(stdout) end @@ -83,7 +83,7 @@ end function applyexp(H, init_state, exponents; kws...) init_prob = ApplyExpProblem(; - state=permute_indices(init_state), operator=itn.ProjTTN(permute_indices(H)) + state=permute_indices(init_state), operator=ProjTTN(permute_indices(H)) ) return applyexp(init_prob, exponents; kws...) end diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index d2b0b042..0edc26dd 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,5 +1,5 @@ -import ITensorNetworks as itn -using Printf +using Printf: @printf +import ConstructionBase: setproperties @kwdef mutable struct EigsolveProblem{State,Operator} state::State @@ -8,7 +8,7 @@ using Printf end eigenvalue(E::EigsolveProblem) = E.eigenvalue -state(E::EigsolveProblem) = E.state +ITensorNetworks.state(E::EigsolveProblem) = E.state operator(E::EigsolveProblem) = E.operator function updater( @@ -36,7 +36,7 @@ function eigsolve_sweep_printer(region_iterator; outputlevel, sweep, nsweeps, kw end E = problem(region_iterator) @printf("eigenvalue=%.12f ", eigenvalue(E)) - @printf("maxlinkdim=%d", itn.maxlinkdim(state(E))) + @printf("maxlinkdim=%d", maxlinkdim(state(E))) println() flush(stdout) end @@ -68,7 +68,7 @@ end function eigsolve(H, init_state; kws...) init_prob = EigsolveProblem(; - state=permute_indices(init_state), operator=itn.ProjTTN(permute_indices(H)) + state=permute_indices(init_state), operator=ProjTTN(permute_indices(H)) ) return eigsolve(init_prob; kws...) end diff --git a/src/solvers/extracter.jl b/src/solvers/extracter.jl index 50a396f6..1bf9f127 100644 --- a/src/solvers/extracter.jl +++ b/src/solvers/extracter.jl @@ -3,7 +3,7 @@ import ConstructionBase: setproperties function extracter(problem, region_iterator; sweep, trunc=(;), kws...) trunc = truncation_parameters(sweep; trunc...) region = current_region(region_iterator) - psi = itn.orthogonalize(state(problem), region) + psi = orthogonalize(state(problem), region) local_state = prod(psi[v] for v in region) problem = setproperties(problem; state=psi) @@ -11,7 +11,7 @@ function extracter(problem, region_iterator; sweep, trunc=(;), kws...) problem, local_state, region_iterator; sweep, trunc, kws... ) - shifted_operator = itn.position(operator(problem), state(problem), region) + shifted_operator = position(operator(problem), state(problem), region) return setproperties(problem; operator=shifted_operator), local_state end diff --git a/src/solvers/fitting.jl b/src/solvers/fitting.jl deleted file mode 100644 index 0995124b..00000000 --- a/src/solvers/fitting.jl +++ /dev/null @@ -1,112 +0,0 @@ -import ITensorNetworks as itn -import ITensors as it -import Graphs: vertices -import NamedGraphs.PartitionedGraphs as npg -using NamedGraphs: AbstractNamedGraph, NamedEdge -using Printf - -@kwdef mutable struct FittingProblem{State<:itn.AbstractBeliefPropagationCache} - state::State - ket_graph::AbstractNamedGraph - overlap::Number = 0 - gauge_region -end - -overlap(F::FittingProblem) = F.overlap -state(F::FittingProblem) = F.state -ket_graph(F::FittingProblem) = F.ket_graph -gauge_region(F::FittingProblem) = F.gauge_region - -function ket(F::FittingProblem) - ket_vertices = vertices(ket_graph(F)) - return first(itn.induced_subgraph(itn.tensornetwork(state(F)), ket_vertices)) -end - -function extracter(problem::FittingProblem, region_iterator; sweep, kws...) - region = current_region(region_iterator) - prev_region = gauge_region(problem) - tn = state(problem) - path = itn.edge_sequence_between_regions(ket_graph(problem), prev_region, region) - tn = itn.gauge_walk(itn.Algorithm("orthogonalize"), tn, path) - pe_path = npg.partitionedges(itn.partitioned_tensornetwork(tn), path) - tn = itn.update( - itn.Algorithm("bp"), tn, pe_path; message_update_function_kwargs=(; normalize=false) - ) - local_tensor = itn.environment(tn, region) - sequence = itn.contraction_sequence(local_tensor; alg="optimal") - local_tensor = dag(it.contract(local_tensor; sequence)) - #problem, local_tensor = subspace_expand(problem, local_tensor, region; sweep, kws...) - return setproperties(problem; state=tn, gauge_region=region), local_tensor -end - -function updater(F::FittingProblem, local_tensor, region; outputlevel, kws...) - n = (local_tensor * dag(local_tensor))[] - F = setproperties(F; overlap=n / sqrt(n)) - if outputlevel >= 2 - @printf(" Region %s: squared overlap = %.12f\n", region, overlap(F)) - end - return F, local_tensor -end - -function region_plan(F::FittingProblem; nsites, sweep_kwargs...) - return euler_sweep(ket_graph(F); nsites, sweep_kwargs...) -end - -function fit_tensornetwork( - overlap_network, - args...; - nsweeps=25, - nsites=1, - outputlevel=0, - extracter_kwargs=(;), - updater_kwargs=(;), - inserter_kwargs=(;), - normalize=true, - kws..., -) - bpc = itn.BeliefPropagationCache(overlap_network, args...) - ket_vertices = itn.ket_vertices(overlap_network) - ket_graph = first( - itn.induced_subgraph( - itn.underlying_graph(overlap_network), itn.ket_vertices(overlap_network) - ), - ) - init_prob = FittingProblem(; - ket_graph=ket_graph, state=bpc, gauge_region=collect(vertices(ket_graph)) - ) - - inserter_kwargs = (; inserter_kwargs..., normalize, set_orthogonal_region=false) - common_sweep_kwargs = (; nsites, outputlevel, updater_kwargs, inserter_kwargs) - kwargs_array = [(; common_sweep_kwargs..., sweep=s) for s in 1:nsweeps] - sweep_iter = sweep_iterator(init_prob, kwargs_array) - converged_prob = sweep_solve(sweep_iter; outputlevel, kws...) - return itn.rename_vertices(itn.inv_vertex_map(overlap_network), ket(converged_prob)) -end - -function fit_tensornetwork(tn, init_state, args...; kwargs...) - return fit_tensornetwork(itn.inner_network(tn, init_state), args; kwargs...) -end - -function itn.truncate(tn; maxdim=default_maxdim(), cutoff=default_cutoff(), kwargs...) - init_state = itn.ITensorNetwork( - v -> inds -> it.delta(inds), itn.siteinds(tn); link_space=maxdim - ) - overlap_network = itn.inner_network(tn, init_state) - inserter_kwargs = (; trunc=(; cutoff, maxdim)) - return fit_tensornetwork(overlap_network; inserter_kwargs, kwargs...) -end - -function itn.apply( - A::itn.ITensorNetwork, - x::itn.ITensorNetwork; - maxdim=default_maxdim(), - cutoff=default_cutoff(), - kwargs..., -) - init_state = itn.ITensorNetwork( - v -> inds -> it.delta(inds), itn.siteinds(x); link_space=maxdim - ) - overlap_network = itn.inner_network(x, A, init_state) - inserter_kwargs = (; trunc=(; cutoff, maxdim)) - return fit_tensornetwork(overlap_network; inserter_kwargs, kwargs...) -end diff --git a/src/solvers/inserter.jl b/src/solvers/inserter.jl index 64a03921..770f231b 100644 --- a/src/solvers/inserter.jl +++ b/src/solvers/inserter.jl @@ -21,13 +21,13 @@ function inserter( indsTe = it.inds(psi[first(region)]) tags = it.tags(psi, e) U, C, _ = it.factorize(local_tensor, indsTe; tags, trunc...) - itn.@preserve_graph psi[first(region)] = U + @preserve_graph psi[first(region)] = U else error("Region of length $(length(region)) not currently supported") end v = last(region) - itn.@preserve_graph psi[v] = C - psi = set_orthogonal_region ? itn.set_ortho_region(psi, [v]) : psi - normalize && itn.@preserve_graph psi[v] = psi[v] / norm(psi[v]) + @preserve_graph psi[v] = C + psi = set_orthogonal_region ? set_ortho_region(psi, [v]) : psi + normalize && @preserve_graph psi[v] = psi[v] / norm(psi[v]) return setproperties(problem; state=psi) end diff --git a/src/solvers/operator_map.jl b/src/solvers/operator_map.jl index 309d85ad..78a6f264 100644 --- a/src/solvers/operator_map.jl +++ b/src/solvers/operator_map.jl @@ -1,42 +1,9 @@ -import ITensorNetworks as itn -function optimal_map(P::itn.ProjTTN, ψ) - envs = [itn.environment(P, e) for e in itn.incident_edges(P)] - site_ops = [itn.operator(P)[s] for s in itn.sites(P)] +function optimal_map(P::ProjTTN, ψ) + envs = [environment(P, e) for e in incident_edges(P)] + site_ops = [operator(P)[s] for s in sites(P)] contract_list = [envs..., site_ops..., ψ] - sequence = itn.contraction_sequence(contract_list; alg="optimal") - Pψ = itn.contract(contract_list; sequence) + sequence = contraction_sequence(contract_list; alg="optimal") + Pψ = contract(contract_list; sequence) return noprime(Pψ) end - -# This function is a workaround for the slow contraction order -# heuristic in ITensorNetworks/src/treetensornetworks/projttns/projttn.jl -# in the projected_operator_tensors(P::ProjTTN) function (line 97 or so) -function operator_map(P::itn.ProjTTN, ψ) - ψ = copy(ψ) - if itn.on_edge(P) - for edge in itn.incident_edges(P) - ψ *= itn.environment(P, edge) - end - else - region = itn.sites(P) - ie = itn.incident_edges(P) - # TODO: improvement ideas - # - check which vertex (first(region) vs. last(region) - # has more incident edges and contract those environments first - for edge in ie - if itn.dst(edge) == first(region) - ψ *= itn.environment(P, edge) - end - end - for s in itn.sites(P) - ψ *= itn.operator(P)[s] - end - for edge in ie - if itn.dst(edge) != first(region) - ψ *= itn.environment(P, edge) - end - end - end - return noprime(ψ) -end diff --git a/src/solvers/permute_indices.jl b/src/solvers/permute_indices.jl index 927f6b87..90526965 100644 --- a/src/solvers/permute_indices.jl +++ b/src/solvers/permute_indices.jl @@ -1,10 +1,8 @@ -import ITensorNetworks as itn -using ITensors function permute_indices(tn) - si = itn.siteinds(tn) + si = siteinds(tn) ptn = copy(tn) - for v in itn.vertices(tn) + for v in vertices(tn) is = inds(tn[v]) ls = setdiff(is, si[v]) isempty(ls) && continue diff --git a/src/solvers/region_plans.jl b/src/solvers/region_plans.jl index ac4ab733..be858dfd 100644 --- a/src/solvers/region_plans.jl +++ b/src/solvers/region_plans.jl @@ -6,7 +6,9 @@ import NamedGraphs: GraphsExtensions # return [fwd_sweep..., reverse(fwd_sweep)...] #end -function tdvp_regions(g::AbstractGraph, time_step; nsites=1, updater_kwargs, sweep_kwargs...) +function tdvp_regions( + g::AbstractGraph, time_step; nsites=1, updater_kwargs, sweep_kwargs... +) @assert nsites==1 fwd_up_args = (; time=(time_step / 2), updater_kwargs...) rev_up_args = (; time=(-time_step / 2), updater_kwargs...) @@ -25,7 +27,7 @@ function tdvp_regions(g::AbstractGraph, time_step; nsites=1, updater_kwargs, swe end function overlap(ea::AbstractEdge, eb::AbstractEdge) - return intersect([src(ea),dst(ea)], [src(eb),dst(eb)]) + return intersect([src(ea), dst(ea)], [src(eb), dst(eb)]) end function forward_region(edges, which_edge; nsites=1, region_kwargs=(;)) diff --git a/src/solvers/subspace.jl b/src/solvers/subspace.jl index 43eb80bb..7fccc8f0 100644 --- a/src/solvers/subspace.jl +++ b/src/solvers/subspace.jl @@ -1,9 +1,21 @@ -using ITensors: commonind, dag, dim, directsum, dot, hascommoninds, Index, norm, onehot, uniqueinds, random_itensor +using ITensors: + commonind, + dag, + dim, + directsum, + dot, + hascommoninds, + Index, + norm, + onehot, + uniqueinds, + random_itensor # TODO: hoist num_expand default value out to a function or similar -function subspace_expand!(problem::EigsolveProblem, local_tensor, region; prev_region, num_expand=4, kws...) - - if isnothing(prev_region) || isa(region, AbstractEdge) +function subspace_expand!( + problem::EigsolveProblem, local_tensor, region; prev_region, num_expand=4, kws... +) + if isnothing(prev_region) || isa(region, AbstractEdge) return local_tensor end @@ -49,4 +61,4 @@ function subspace_expand!(problem::EigsolveProblem, local_tensor, region; prev_r local_tensor = prod(psi[v] for v in region) return local_tensor -end +end diff --git a/src/solvers/subspace/densitymatrix.jl b/src/solvers/subspace/densitymatrix.jl index 1944cf4d..7511901d 100644 --- a/src/solvers/subspace/densitymatrix.jl +++ b/src/solvers/subspace/densitymatrix.jl @@ -1,6 +1,5 @@ -import ITensors as it -import ITensorNetworks as itn using NamedGraphs.GraphsExtensions: incident_edges +using Printf: @printf function subspace_expand( ::Backend"densitymatrix", @@ -16,12 +15,12 @@ function subspace_expand( region = current_region(region_iterator) psi = copy(state(problem)) - prev_vertex_set = setdiff(itn.pos(operator(problem)), region) + prev_vertex_set = setdiff(pos(operator(problem)), region) (length(prev_vertex_set) != 1) && return problem, local_state prev_vertex = only(prev_vertex_set) A = psi[prev_vertex] - next_vertices = filter(v -> (it.hascommoninds(psi[v], A)), region) + next_vertices = filter(v -> (hascommoninds(psi[v], A)), region) isempty(next_vertices) && return problem, local_state next_vertex = only(next_vertices) C = psi[next_vertex] @@ -36,10 +35,10 @@ function subspace_expand( expanded_maxdim <= 0 && return problem, local_state trunc = (; trunc..., maxdim=expanded_maxdim) - envs = itn.environments(operator(problem)) - H = itn.operator(operator(problem)) + envs = environments(operator(problem)) + H = operator(operator(problem)) sqrt_rho = A - for e in itn.incident_edges(operator(problem)) + for e in incident_edges(operator(problem)) (src(e) ∈ region || dst(e) ∈ region) && continue sqrt_rho *= envs[e] end diff --git a/src/solvers/tdvp.jl b/src/solvers/tdvp.jl deleted file mode 100644 index db5a36e8..00000000 --- a/src/solvers/tdvp.jl +++ /dev/null @@ -1,36 +0,0 @@ -import ITensorNetworks as itn -using Printf - -@kwdef mutable struct TDVPProblem{State} - state::State - operator - #current_time::Number = 0.0 -end - -state(tdvp::TDVPProblem) = tdvp.state -operator(tdvp::TDVPProblem) = tdvp.operator - -function set(tdvp::TDVPProblem; state=state(tdvp), operator=operator(tdvp)) - return TDVPProblem(; state, operator) -end - -function region_plan(tdvp::TDVPProblem; nsites, time_step, sweep_kwargs...) - return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) -end - -function updater!(tdvp::TDVPProblem, local_tensor, region; outputlevel, kws...) - local_tensor, info = exponentiate_updater(operator(tdvp), local_tensor; kws...) - return local_tensor -end - -function applyexp( - H, init_state, time_points; updater_kwargs=(;), inserter_kwargs=(;), outputlevel=0, kws... -) - init_prob = TDVPProblem(; state=copy(init_state), operator=itn.ProjTTN(H)) - time_steps = diff([0.0, time_points...]) - common_sweep_kwargs = (; outputlevel, updater_kwargs, inserter_kwargs) - kwargs_array = [(; common_sweep_kwargs..., time_step=t) for t in time_steps] - sweep_iter = sweep_iterator(init_prob, kwargs_array) - converged_prob = alternating_update(sweep_iter; outputlevel, kws...) - return state(converged_prob) -end diff --git a/test/solvers/test_tree_dmrg.jl b/test/solvers/test_tree_dmrg.jl new file mode 100644 index 00000000..d5f610d0 --- /dev/null +++ b/test/solvers/test_tree_dmrg.jl @@ -0,0 +1,68 @@ +import NetworkSolvers as ns +using Test: @test, @testset + +using ITensors +import ITensorNetworks as itn +import Graphs as gr +import NamedGraphs as ng +import ITensorMPS as itm + +include("utilities/simple_ed_methods.jl") +include("utilities/tree_graphs.jl") + +@testset "Tree DMRG" begin + outputlevel = 1 + + g = build_tree(; nbranch=3, nbranch_sites=3) + + sites = itn.siteinds("S=1/2", g) + + # Make Heisenberg model Hamiltonian + h = itm.OpSum() + for edge in gr.edges(sites) + i, j = gr.src(edge), gr.dst(edge) + h += "Sz", i, "Sz", j + h += 1/2, "S+", i, "S-", j + h += 1/2, "S-", i, "S+", j + end + H = itn.ttn(h, sites) + + # Make initial product state + state = Dict{Tuple{Int,Int},String}() + for (j, v) in enumerate(gr.vertices(sites)) + state[v] = iseven(j) ? "Up" : "Dn" + end + psi0 = itn.ttn(state, sites) + + (outputlevel >= 1) && println("Computing exact ground state") + Ex, psix = ed_ground_state(H, psi0) + (outputlevel >= 1) && println("Ex = ", Ex) + + cutoff = 1E-5 + maxdim = 40 + nsweeps = 5 + + # + # Test 2-site DMRG without subspace expansion + # + nsites = 2 + trunc = (; cutoff, maxdim) + inserter_kwargs = (; trunc) + E, psi = ns.dmrg(H, psi0; inserter_kwargs, nsites, nsweeps, outputlevel) + (outputlevel >= 1) && println("2-site DMRG energy = ", E) + @test abs(E-Ex) < 1E-5 + + # + # Test 1-site DMRG with subspace expansion + # + nsites = 1 + nsweeps = 5 + trunc = (; cutoff, maxdim) + extracter_kwargs = (; trunc, subspace_algorithm="densitymatrix") + inserter_kwargs = (; trunc) + cutoff = 1E-10 + maxdim = 200 + E, psi = ns.dmrg(H, psi0; extracter_kwargs, inserter_kwargs, nsites, nsweeps, outputlevel) + (outputlevel >= 1) && println("1-site+subspace DMRG energy = ", E) + @test abs(E-Ex) < 1E-5 +end diff --git a/test/solvers/utilities/simple_ed_methods.jl b/test/solvers/utilities/simple_ed_methods.jl new file mode 100644 index 00000000..19e58109 --- /dev/null +++ b/test/solvers/utilities/simple_ed_methods.jl @@ -0,0 +1,42 @@ +import ITensorNetworks as itn +import NetworkSolvers as ns + +function ed_ground_state(H, psi0) + ITensors.disable_warn_order() + H = prod(H) + psi = prod(psi0) + # TODO: call KrylovKit instead + expH = exp(H*(-20.0)) + for napply in 1:10 + psi = noprime(expH*psi) + psi /= norm(psi) + end + E = scalar(prime(psi)*H*psi) + return E, psi +end + +function ed_time_evolution( + H::itn.AbstractITensorNetwork, + psi::itn.AbstractITensorNetwork, + time_points; + normalize=false, +) + ITensors.disable_warn_order() + H = prod(H) + psi = prod(psi) + exponents = [-im*t for t in time_points] + steps = diff([0.0, exponents...])[2:end] + H_map = ψ -> noprime(H*psi) + for step in steps + expH = exp(H * step) + psi = noprime(expH * psi) + #psir, _ = ns.runge_kutta_solver(H_map,step,psi; order=4) + #psix, _ = ns.exponentiate_solver(H_map,step,psi) + #@show abs(scalar(dag(psi)*psix)) + #@show abs(scalar(dag(psir)*psix)) + if normalize + psi /= norm(psi) + end + end + return psi +end diff --git a/test/solvers/utilities/tree_graphs.jl b/test/solvers/utilities/tree_graphs.jl new file mode 100644 index 00000000..4fdff761 --- /dev/null +++ b/test/solvers/utilities/tree_graphs.jl @@ -0,0 +1,23 @@ +import Graphs as gr +import NamedGraphs as ng + +""" + build_tree + + Make a tree with central vertex (0,0) and + nbranch branches of nbranch_sites each. +""" +function build_tree(; nbranch=3, nbranch_sites=3) + g = ng.NamedGraph() + gr.add_vertex!(g, (0, 0)) + for branch in 1:nbranch, site in 1:nbranch_sites + gr.add_vertex!(g, (branch, site)) + end + for branch in 1:nbranch + gr.add_edge!(g, (0, 0)=>(branch, 1)) + for site in 2:nbranch_sites + gr.add_edge!(g, (branch, site-1)=>(branch, site)) + end + end + return g +end From c081dfb6efbb68ac63424749589e2da585d906ed Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Jul 2025 22:41:49 -0400 Subject: [PATCH 07/28] Continue adapting code and improve DMRG test --- src/ITensorNetworks.jl | 25 ++++++++++-------- src/solvers/inserter.jl | 9 ++++--- test/solvers/test_tree_dmrg.jl | 28 +++++++++------------ test/solvers/utilities/simple_ed_methods.jl | 13 ++-------- 4 files changed, 33 insertions(+), 42 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 02d4c25b..582a51e2 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -35,12 +35,23 @@ include("gauging.jl") include("utils.jl") include("update_observer.jl") +include("treetensornetworks/abstracttreetensornetwork.jl") +include("treetensornetworks/treetensornetwork.jl") +include("treetensornetworks/opsum_to_ttn/matelem.jl") +include("treetensornetworks/opsum_to_ttn/qnarrelem.jl") +include("treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl") +include("treetensornetworks/projttns/abstractprojttn.jl") +include("treetensornetworks/projttns/projttn.jl") +include("treetensornetworks/projttns/projttnsum.jl") +include("treetensornetworks/projttns/projouterprodttn.jl") + include("solvers/local_solvers/eigsolve.jl") include("solvers/local_solvers/exponentiate.jl") include("solvers/local_solvers/runge_kutta.jl") include("solvers/truncation_parameters.jl") -include("solvers/sweep_solve.jl") include("solvers/iterators.jl") +include("solvers/adapters.jl") +include("solvers/sweep_solve.jl") include("solvers/region_plans/dfs_plans.jl") include("solvers/region_plans/euler_tour.jl") include("solvers/region_plans/euler_plans.jl") @@ -49,19 +60,11 @@ include("solvers/extracter.jl") include("solvers/inserter.jl") include("solvers/subspace/subspace.jl") include("solvers/subspace/densitymatrix.jl") +include("solvers/permute_indices.jl") +include("solvers/operator_map.jl") include("solvers/eigsolve.jl") include("solvers/applyexp.jl") -include("treetensornetworks/abstracttreetensornetwork.jl") -include("treetensornetworks/treetensornetwork.jl") -include("treetensornetworks/opsum_to_ttn/matelem.jl") -include("treetensornetworks/opsum_to_ttn/qnarrelem.jl") -include("treetensornetworks/opsum_to_ttn/opsum_to_ttn.jl") -include("treetensornetworks/projttns/abstractprojttn.jl") -include("treetensornetworks/projttns/projttn.jl") -include("treetensornetworks/projttns/projttnsum.jl") -include("treetensornetworks/projttns/projouterprodttn.jl") - include("apply.jl") include("inner.jl") include("normalize.jl") diff --git a/src/solvers/inserter.jl b/src/solvers/inserter.jl index 770f231b..2180fb14 100644 --- a/src/solvers/inserter.jl +++ b/src/solvers/inserter.jl @@ -1,4 +1,5 @@ import ConstructionBase: setproperties +import NamedGraphs: edgetype function inserter( problem, @@ -17,10 +18,10 @@ function inserter( if length(region) == 1 C = local_tensor elseif length(region) == 2 - e = ng.edgetype(psi)(first(region), last(region)) - indsTe = it.inds(psi[first(region)]) - tags = it.tags(psi, e) - U, C, _ = it.factorize(local_tensor, indsTe; tags, trunc...) + e = edgetype(psi)(first(region), last(region)) + indsTe = inds(psi[first(region)]) + tags = ITensors.tags(psi, e) + U, C, _ = factorize(local_tensor, indsTe; tags, trunc...) @preserve_graph psi[first(region)] = U else error("Region of length $(length(region)) not currently supported") diff --git a/test/solvers/test_tree_dmrg.jl b/test/solvers/test_tree_dmrg.jl index d5f610d0..57459056 100644 --- a/test/solvers/test_tree_dmrg.jl +++ b/test/solvers/test_tree_dmrg.jl @@ -1,11 +1,9 @@ -import NetworkSolvers as ns using Test: @test, @testset - using ITensors -import ITensorNetworks as itn -import Graphs as gr -import NamedGraphs as ng -import ITensorMPS as itm +using TensorOperations # Needed to use contraction order finding +using ITensorNetworks: siteinds, ttn, dmrg +import Graphs: dst, edges, src +import ITensorMPS: OpSum include("utilities/simple_ed_methods.jl") include("utilities/tree_graphs.jl") @@ -15,24 +13,24 @@ include("utilities/tree_graphs.jl") g = build_tree(; nbranch=3, nbranch_sites=3) - sites = itn.siteinds("S=1/2", g) + sites = siteinds("S=1/2", g) # Make Heisenberg model Hamiltonian - h = itm.OpSum() - for edge in gr.edges(sites) - i, j = gr.src(edge), gr.dst(edge) + h = OpSum() + for edge in edges(sites) + i, j = src(edge), dst(edge) h += "Sz", i, "Sz", j h += 1/2, "S+", i, "S-", j h += 1/2, "S-", i, "S+", j end - H = itn.ttn(h, sites) + H = ttn(h, sites) # Make initial product state state = Dict{Tuple{Int,Int},String}() for (j, v) in enumerate(gr.vertices(sites)) state[v] = iseven(j) ? "Up" : "Dn" end - psi0 = itn.ttn(state, sites) + psi0 = ttn(state, sites) (outputlevel >= 1) && println("Computing exact ground state") Ex, psix = ed_ground_state(H, psi0) @@ -48,7 +46,7 @@ include("utilities/tree_graphs.jl") nsites = 2 trunc = (; cutoff, maxdim) inserter_kwargs = (; trunc) - E, psi = ns.dmrg(H, psi0; inserter_kwargs, nsites, nsweeps, outputlevel) + E, psi = dmrg(H, psi0; inserter_kwargs, nsites, nsweeps, outputlevel) (outputlevel >= 1) && println("2-site DMRG energy = ", E) @test abs(E-Ex) < 1E-5 @@ -60,9 +58,7 @@ include("utilities/tree_graphs.jl") trunc = (; cutoff, maxdim) extracter_kwargs = (; trunc, subspace_algorithm="densitymatrix") inserter_kwargs = (; trunc) - cutoff = 1E-10 - maxdim = 200 - E, psi = ns.dmrg(H, psi0; extracter_kwargs, inserter_kwargs, nsites, nsweeps, outputlevel) + E, psi = dmrg(H, psi0; extracter_kwargs, inserter_kwargs, nsites, nsweeps, outputlevel) (outputlevel >= 1) && println("1-site+subspace DMRG energy = ", E) @test abs(E-Ex) < 1E-5 end diff --git a/test/solvers/utilities/simple_ed_methods.jl b/test/solvers/utilities/simple_ed_methods.jl index 19e58109..c80d3be8 100644 --- a/test/solvers/utilities/simple_ed_methods.jl +++ b/test/solvers/utilities/simple_ed_methods.jl @@ -1,11 +1,9 @@ -import ITensorNetworks as itn -import NetworkSolvers as ns +import ITensorNetworks: AbstractITensorNetwork function ed_ground_state(H, psi0) ITensors.disable_warn_order() H = prod(H) psi = prod(psi0) - # TODO: call KrylovKit instead expH = exp(H*(-20.0)) for napply in 1:10 psi = noprime(expH*psi) @@ -16,10 +14,7 @@ function ed_ground_state(H, psi0) end function ed_time_evolution( - H::itn.AbstractITensorNetwork, - psi::itn.AbstractITensorNetwork, - time_points; - normalize=false, + H::AbstractITensorNetwork, psi::AbstractITensorNetwork, time_points; normalize=false ) ITensors.disable_warn_order() H = prod(H) @@ -30,10 +25,6 @@ function ed_time_evolution( for step in steps expH = exp(H * step) psi = noprime(expH * psi) - #psir, _ = ns.runge_kutta_solver(H_map,step,psi; order=4) - #psix, _ = ns.exponentiate_solver(H_map,step,psi) - #@show abs(scalar(dag(psi)*psix)) - #@show abs(scalar(dag(psir)*psix)) if normalize psi /= norm(psi) end From e5047fc5ea5cecdc4cbeda1b27b5e2c06c0f7a3b Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 1 Jul 2025 22:56:20 -0400 Subject: [PATCH 08/28] Rename tdvp to time_evolve. Add tests. --- src/solvers/applyexp.jl | 24 +++++------ test/solvers/test_tree_tdvp.jl | 76 ++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 12 deletions(-) create mode 100644 test/solvers/test_tree_tdvp.jl diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 2cb76024..01c042d7 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -7,16 +7,16 @@ import ConstructionBase: setproperties current_time::Number = 0.0 end -ITensorNetworks.state(tdvp::ApplyExpProblem) = tdvp.state -operator(tdvp::ApplyExpProblem) = tdvp.operator -current_time(tdvp::ApplyExpProblem) = tdvp.current_time +ITensorNetworks.state(A::ApplyExpProblem) = A.state +operator(A::ApplyExpProblem) = A.operator +current_time(A::ApplyExpProblem) = A.current_time function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...) return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) end function updater( - T::ApplyExpProblem, + A::ApplyExpProblem, local_state, region_iterator; nsites, @@ -25,26 +25,26 @@ function updater( outputlevel, kws..., ) - local_state, info = solver(x->optimal_map(operator(T), x), time_step, local_state; kws...) + local_state, info = solver(x->optimal_map(operator(A), x), time_step, local_state; kws...) if nsites==1 curr_reg = current_region(region_iterator) next_reg = next_region(region_iterator) if !isnothing(next_reg) && next_reg != curr_reg - next_edge = first(edge_sequence_between_regions(state(T), curr_reg, next_reg)) + next_edge = first(edge_sequence_between_regions(state(A), curr_reg, next_reg)) v1, v2 = src(next_edge), dst(next_edge) - psi = copy(state(T)) + psi = copy(state(A)) psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2])) - shifted_operator = position(operator(T), psi, NamedEdge(v1=>v2)) + shifted_operator = position(operator(A), psi, NamedEdge(v1=>v2)) R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...) local_state = psi[v1]*R_t end end - curr_time = current_time(T) + time_step - T = setproperties(T; current_time=curr_time) + curr_time = current_time(A) + time_step + A = setproperties(A; current_time=curr_time) - return T, local_state + return A, local_state end function applyexp_sweep_printer( @@ -90,7 +90,7 @@ end process_real_times(z) = round(-imag(z); digits=10) -function tdvp( +function time_evolve( H, init_state, time_points; diff --git a/test/solvers/test_tree_tdvp.jl b/test/solvers/test_tree_tdvp.jl new file mode 100644 index 00000000..04bdece4 --- /dev/null +++ b/test/solvers/test_tree_tdvp.jl @@ -0,0 +1,76 @@ +using Test: @test, @testset +using ITensors +using TensorOperations # Needed to use contraction order finding +import ITensorNetworks: dmrg, maxlinkdim, siteinds, time_evolve, ttn +import Graphs: add_vertex!, add_edge!, vertices +import NamedGraphs: NamedGraph +import ITensorMPS: OpSum + +function chain_plus_ancilla(; nchain) + g = NamedGraph() + for j in 1:nchain + add_vertex!(g, j) + end + for j in 1:(nchain - 1) + add_edge!(g, j=>j+1) + end + # Add ancilla vertex near middle of chain + add_vertex!(g, 0) + add_edge!(g, 0=>nchain÷2) + return g +end + +@testset "Tree TDVP on chain plus ancilla" begin + outputlevel = 1 + + N = 10 + g = chain_plus_ancilla(; nchain=N) + + sites = siteinds("S=1/2", g) + + # Make Heisenberg model Hamiltonian + h = OpSum() + for j in 1:(N - 1) + h += "Sz", j, "Sz", j+1 + h += 1/2, "S+", j, "S-", j+1 + h += 1/2, "S-", j, "S+", j+1 + end + H = ttn(h, sites) + + # Make initial product state + state = Dict{Int,String}() + for (j, v) in enumerate(vertices(sites)) + state[v] = iseven(j) ? "Up" : "Dn" + end + psi0 = ttn(state, sites) + + cutoff = 1E-10 + maxdim = 100 + nsweeps = 5 + + nsites = 2 + trunc = (; cutoff, maxdim) + E, gs_psi = dmrg(H, psi0; inserter_kwargs=(; trunc), nsites, nsweeps, outputlevel) + (outputlevel >= 1) && println("2-site DMRG energy = ", E) + + inserter_kwargs=(; trunc) + nsites = 1 + tmax = 0.10 + time_range = 0.0:0.02:tmax + psi1_t = time_evolve(H, gs_psi, time_range; inserter_kwargs, nsites, outputlevel) + (outputlevel >= 1) && println("Done with $nsites-site TDVP") + + @test norm(psi1_t) > 0.999 + + nsites = 2 + psi2_t = time_evolve(H, gs_psi, time_range; inserter_kwargs, nsites, outputlevel) + (outputlevel >= 1) && println("Done with $nsites-site TDVP") + @test norm(psi2_t) > 0.999 + + @test abs(inner(psi1_t, gs_psi)) > 0.99 + @test abs(inner(psi1_t, psi2_t)) > 0.99 + + # Test that accumulated phase angle is E*tmax + z = inner(psi1_t, gs_psi) + @test abs(atan(imag(z)/real(z)) - E*tmax) < 1E-4 +end From 3f7442f55b3c9831838051412b4b32c875c8cf4f Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 2 Jul 2025 10:50:21 -0400 Subject: [PATCH 09/28] Change default outputlevel in tests --- test/solvers/test_tree_dmrg.jl | 2 +- test/solvers/test_tree_tdvp.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/solvers/test_tree_dmrg.jl b/test/solvers/test_tree_dmrg.jl index 57459056..c5c8ab3d 100644 --- a/test/solvers/test_tree_dmrg.jl +++ b/test/solvers/test_tree_dmrg.jl @@ -9,7 +9,7 @@ include("utilities/simple_ed_methods.jl") include("utilities/tree_graphs.jl") @testset "Tree DMRG" begin - outputlevel = 1 + outputlevel = 0 g = build_tree(; nbranch=3, nbranch_sites=3) diff --git a/test/solvers/test_tree_tdvp.jl b/test/solvers/test_tree_tdvp.jl index 04bdece4..ee6731fe 100644 --- a/test/solvers/test_tree_tdvp.jl +++ b/test/solvers/test_tree_tdvp.jl @@ -21,7 +21,7 @@ function chain_plus_ancilla(; nchain) end @testset "Tree TDVP on chain plus ancilla" begin - outputlevel = 1 + outputlevel = 0 N = 10 g = chain_plus_ancilla(; nchain=N) From 8ef353f8ccb009c9d5e83b2f2c8456544d891295 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 2 Jul 2025 15:22:04 -0400 Subject: [PATCH 10/28] Rename and reorder some arguments --- src/ITensorNetworks.jl | 1 + src/apply.jl | 2 +- src/solvers/applyexp.jl | 23 ++++++++++------------- src/solvers/eigsolve.jl | 4 ++-- test/solvers/test_tree_tdvp.jl | 4 ++-- 5 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 582a51e2..ad209f33 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -64,6 +64,7 @@ include("solvers/permute_indices.jl") include("solvers/operator_map.jl") include("solvers/eigsolve.jl") include("solvers/applyexp.jl") +include("solvers/fitting.jl") include("apply.jl") include("inner.jl") diff --git a/src/apply.jl b/src/apply.jl index 6be06756..7da4bd33 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -186,7 +186,7 @@ function simple_update_bp(o, ψ, v⃗; envs, callback=Returns(nothing), apply_kw end function ITensors.apply( - o, + o::Union{ITensor,NamedEdge}, ψ::AbstractITensorNetwork; envs=ITensor[], normalize=false, diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 01c042d7..8b4b42a8 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -60,8 +60,9 @@ function applyexp_sweep_printer( end function applyexp( - init_prob, - exponents; + operator, + exponents, + state; extracter_kwargs=(;), updater_kwargs=(;), inserter_kwargs=(;), @@ -71,6 +72,9 @@ function applyexp( sweep_printer=applyexp_sweep_printer, kws..., ) + init_prob = ApplyExpProblem(; + state=permute_indices(init_state), operator=ProjTTN(permute_indices(operator)) + ) time_steps = diff([0.0, exponents...])[2:end] sweep_kws = (; outputlevel, extracter_kwargs, inserter_kwargs, nsites, tdvp_order, updater_kwargs @@ -81,23 +85,16 @@ function applyexp( return state(converged_prob) end -function applyexp(H, init_state, exponents; kws...) - init_prob = ApplyExpProblem(; - state=permute_indices(init_state), operator=ProjTTN(permute_indices(H)) - ) - return applyexp(init_prob, exponents; kws...) -end - process_real_times(z) = round(-imag(z); digits=10) function time_evolve( - H, - init_state, - time_points; + operator, + time_points, + init_state; process_time=process_real_times, sweep_printer=(a...; k...)->applyexp_sweep_printer(a...; process_time, k...), kws..., ) exponents = [-im*t for t in time_points] - return applyexp(H, init_state, exponents; sweep_printer, kws...) + return applyexp(operator, exponents, init_state; sweep_printer, kws...) end diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 0edc26dd..9e5305d4 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -66,9 +66,9 @@ function eigsolve( return eigenvalue(prob), state(prob) end -function eigsolve(H, init_state; kws...) +function eigsolve(operator, init_state; kws...) init_prob = EigsolveProblem(; - state=permute_indices(init_state), operator=ProjTTN(permute_indices(H)) + state=permute_indices(init_state), operator=ProjTTN(permute_indices(operator)) ) return eigsolve(init_prob; kws...) end diff --git a/test/solvers/test_tree_tdvp.jl b/test/solvers/test_tree_tdvp.jl index ee6731fe..676f49f3 100644 --- a/test/solvers/test_tree_tdvp.jl +++ b/test/solvers/test_tree_tdvp.jl @@ -57,13 +57,13 @@ end nsites = 1 tmax = 0.10 time_range = 0.0:0.02:tmax - psi1_t = time_evolve(H, gs_psi, time_range; inserter_kwargs, nsites, outputlevel) + psi1_t = time_evolve(H, time_range, gs_psi; inserter_kwargs, nsites, outputlevel) (outputlevel >= 1) && println("Done with $nsites-site TDVP") @test norm(psi1_t) > 0.999 nsites = 2 - psi2_t = time_evolve(H, gs_psi, time_range; inserter_kwargs, nsites, outputlevel) + psi2_t = time_evolve(H, time_range, gs_psi; inserter_kwargs, nsites, outputlevel) (outputlevel >= 1) && println("Done with $nsites-site TDVP") @test norm(psi2_t) > 0.999 From 7391c0c0fbf147a1364ecb9d204156e030ab9a62 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 2 Jul 2025 15:24:22 -0400 Subject: [PATCH 11/28] Rename test files --- test/solvers/{test_tree_tdvp.jl => test_tree_applyexp.jl} | 0 test/solvers/{test_tree_dmrg.jl => test_tree_eigsolve.jl} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename test/solvers/{test_tree_tdvp.jl => test_tree_applyexp.jl} (100%) rename test/solvers/{test_tree_dmrg.jl => test_tree_eigsolve.jl} (100%) diff --git a/test/solvers/test_tree_tdvp.jl b/test/solvers/test_tree_applyexp.jl similarity index 100% rename from test/solvers/test_tree_tdvp.jl rename to test/solvers/test_tree_applyexp.jl diff --git a/test/solvers/test_tree_dmrg.jl b/test/solvers/test_tree_eigsolve.jl similarity index 100% rename from test/solvers/test_tree_dmrg.jl rename to test/solvers/test_tree_eigsolve.jl From 735fa7e44e8e7b4c0f657cd4cdd9c90c7f252ac7 Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 3 Jul 2025 13:43:50 -0400 Subject: [PATCH 12/28] Bump version number to 0.14.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 12c8898f..edffe416 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" authors = ["Matthew Fishman , Joseph Tindall and contributors"] -version = "0.13.12" +version = "0.14.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" From 1d1dbc09ca81ced3be3ba23e3c5a1fcd9d4c9105 Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 4 Jul 2025 13:19:06 -0400 Subject: [PATCH 13/28] Rename permute_indices to align_indices --- src/solvers/applyexp.jl | 6 +++--- src/solvers/eigsolve.jl | 13 +++++-------- src/solvers/permute_indices.jl | 2 +- 3 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 8b4b42a8..e2511da2 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -60,9 +60,9 @@ function applyexp_sweep_printer( end function applyexp( - operator, + operator::AbstractITensorNetwork, exponents, - state; + init_state::AbstractITensorNetwork; extracter_kwargs=(;), updater_kwargs=(;), inserter_kwargs=(;), @@ -73,7 +73,7 @@ function applyexp( kws..., ) init_prob = ApplyExpProblem(; - state=permute_indices(init_state), operator=ProjTTN(permute_indices(operator)) + state=align_indices(init_state), operator=ProjTTN(align_indices(operator)) ) time_steps = diff([0.0, exponents...])[2:end] sweep_kws = (; diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 9e5305d4..a62cfe0c 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -43,7 +43,8 @@ function eigsolve_sweep_printer(region_iterator; outputlevel, sweep, nsweeps, kw end function eigsolve( - init_prob; + operator, + init_state; nsweeps, nsites=1, outputlevel=0, @@ -53,6 +54,9 @@ function eigsolve( sweep_printer=eigsolve_sweep_printer, kws..., ) + init_prob = EigsolveProblem(; + state=align_indices(init_state), operator=ProjTTN(align_indices(operator)) + ) sweep_iter = sweep_iterator( init_prob, nsweeps; @@ -66,11 +70,4 @@ function eigsolve( return eigenvalue(prob), state(prob) end -function eigsolve(operator, init_state; kws...) - init_prob = EigsolveProblem(; - state=permute_indices(init_state), operator=ProjTTN(permute_indices(operator)) - ) - return eigsolve(init_prob; kws...) -end - dmrg(args...; kws...) = eigsolve(args...; kws...) diff --git a/src/solvers/permute_indices.jl b/src/solvers/permute_indices.jl index 90526965..9a8e33cf 100644 --- a/src/solvers/permute_indices.jl +++ b/src/solvers/permute_indices.jl @@ -1,5 +1,5 @@ -function permute_indices(tn) +function align_indices(tn) si = siteinds(tn) ptn = copy(tn) for v in vertices(tn) From b918731697850bf97d29593d8e1ca26aa0b71d6a Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 4 Jul 2025 13:19:41 -0400 Subject: [PATCH 14/28] Rename file --- src/ITensorNetworks.jl | 2 +- src/solvers/{permute_indices.jl => align_indices.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/solvers/{permute_indices.jl => align_indices.jl} (100%) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index ad209f33..0c80c198 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -60,7 +60,7 @@ include("solvers/extracter.jl") include("solvers/inserter.jl") include("solvers/subspace/subspace.jl") include("solvers/subspace/densitymatrix.jl") -include("solvers/permute_indices.jl") +include("solvers/align_indices.jl") include("solvers/operator_map.jl") include("solvers/eigsolve.jl") include("solvers/applyexp.jl") diff --git a/src/solvers/permute_indices.jl b/src/solvers/align_indices.jl similarity index 100% rename from src/solvers/permute_indices.jl rename to src/solvers/align_indices.jl From 8928162411ec73800560135c222bd7e6804508df Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 4 Jul 2025 13:30:23 -0400 Subject: [PATCH 15/28] Implement many style changes --- src/solvers/applyexp.jl | 28 +++++++++++++++--------- src/solvers/eigsolve.jl | 29 ++++++++++++++++--------- src/solvers/extracter.jl | 9 ++------ src/solvers/inserter.jl | 6 ++--- src/solvers/region_plans.jl | 4 ++-- src/solvers/region_plans/dfs_plans.jl | 4 ++-- src/solvers/region_plans/euler_plans.jl | 4 ++-- src/solvers/region_plans/euler_tour.jl | 10 ++++----- src/solvers/subspace/densitymatrix.jl | 5 +---- src/solvers/subspace/subspace.jl | 1 - 10 files changed, 53 insertions(+), 47 deletions(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index e2511da2..11225626 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -1,9 +1,8 @@ using Printf: @printf -import ConstructionBase: setproperties @kwdef mutable struct ApplyExpProblem{State} - state::State operator + state::State current_time::Number = 0.0 end @@ -11,12 +10,20 @@ ITensorNetworks.state(A::ApplyExpProblem) = A.state operator(A::ApplyExpProblem) = A.operator current_time(A::ApplyExpProblem) = A.current_time +function set_operator(A::ApplyExpProblem, operator) + ApplyExpProblem(operator, A.state, A.current_time) +end +set_state(A::ApplyExpProblem, state) = ApplyExpProblem(A.operator, state, A.current_time) +function set_current_time(A::ApplyExpProblem, current_time) + ApplyExpProblem(A.operator, A.state, current_time) +end + function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...) return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) end function updater( - A::ApplyExpProblem, + prob::ApplyExpProblem, local_state, region_iterator; nsites, @@ -25,26 +32,27 @@ function updater( outputlevel, kws..., ) - local_state, info = solver(x->optimal_map(operator(A), x), time_step, local_state; kws...) + local_state, info = solver( + x->optimal_map(operator(prob), x), time_step, local_state; kws... + ) if nsites==1 curr_reg = current_region(region_iterator) next_reg = next_region(region_iterator) if !isnothing(next_reg) && next_reg != curr_reg - next_edge = first(edge_sequence_between_regions(state(A), curr_reg, next_reg)) + next_edge = first(edge_sequence_between_regions(state(prob), curr_reg, next_reg)) v1, v2 = src(next_edge), dst(next_edge) - psi = copy(state(A)) + psi = copy(state(prob)) psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2])) - shifted_operator = position(operator(A), psi, NamedEdge(v1=>v2)) + shifted_operator = position(operator(prob), psi, NamedEdge(v1=>v2)) R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...) local_state = psi[v1]*R_t end end - curr_time = current_time(A) + time_step - A = setproperties(A; current_time=curr_time) + prob = set_current_time(prob, current_time(prob)+time_step) - return A, local_state + return prob, local_state end function applyexp_sweep_printer( diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index a62cfe0c..0d773ef1 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,9 +1,8 @@ using Printf: @printf -import ConstructionBase: setproperties @kwdef mutable struct EigsolveProblem{State,Operator} - state::State operator::Operator + state::State eigenvalue::Number = Inf end @@ -11,20 +10,30 @@ eigenvalue(E::EigsolveProblem) = E.eigenvalue ITensorNetworks.state(E::EigsolveProblem) = E.state operator(E::EigsolveProblem) = E.operator +function set_operator(E::EigsolveProblem, operator) + EigsolveProblem(operator, E.state, E.eigenvalue) +end +function set_eigenvalue(E::EigsolveProblem, eigenvalue) + EigsolveProblem(E.operator, E.state, eigenvalue) +end +set_state(E::EigsolveProblem, state) = EigsolveProblem(E.operator, state, E.eigenvalue) + function updater( - E::EigsolveProblem, + prob::EigsolveProblem, local_state, region_iterator; outputlevel, solver=eigsolve_solver, kws..., ) - eigval, local_state = solver(ψ->optimal_map(operator(E), ψ), local_state; kws...) - E = setproperties(E; eigenvalue=eigval) + eigval, local_state = solver(ψ->optimal_map(operator(prob), ψ), local_state; kws...) + prob = set_eigenvalue(prob, eigval) if outputlevel >= 2 - @printf(" Region %s: energy = %.12f\n", current_region(region_iterator), eigenvalue(E)) + @printf( + " Region %s: energy = %.12f\n", current_region(region_iterator), eigenvalue(prob) + ) end - return E, local_state + return prob, local_state end function eigsolve_sweep_printer(region_iterator; outputlevel, sweep, nsweeps, kws...) @@ -34,9 +43,9 @@ function eigsolve_sweep_printer(region_iterator; outputlevel, sweep, nsweeps, kw else @printf("After sweep %d/%d ", sweep, nsweeps) end - E = problem(region_iterator) - @printf("eigenvalue=%.12f ", eigenvalue(E)) - @printf("maxlinkdim=%d", maxlinkdim(state(E))) + prob = problem(region_iterator) + @printf("eigenvalue=%.12f ", eigenvalue(prob)) + @printf("maxlinkdim=%d", maxlinkdim(state(prob))) println() flush(stdout) end diff --git a/src/solvers/extracter.jl b/src/solvers/extracter.jl index 1bf9f127..abce03c8 100644 --- a/src/solvers/extracter.jl +++ b/src/solvers/extracter.jl @@ -1,17 +1,12 @@ -import ConstructionBase: setproperties - function extracter(problem, region_iterator; sweep, trunc=(;), kws...) trunc = truncation_parameters(sweep; trunc...) region = current_region(region_iterator) psi = orthogonalize(state(problem), region) local_state = prod(psi[v] for v in region) - problem = setproperties(problem; state=psi) - + problem = set_state(problem, psi) problem, local_state = subspace_expand( problem, local_state, region_iterator; sweep, trunc, kws... ) - shifted_operator = position(operator(problem), state(problem), region) - - return setproperties(problem; operator=shifted_operator), local_state + return set_operator(problem, shifted_operator), local_state end diff --git a/src/solvers/inserter.jl b/src/solvers/inserter.jl index 2180fb14..60274ff2 100644 --- a/src/solvers/inserter.jl +++ b/src/solvers/inserter.jl @@ -1,5 +1,4 @@ -import ConstructionBase: setproperties -import NamedGraphs: edgetype +using NamedGraphs: edgetype function inserter( problem, @@ -12,7 +11,6 @@ function inserter( kws..., ) trunc = truncation_parameters(sweep; trunc...) - region = current_region(region_iterator) psi = copy(state(problem)) if length(region) == 1 @@ -30,5 +28,5 @@ function inserter( @preserve_graph psi[v] = C psi = set_orthogonal_region ? set_ortho_region(psi, [v]) : psi normalize && @preserve_graph psi[v] = psi[v] / norm(psi[v]) - return setproperties(problem; state=psi) + return set_state(problem, psi) end diff --git a/src/solvers/region_plans.jl b/src/solvers/region_plans.jl index be858dfd..53396aed 100644 --- a/src/solvers/region_plans.jl +++ b/src/solvers/region_plans.jl @@ -1,5 +1,5 @@ -import Graphs: AbstractGraph, AbstractEdge, edges, dst, src, vertices -import NamedGraphs: GraphsExtensions +using Graphs: AbstractGraph, AbstractEdge, edges, dst, src, vertices +using NamedGraphs: GraphsExtensions #function basic_path_regions(g::AbstractGraph; sweep_kwargs...) # fwd_sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges(g)] diff --git a/src/solvers/region_plans/dfs_plans.jl b/src/solvers/region_plans/dfs_plans.jl index facbaca6..074fa94a 100644 --- a/src/solvers/region_plans/dfs_plans.jl +++ b/src/solvers/region_plans/dfs_plans.jl @@ -1,5 +1,5 @@ -import Graphs: dst, src -import NamedGraphs.GraphsExtensions: +using Graphs: dst, src +using NamedGraphs.GraphsExtensions: default_root_vertex, post_order_dfs_edges, post_order_dfs_vertices function post_order_dfs_plan( diff --git a/src/solvers/region_plans/euler_plans.jl b/src/solvers/region_plans/euler_plans.jl index f3b9a96e..cf661d0d 100644 --- a/src/solvers/region_plans/euler_plans.jl +++ b/src/solvers/region_plans/euler_plans.jl @@ -1,5 +1,5 @@ -import Graphs: dst, src -import NamedGraphs.GraphsExtensions: default_root_vertex +using Graphs: dst, src +using NamedGraphs.GraphsExtensions: default_root_vertex function euler_sweep(graph; nsites, root_vertex=default_root_vertex(graph), sweep_kwargs...) if nsites == 1 diff --git a/src/solvers/region_plans/euler_tour.jl b/src/solvers/region_plans/euler_tour.jl index 6aeb0029..e5010d89 100644 --- a/src/solvers/region_plans/euler_tour.jl +++ b/src/solvers/region_plans/euler_tour.jl @@ -1,8 +1,8 @@ -import Graphs: dst, edges, src, vertices -import NamedGraphs as ng +using Graphs: dst, edges, src, vertices +using NamedGraphs: edgetype, vertextype function compute_adjacencies(G) - adj = Dict(v => Vector{ng.vertextype(G)}() for v in vertices(G)) + adj = Dict(v => Vector{vertextype(G)}() for v in vertices(G)) for e in edges(G) push!(adj[src(e)], dst(e)) push!(adj[dst(e)], src(e)) @@ -12,8 +12,8 @@ end function euler_tour_edges(G, start_vertex) adj = compute_adjacencies(G) - etype = ng.edgetype(G) - vtype = ng.vertextype(G) + etype = edgetype(G) + vtype = vertextype(G) visited = Set{Tuple{vtype,vtype}}() tour = Vector{etype}() stack = [start_vertex] diff --git a/src/solvers/subspace/densitymatrix.jl b/src/solvers/subspace/densitymatrix.jl index 7511901d..f70ad9eb 100644 --- a/src/solvers/subspace/densitymatrix.jl +++ b/src/solvers/subspace/densitymatrix.jl @@ -61,13 +61,10 @@ function subspace_expand( end Ax, ax = directsum(A=>a, U=>commonind(U, D)) - #println("Old space: ", space(a)) - #println("New space: ", space(ax)) - #ITensors.pause() expander = dag(Ax) * A psi[prev_vertex] = Ax psi[next_vertex] = expander * C local_state = expander*local_state - return setproperties(problem; state=psi), local_state + return set_state(problem, psi), local_state end diff --git a/src/solvers/subspace/subspace.jl b/src/solvers/subspace/subspace.jl index 7a332c3b..03102558 100644 --- a/src/solvers/subspace/subspace.jl +++ b/src/solvers/subspace/subspace.jl @@ -1,6 +1,5 @@ using NDTensors: NDTensors using NDTensors.BackendSelection: Backend, @Backend_str -import ConstructionBase: setproperties default_expansion_factor() = 1.5 default_max_expand() = typemax(Int) From 0a4db7718827a3fdcd33a364d625f8f9b58db783 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 11:37:40 -0400 Subject: [PATCH 16/28] Improve time_steps calculation code Better handling of zero and vector creation Co-authored-by: Matt Fishman --- src/solvers/applyexp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 01c042d7..dfaf6e99 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -71,7 +71,7 @@ function applyexp( sweep_printer=applyexp_sweep_printer, kws..., ) - time_steps = diff([0.0, exponents...])[2:end] + time_steps = diff([zero(eltype(exponents)); exponents])[2:end] sweep_kws = (; outputlevel, extracter_kwargs, inserter_kwargs, nsites, tdvp_order, updater_kwargs ) From 9abc7320eecdeb764d33cc2fbb2a5f91ca2aa4d9 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 22:54:56 -0400 Subject: [PATCH 17/28] Change name of callbacks and use overloading --- src/solvers/applyexp.jl | 12 +++++------- src/solvers/eigsolve.jl | 10 ++++------ src/solvers/sweep_solve.jl | 12 ++++++------ 3 files changed, 15 insertions(+), 19 deletions(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 11225626..d3547fe8 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -55,13 +55,12 @@ function updater( return prob, local_state end -function applyexp_sweep_printer( - region_iterator; outputlevel, sweep, nsweeps, process_time=identity, kws... +function sweep_printer( + problem::ApplyExpProblem; outputlevel, sweep, nsweeps, process_time=identity, kws... ) if outputlevel >= 1 - T = problem(region_iterator) - @printf(" Current time = %s, ", process_time(current_time(T))) - @printf("maxlinkdim=%d", maxlinkdim(state(T))) + @printf(" Current time = %s, ", process_time(current_time(problem))) + @printf("maxlinkdim=%d", maxlinkdim(state(problem))) println() flush(stdout) end @@ -77,7 +76,6 @@ function applyexp( outputlevel=0, nsites=1, tdvp_order=4, - sweep_printer=applyexp_sweep_printer, kws..., ) init_prob = ApplyExpProblem(; @@ -89,7 +87,7 @@ function applyexp( ) kws_array = [(; sweep_kws..., time_step=t) for t in time_steps] sweep_iter = sweep_iterator(init_prob, kws_array) - converged_prob = sweep_solve(sweep_iter; outputlevel, sweep_printer, kws...) + converged_prob = sweep_solve(sweep_iter; outputlevel, kws...) return state(converged_prob) end diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 0d773ef1..a3876829 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -36,16 +36,15 @@ function updater( return prob, local_state end -function eigsolve_sweep_printer(region_iterator; outputlevel, sweep, nsweeps, kws...) +function sweep_printer(problem::EigsolveProblem; outputlevel, sweep, nsweeps, kws...) if outputlevel >= 1 if nsweeps >= 10 @printf("After sweep %02d/%d ", sweep, nsweeps) else @printf("After sweep %d/%d ", sweep, nsweeps) end - prob = problem(region_iterator) - @printf("eigenvalue=%.12f ", eigenvalue(prob)) - @printf("maxlinkdim=%d", maxlinkdim(state(prob))) + @printf("eigenvalue=%.12f ", eigenvalue(problem)) + @printf("maxlinkdim=%d", maxlinkdim(state(problem))) println() flush(stdout) end @@ -60,7 +59,6 @@ function eigsolve( extracter_kwargs=(;), updater_kwargs=(;), inserter_kwargs=(;), - sweep_printer=eigsolve_sweep_printer, kws..., ) init_prob = EigsolveProblem(; @@ -75,7 +73,7 @@ function eigsolve( updater_kwargs, inserter_kwargs, ) - prob = sweep_solve(sweep_iter; outputlevel, sweep_printer, kws...) + prob = sweep_solve(sweep_iter; outputlevel, kws...) return eigenvalue(prob), state(prob) end diff --git a/src/solvers/sweep_solve.jl b/src/solvers/sweep_solve.jl index f5c5862d..b5496e5e 100644 --- a/src/solvers/sweep_solve.jl +++ b/src/solvers/sweep_solve.jl @@ -1,9 +1,9 @@ -default_region_callback(problem; kws...) = nothing +region_callback(problem; kws...) = nothing -default_sweep_callback(problem; kws...) = nothing +sweep_callback(problem; kws...) = nothing -function default_sweep_printer(problem; outputlevel, sweep, nsweeps, kws...) +function sweep_printer(problem; outputlevel, sweep, nsweeps, kws...) if outputlevel >= 1 println("Done with sweep $sweep/$nsweeps") end @@ -12,9 +12,9 @@ end function sweep_solve( sweep_iterator; outputlevel=0, - region_callback=default_region_callback, - sweep_callback=default_sweep_callback, - sweep_printer=default_sweep_printer, + region_callback=region_callback, + sweep_callback=sweep_callback, + sweep_printer=sweep_printer, kwargs..., ) for (sweep, region_iter) in enumerate(sweep_iterator) From 7719cc24b3b9f1208e34b0757cc23d7e38586b1b Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 22:55:09 -0400 Subject: [PATCH 18/28] Improve import/using code style --- test/solvers/test_tree_applyexp.jl | 10 +++++----- test/solvers/test_tree_eigsolve.jl | 8 ++++---- test/solvers/utilities/simple_ed_methods.jl | 2 +- test/solvers/utilities/tree_graphs.jl | 14 +++++++------- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/test/solvers/test_tree_applyexp.jl b/test/solvers/test_tree_applyexp.jl index 676f49f3..4d3dde40 100644 --- a/test/solvers/test_tree_applyexp.jl +++ b/test/solvers/test_tree_applyexp.jl @@ -1,10 +1,10 @@ using Test: @test, @testset using ITensors -using TensorOperations # Needed to use contraction order finding -import ITensorNetworks: dmrg, maxlinkdim, siteinds, time_evolve, ttn -import Graphs: add_vertex!, add_edge!, vertices -import NamedGraphs: NamedGraph -import ITensorMPS: OpSum +using ITensorNetworks: dmrg, maxlinkdim, siteinds, time_evolve, ttn +using Graphs: add_vertex!, add_edge!, vertices +using NamedGraphs: NamedGraph +using ITensorMPS: OpSum +using TensorOperations: TensorOperations #For contraction order finding function chain_plus_ancilla(; nchain) g = NamedGraph() diff --git a/test/solvers/test_tree_eigsolve.jl b/test/solvers/test_tree_eigsolve.jl index c5c8ab3d..1beecf3b 100644 --- a/test/solvers/test_tree_eigsolve.jl +++ b/test/solvers/test_tree_eigsolve.jl @@ -1,9 +1,9 @@ using Test: @test, @testset using ITensors -using TensorOperations # Needed to use contraction order finding using ITensorNetworks: siteinds, ttn, dmrg -import Graphs: dst, edges, src -import ITensorMPS: OpSum +using Graphs: dst, edges, src, vertices +using ITensorMPS: OpSum +using TensorOperations: TensorOperations #For contraction order finding include("utilities/simple_ed_methods.jl") include("utilities/tree_graphs.jl") @@ -27,7 +27,7 @@ include("utilities/tree_graphs.jl") # Make initial product state state = Dict{Tuple{Int,Int},String}() - for (j, v) in enumerate(gr.vertices(sites)) + for (j, v) in enumerate(vertices(sites)) state[v] = iseven(j) ? "Up" : "Dn" end psi0 = ttn(state, sites) diff --git a/test/solvers/utilities/simple_ed_methods.jl b/test/solvers/utilities/simple_ed_methods.jl index c80d3be8..e21afb9a 100644 --- a/test/solvers/utilities/simple_ed_methods.jl +++ b/test/solvers/utilities/simple_ed_methods.jl @@ -1,4 +1,4 @@ -import ITensorNetworks: AbstractITensorNetwork +using ITensorNetworks: AbstractITensorNetwork function ed_ground_state(H, psi0) ITensors.disable_warn_order() diff --git a/test/solvers/utilities/tree_graphs.jl b/test/solvers/utilities/tree_graphs.jl index 4fdff761..4021da63 100644 --- a/test/solvers/utilities/tree_graphs.jl +++ b/test/solvers/utilities/tree_graphs.jl @@ -1,5 +1,5 @@ -import Graphs as gr -import NamedGraphs as ng +using Graphs: add_edge!, add_vertex! +using NamedGraphs: NamedGraph """ build_tree @@ -8,15 +8,15 @@ import NamedGraphs as ng nbranch branches of nbranch_sites each. """ function build_tree(; nbranch=3, nbranch_sites=3) - g = ng.NamedGraph() - gr.add_vertex!(g, (0, 0)) + g = NamedGraph() + add_vertex!(g, (0, 0)) for branch in 1:nbranch, site in 1:nbranch_sites - gr.add_vertex!(g, (branch, site)) + add_vertex!(g, (branch, site)) end for branch in 1:nbranch - gr.add_edge!(g, (0, 0)=>(branch, 1)) + add_edge!(g, (0, 0)=>(branch, 1)) for site in 2:nbranch_sites - gr.add_edge!(g, (branch, site-1)=>(branch, site)) + add_edge!(g, (branch, site-1)=>(branch, site)) end end return g From 751902cfa17dca38fdf84cd99be7fb04a914aa92 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 22:55:25 -0400 Subject: [PATCH 19/28] Rename region_iterator_action to region_step --- src/solvers/iterators.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/iterators.jl b/src/solvers/iterators.jl index aaf530cf..413b994f 100644 --- a/src/solvers/iterators.jl +++ b/src/solvers/iterators.jl @@ -66,7 +66,7 @@ function Base.iterate(R::RegionIterator, which=1) region_plan_state = iterate(R.region_plan, which) isnothing(region_plan_state) && return nothing (current_region, region_kwargs), next = region_plan_state - R.problem = region_iterator_action(problem(R), R; region_kwargs...) + R.problem = region_step(problem(R), R; region_kwargs...) return R, next end @@ -78,7 +78,7 @@ function region_iterator(problem; sweep_kwargs...) return RegionIterator(; problem, region_plan=region_plan(problem; sweep_kwargs...)) end -function region_iterator_action( +function region_step( problem, region_iterator; extracter_kwargs=(;), From ddafd9f58c38d14771ccf831c5b93ae2f3162914 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 22:55:47 -0400 Subject: [PATCH 20/28] Change callbacks and printer to take problem --- src/solvers/sweep_solve.jl | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/solvers/sweep_solve.jl b/src/solvers/sweep_solve.jl index b5496e5e..ecad9730 100644 --- a/src/solvers/sweep_solve.jl +++ b/src/solvers/sweep_solve.jl @@ -23,6 +23,7 @@ function sweep_solve( problem(region_iter); nsweeps=length(sweep_iterator), outputlevel, + region_iterator=region_iter, region, region_kwargs, sweep, @@ -30,10 +31,20 @@ function sweep_solve( ) end sweep_callback( - region_iter; nsweeps=length(sweep_iterator), outputlevel, sweep, kwargs... + problem(region_iter); + nsweeps=length(sweep_iterator), + outputlevel, + region_iterator=region_iter, + sweep, + kwargs..., ) sweep_printer( - region_iter; nsweeps=length(sweep_iterator), outputlevel, sweep, kwargs... + problem(region_iter); + nsweeps=length(sweep_iterator), + outputlevel, + region_iterator=region_iter, + sweep, + kwargs..., ) end return problem(sweep_iterator) From 31a86d3f6aa264830c12b41292f4920c40079dfa Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 23:05:06 -0400 Subject: [PATCH 21/28] Remove older region_plans.jl file --- src/solvers/region_plans.jl | 73 ------------------------------------- 1 file changed, 73 deletions(-) delete mode 100644 src/solvers/region_plans.jl diff --git a/src/solvers/region_plans.jl b/src/solvers/region_plans.jl deleted file mode 100644 index 53396aed..00000000 --- a/src/solvers/region_plans.jl +++ /dev/null @@ -1,73 +0,0 @@ -using Graphs: AbstractGraph, AbstractEdge, edges, dst, src, vertices -using NamedGraphs: GraphsExtensions - -#function basic_path_regions(g::AbstractGraph; sweep_kwargs...) -# fwd_sweep = [([src(e), dst(e)], sweep_kwargs) for e in edges(g)] -# return [fwd_sweep..., reverse(fwd_sweep)...] -#end - -function tdvp_regions( - g::AbstractGraph, time_step; nsites=1, updater_kwargs, sweep_kwargs... -) - @assert nsites==1 - fwd_up_args = (; time=(time_step / 2), updater_kwargs...) - rev_up_args = (; time=(-time_step / 2), updater_kwargs...) - - fwd_sweep = [] - for e in edges(g) - push!(fwd_sweep, ([src(e)], (; updater_kwargs=fwd_up_args, sweep_kwargs...))) - push!(fwd_sweep, (e, (; updater_kwargs=rev_up_args, sweep_kwargs...))) - end - push!(fwd_sweep, ([dst(last(edges(g)))], (; updater_kwargs=fwd_up_args, sweep_kwargs...))) - - # Reverse regions as well as ordering of regions: - rev_sweep = [(reverse(rk[1]), rk[2]) for rk in reverse(fwd_sweep)] - - return [fwd_sweep..., rev_sweep...] -end - -function overlap(ea::AbstractEdge, eb::AbstractEdge) - return intersect([src(ea), dst(ea)], [src(eb), dst(eb)]) -end - -function forward_region(edges, which_edge; nsites=1, region_kwargs=(;)) - current_edge = edges[which_edge] - if nsites == 1 - #handle edge case - if current_edge == last(edges) - overlapping_vertex = only( - union([overlap(e, current_edge) for e in edges[1:(which_edge - 1)]]...) - ) - nonoverlapping_vertex = only( - setdiff([src(current_edge), dst(current_edge)], [overlapping_vertex]) - ) - return [ - ([overlapping_vertex], region_kwargs), ([nonoverlapping_vertex], region_kwargs) - ] - else - future_edges = edges[(which_edge + 1):end] - future_edges = isa(future_edges, AbstractEdge) ? [future_edges] : future_edges - overlapping_vertex = only(union([overlap(e, current_edge) for e in future_edges]...)) - nonoverlapping_vertex = only( - setdiff([src(current_edge), dst(current_edge)], [overlapping_vertex]) - ) - return [([nonoverlapping_vertex], region_kwargs)] - end - elseif nsites == 2 - return [([src(current_edge), dst(current_edge)], region_kwargs)] - end -end - -function basic_region_plan( - graph::AbstractGraph; - nsites, - root_vertex=GraphsExtensions.default_root_vertex(graph), - sweep_kwargs..., -) - edges = GraphsExtensions.post_order_dfs_edges(graph, root_vertex) - fwd_sweep = [ - forward_region(edges, i; nsites, region_kwargs=sweep_kwargs) for i in 1:length(edges) - ] - fwd_sweep = collect(Iterators.flatten(fwd_sweep)) - return [fwd_sweep..., reverse(fwd_sweep)...] -end From 3882692746de9cc37f0b372e99dc74e2185104b6 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 23:08:13 -0400 Subject: [PATCH 22/28] Change extracter etc to extract, update, insert. Change related kwarg names. --- src/solvers/applyexp.jl | 10 +++++----- src/solvers/eigsolve.jl | 16 +++++----------- src/solvers/extracter.jl | 2 +- src/solvers/inserter.jl | 2 +- src/solvers/iterators.jl | 18 +++++++----------- src/solvers/region_plans/tdvp_region_plans.jl | 14 +++++++------- src/solvers/subspace/subspace.jl | 2 +- test/solvers/test_tree_applyexp.jl | 8 ++++---- test/solvers/test_tree_eigsolve.jl | 10 +++++----- 9 files changed, 36 insertions(+), 46 deletions(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index d3547fe8..5fca915c 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -22,7 +22,7 @@ function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...) return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) end -function updater( +function update( prob::ApplyExpProblem, local_state, region_iterator; @@ -70,9 +70,9 @@ function applyexp( operator::AbstractITensorNetwork, exponents, init_state::AbstractITensorNetwork; - extracter_kwargs=(;), - updater_kwargs=(;), - inserter_kwargs=(;), + extract_kwargs=(;), + update_kwargs=(;), + insert_kwargs=(;), outputlevel=0, nsites=1, tdvp_order=4, @@ -83,7 +83,7 @@ function applyexp( ) time_steps = diff([0.0, exponents...])[2:end] sweep_kws = (; - outputlevel, extracter_kwargs, inserter_kwargs, nsites, tdvp_order, updater_kwargs + outputlevel, extract_kwargs, insert_kwargs, nsites, tdvp_order, update_kwargs ) kws_array = [(; sweep_kws..., time_step=t) for t in time_steps] sweep_iter = sweep_iterator(init_prob, kws_array) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index a3876829..438f7229 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -18,7 +18,7 @@ function set_eigenvalue(E::EigsolveProblem, eigenvalue) end set_state(E::EigsolveProblem, state) = EigsolveProblem(E.operator, state, E.eigenvalue) -function updater( +function update( prob::EigsolveProblem, local_state, region_iterator; @@ -56,22 +56,16 @@ function eigsolve( nsweeps, nsites=1, outputlevel=0, - extracter_kwargs=(;), - updater_kwargs=(;), - inserter_kwargs=(;), + extract_kwargs=(;), + update_kwargs=(;), + insert_kwargs=(;), kws..., ) init_prob = EigsolveProblem(; state=align_indices(init_state), operator=ProjTTN(align_indices(operator)) ) sweep_iter = sweep_iterator( - init_prob, - nsweeps; - nsites, - outputlevel, - extracter_kwargs, - updater_kwargs, - inserter_kwargs, + init_prob, nsweeps; nsites, outputlevel, extract_kwargs, update_kwargs, insert_kwargs ) prob = sweep_solve(sweep_iter; outputlevel, kws...) return eigenvalue(prob), state(prob) diff --git a/src/solvers/extracter.jl b/src/solvers/extracter.jl index abce03c8..011058af 100644 --- a/src/solvers/extracter.jl +++ b/src/solvers/extracter.jl @@ -1,4 +1,4 @@ -function extracter(problem, region_iterator; sweep, trunc=(;), kws...) +function extract(problem, region_iterator; sweep, trunc=(;), kws...) trunc = truncation_parameters(sweep; trunc...) region = current_region(region_iterator) psi = orthogonalize(state(problem), region) diff --git a/src/solvers/inserter.jl b/src/solvers/inserter.jl index 60274ff2..02ae8d73 100644 --- a/src/solvers/inserter.jl +++ b/src/solvers/inserter.jl @@ -1,6 +1,6 @@ using NamedGraphs: edgetype -function inserter( +function insert( problem, local_tensor, region_iterator; diff --git a/src/solvers/iterators.jl b/src/solvers/iterators.jl index 413b994f..023818ba 100644 --- a/src/solvers/iterators.jl +++ b/src/solvers/iterators.jl @@ -81,21 +81,17 @@ end function region_step( problem, region_iterator; - extracter_kwargs=(;), - updater_kwargs=(;), - inserter_kwargs=(;), + extract_kwargs=(;), + update_kwargs=(;), + insert_kwargs=(;), sweep, kws..., ) - problem, local_state = extracter( - problem, region_iterator; extracter_kwargs..., sweep, kws... - ) - problem, local_state = updater( - problem, local_state, region_iterator; updater_kwargs..., kws... - ) - problem = inserter( - problem, local_state, region_iterator; sweep, inserter_kwargs..., kws... + problem, local_state = extract(problem, region_iterator; extract_kwargs..., sweep, kws...) + problem, local_state = update( + problem, local_state, region_iterator; update_kwargs..., kws... ) + problem = insert(problem, local_state, region_iterator; sweep, insert_kwargs..., kws...) return problem end diff --git a/src/solvers/region_plans/tdvp_region_plans.jl b/src/solvers/region_plans/tdvp_region_plans.jl index 0624d666..cfe31555 100644 --- a/src/solvers/region_plans/tdvp_region_plans.jl +++ b/src/solvers/region_plans/tdvp_region_plans.jl @@ -12,16 +12,16 @@ function tdvp_sub_time_steps(tdvp_order) end function first_order_sweep( - graph, time_step, dir=Base.Forward; updater_kwargs, nsites, kws... + graph, time_step, dir=Base.Forward; update_kwargs, nsites, kws... ) basic_fwd_sweep = post_order_dfs_plan(graph; nsites, kws...) - updater_kwargs = (; nsites, time_step, updater_kwargs...) + update_kwargs = (; nsites, time_step, update_kwargs...) sweep = [] for (j, (region, region_kws)) in enumerate(basic_fwd_sweep) - push!(sweep, (region, (; nsites, updater_kwargs, region_kws...))) + push!(sweep, (region, (; nsites, update_kwargs, region_kws...))) if length(region) == 2 && j < length(basic_fwd_sweep) - rev_kwargs = (; updater_kwargs..., time_step=(-updater_kwargs.time_step)) - push!(sweep, ([last(region)], (; updater_kwargs=rev_kwargs, region_kws...))) + rev_kwargs = (; update_kwargs..., time_step=(-update_kwargs.time_step)) + push!(sweep, ([last(region)], (; update_kwargs=rev_kwargs, region_kws...))) end end if dir==Base.Reverse @@ -31,13 +31,13 @@ function first_order_sweep( return sweep end -function tdvp_regions(graph, time_step; updater_kwargs, tdvp_order, nsites, kws...) +function tdvp_regions(graph, time_step; update_kwargs, tdvp_order, nsites, kws...) sweep_plan = [] for (step, weight) in enumerate(tdvp_sub_time_steps(tdvp_order)) dir = isodd(step) ? Base.Forward : Base.Reverse append!( sweep_plan, - first_order_sweep(graph, weight*time_step, dir; updater_kwargs, nsites, kws...), + first_order_sweep(graph, weight*time_step, dir; update_kwargs, nsites, kws...), ) end return sweep_plan diff --git a/src/solvers/subspace/subspace.jl b/src/solvers/subspace/subspace.jl index 03102558..dad09300 100644 --- a/src/solvers/subspace/subspace.jl +++ b/src/solvers/subspace/subspace.jl @@ -39,7 +39,7 @@ function compute_expansion( # Restrict expand_maxdim below theoretical upper limit expand_maxdim = min(basis_size-current_dim, expand_maxdim) - # Enforce total maxdim setting (e.g. used in inserter step) + # Enforce total maxdim setting (e.g. used in insert step) expand_maxdim = min(maxdim-current_dim, expand_maxdim) # Ensure expand_maxdim is non-negative expand_maxdim = max(0, expand_maxdim) diff --git a/test/solvers/test_tree_applyexp.jl b/test/solvers/test_tree_applyexp.jl index 4d3dde40..5e8229bd 100644 --- a/test/solvers/test_tree_applyexp.jl +++ b/test/solvers/test_tree_applyexp.jl @@ -50,20 +50,20 @@ end nsites = 2 trunc = (; cutoff, maxdim) - E, gs_psi = dmrg(H, psi0; inserter_kwargs=(; trunc), nsites, nsweeps, outputlevel) + E, gs_psi = dmrg(H, psi0; insert_kwargs=(; trunc), nsites, nsweeps, outputlevel) (outputlevel >= 1) && println("2-site DMRG energy = ", E) - inserter_kwargs=(; trunc) + insert_kwargs=(; trunc) nsites = 1 tmax = 0.10 time_range = 0.0:0.02:tmax - psi1_t = time_evolve(H, time_range, gs_psi; inserter_kwargs, nsites, outputlevel) + psi1_t = time_evolve(H, time_range, gs_psi; insert_kwargs, nsites, outputlevel) (outputlevel >= 1) && println("Done with $nsites-site TDVP") @test norm(psi1_t) > 0.999 nsites = 2 - psi2_t = time_evolve(H, time_range, gs_psi; inserter_kwargs, nsites, outputlevel) + psi2_t = time_evolve(H, time_range, gs_psi; insert_kwargs, nsites, outputlevel) (outputlevel >= 1) && println("Done with $nsites-site TDVP") @test norm(psi2_t) > 0.999 diff --git a/test/solvers/test_tree_eigsolve.jl b/test/solvers/test_tree_eigsolve.jl index 1beecf3b..75194cae 100644 --- a/test/solvers/test_tree_eigsolve.jl +++ b/test/solvers/test_tree_eigsolve.jl @@ -45,8 +45,8 @@ include("utilities/tree_graphs.jl") # nsites = 2 trunc = (; cutoff, maxdim) - inserter_kwargs = (; trunc) - E, psi = dmrg(H, psi0; inserter_kwargs, nsites, nsweeps, outputlevel) + insert_kwargs = (; trunc) + E, psi = dmrg(H, psi0; insert_kwargs, nsites, nsweeps, outputlevel) (outputlevel >= 1) && println("2-site DMRG energy = ", E) @test abs(E-Ex) < 1E-5 @@ -56,9 +56,9 @@ include("utilities/tree_graphs.jl") nsites = 1 nsweeps = 5 trunc = (; cutoff, maxdim) - extracter_kwargs = (; trunc, subspace_algorithm="densitymatrix") - inserter_kwargs = (; trunc) - E, psi = dmrg(H, psi0; extracter_kwargs, inserter_kwargs, nsites, nsweeps, outputlevel) + extract_kwargs = (; trunc, subspace_algorithm="densitymatrix") + insert_kwargs = (; trunc) + E, psi = dmrg(H, psi0; extract_kwargs, insert_kwargs, nsites, nsweeps, outputlevel) (outputlevel >= 1) && println("1-site+subspace DMRG energy = ", E) @test abs(E-Ex) < 1E-5 end From ad874a5ac20f296f5ced516e21a3a6b3559e1e1d Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 5 Jul 2025 23:15:44 -0400 Subject: [PATCH 23/28] Fix line of code in time_evolve --- src/solvers/applyexp.jl | 2 +- test/solvers/test_tree_applyexp.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 5fca915c..36123c8b 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -98,7 +98,7 @@ function time_evolve( time_points, init_state; process_time=process_real_times, - sweep_printer=(a...; k...)->applyexp_sweep_printer(a...; process_time, k...), + sweep_printer=(a...; k...)->sweep_printer(a...; process_time, k...), kws..., ) exponents = [-im*t for t in time_points] diff --git a/test/solvers/test_tree_applyexp.jl b/test/solvers/test_tree_applyexp.jl index 5e8229bd..431d5a6f 100644 --- a/test/solvers/test_tree_applyexp.jl +++ b/test/solvers/test_tree_applyexp.jl @@ -20,7 +20,7 @@ function chain_plus_ancilla(; nchain) return g end -@testset "Tree TDVP on chain plus ancilla" begin +@testset "Test Applyexp" begin outputlevel = 0 N = 10 From ba9c1ef97ba927631140e9fe12442c80f8f43616 Mon Sep 17 00:00:00 2001 From: Miles Date: Mon, 7 Jul 2025 16:19:07 -0400 Subject: [PATCH 24/28] Commit fitting code and test --- src/solvers/fitting.jl | 113 +++++++++++++++++++++++++++++++++++ test/solvers/test_fitting.jl | 53 ++++++++++++++++ 2 files changed, 166 insertions(+) create mode 100644 src/solvers/fitting.jl create mode 100644 test/solvers/test_fitting.jl diff --git a/src/solvers/fitting.jl b/src/solvers/fitting.jl new file mode 100644 index 00000000..4cfc02e2 --- /dev/null +++ b/src/solvers/fitting.jl @@ -0,0 +1,113 @@ +using Graphs: vertices +using NamedGraphs: AbstractNamedGraph, NamedEdge +using NamedGraphs.PartitionedGraphs: partitionedges +using Printf: @printf +using ConstructionBase: setproperties + +@kwdef mutable struct FittingProblem{State<:AbstractBeliefPropagationCache} + state::State + ket_graph::AbstractNamedGraph + overlap::Number = 0 + gauge_region +end + +overlap(F::FittingProblem) = F.overlap +ITensorNetworks.state(F::FittingProblem) = F.state +ket_graph(F::FittingProblem) = F.ket_graph +gauge_region(F::FittingProblem) = F.gauge_region + +function set_state(F::FittingProblem, state) + FittingProblem(state, F.ket_graph, F.overlap, F.gauge_region) +end +function set_overlap(F::FittingProblem, overlap) + FittingProblem(F.state, F.ket_graph, overlap, F.gauge_region) +end + +function ket(F::FittingProblem) + ket_vertices = vertices(ket_graph(F)) + return first(induced_subgraph(tensornetwork(state(F)), ket_vertices)) +end + +function extract(problem::FittingProblem, region_iterator; sweep, kws...) + region = current_region(region_iterator) + prev_region = gauge_region(problem) + tn = state(problem) + path = edge_sequence_between_regions(ket_graph(problem), prev_region, region) + tn = gauge_walk(Algorithm("orthogonalize"), tn, path) + pe_path = partitionedges(partitioned_tensornetwork(tn), path) + tn = update( + Algorithm("bp"), tn, pe_path; message_update_function_kwargs=(; normalize=false) + ) + local_tensor = environment(tn, region) + sequence = contraction_sequence(local_tensor; alg="optimal") + local_tensor = dag(contract(local_tensor; sequence)) + #problem, local_tensor = subspace_expand(problem, local_tensor, region; sweep, kws...) + return setproperties(problem; state=tn, gauge_region=region), local_tensor +end + +function update(F::FittingProblem, local_tensor, region; outputlevel, kws...) + n = (local_tensor * dag(local_tensor))[] + F = set_overlap(F, n / sqrt(n)) + if outputlevel >= 2 + @printf(" Region %s: squared overlap = %.12f\n", region, overlap(F)) + end + return F, local_tensor +end + +function region_plan(F::FittingProblem; nsites, sweep_kwargs...) + return euler_sweep(ket_graph(F); nsites, sweep_kwargs...) +end + +function fit_tensornetwork( + overlap_network, + args...; + nsweeps=25, + nsites=1, + outputlevel=0, + extract_kwargs=(;), + update_kwargs=(;), + insert_kwargs=(;), + normalize=true, + kws..., +) + bpc = BeliefPropagationCache(overlap_network, args...) + ket_graph = first( + induced_subgraph(underlying_graph(overlap_network), ket_vertices(overlap_network)) + ) + init_prob = FittingProblem(; + ket_graph, state=bpc, gauge_region=collect(vertices(ket_graph)) + ) + + insert_kwargs = (; insert_kwargs..., normalize, set_orthogonal_region=false) + common_sweep_kwargs = (; nsites, outputlevel, update_kwargs, insert_kwargs) + kwargs_array = [(; common_sweep_kwargs..., sweep=s) for s in 1:nsweeps] + sweep_iter = sweep_iterator(init_prob, kwargs_array) + converged_prob = sweep_solve(sweep_iter; outputlevel, kws...) + return rename_vertices(inv_vertex_map(overlap_network), ket(converged_prob)) +end + +function fit_tensornetwork(tn, init_state, args...; kwargs...) + return fit_tensornetwork(inner_network(tn, init_state), args; kwargs...) +end + +#function truncate(tn; maxdim=default_maxdim(), cutoff=default_cutoff(), kwargs...) +# init_state = ITensorNetwork( +# v -> inds -> delta(inds), siteinds(tn); link_space=maxdim +# ) +# overlap_network = inner_network(tn, init_state) +# insert_kwargs = (; trunc=(; cutoff, maxdim)) +# return fit_tensornetwork(overlap_network; insert_kwargs, kwargs...) +#end + +function ITensors.apply( + A::ITensorNetwork, + x::ITensorNetwork; + maxdim=default_maxdim(), + cutoff=default_cutoff(), + kwargs..., +) + init_state = ITensorNetwork(v -> inds -> delta(inds), siteinds(x); link_space=maxdim) + overlap_network = inner_network(x, A, init_state) + insert_kwargs = (; trunc=(; cutoff, maxdim)) + return fit_tensornetwork(overlap_network; insert_kwargs, kwargs...) +end diff --git a/test/solvers/test_fitting.jl b/test/solvers/test_fitting.jl new file mode 100644 index 00000000..9395ae30 --- /dev/null +++ b/test/solvers/test_fitting.jl @@ -0,0 +1,53 @@ +using ITensors: apply, inner +using ITensorNetworks: ITensorNetwork, siteinds, ttn, random_tensornetwork +using ITensorNetworks.ModelHamiltonians: heisenberg +using NamedGraphs.NamedGraphGenerators: named_comb_tree +using Test: @test, @testset +using Printf +using StableRNGs: StableRNG +using TensorOperations: TensorOperations #For contraction order finding + +@testset "Fitting Tests" begin + outputlevel = 1 + for elt in (Float32, Float64, Complex{Float32}, Complex{Float64}) + (outputlevel >= 1) && println("\nFitting tests with elt = ", elt) + g = named_comb_tree((3, 2)) + s = siteinds("S=1/2", g) + + rng = StableRNG(1234) + + ##One-site truncation + #a = random_tensornetwork(rng, elt, s; link_space=3) + #b = truncate(a; maxdim=3) + #f = + # inner(a, b; alg="exact") / + # sqrt(inner(a, a; alg="exact") * inner(b, b; alg="exact")) + #(outputlevel >= 1) && @printf("One-site truncation. Fidelity = %s\n", f) + #@test abs(abs(f) - 1.0) <= 10*eps(real(elt)) + + ##Two-site truncation + #a = random_tensornetwork(rng, elt, s; link_space=3) + #b = truncate(a; maxdim=3, cutoff=1e-16, nsites=2) + #f = + # inner(a, b; alg="exact") / + # sqrt(inner(a, a; alg="exact") * inner(b, b; alg="exact")) + #(outputlevel >= 1) && @printf("Two-site truncation. Fidelity = %s\n", f) + #@test abs(abs(f) - 1.0) <= 10*eps(real(elt)) + + # #One-site apply (no normalization) + a = random_tensornetwork(rng, elt, s; link_space=2) + H = ITensorNetwork(ttn(heisenberg(g), s)) + Ha = apply(H, a; maxdim=4, nsites=1, normalize=false) + f = inner(Ha, a; alg="exact") / inner(a, H, a; alg="exact") + (outputlevel >= 1) && @printf("One-site apply. Fidelity = %s\n", f) + @test abs(f - 1.0) <= 500*eps(real(elt)) + + # #Two-site apply (no normalization) + a = random_tensornetwork(rng, elt, s; link_space=2) + H = ITensorNetwork(ttn(heisenberg(g), s)) + Ha = apply(H, a; maxdim=4, cutoff=1e-16, nsites=2, normalize=false) + f = inner(Ha, a; alg="exact") / inner(a, H, a; alg="exact") + (outputlevel >= 1) && @printf("Two-site apply. Fidelity = %s\n", f) + @test abs(f - 1.0) <= 500*eps(real(elt)) + end +end From 33d037a3bf65b63255fd11c674b8e67a6f44881c Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 8 Jul 2025 14:01:23 -0400 Subject: [PATCH 25/28] Update ITensorNetworks compat in Project.toml files --- docs/Project.toml | 2 +- examples/Project.toml | 2 +- test/Project.toml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 55b577dd..8fbc4729 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,5 +5,5 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" [compat] Documenter = "1.10.0" -ITensorNetworks = "0.13.0" +ITensorNetworks = "0.14.0" Literate = "2.20.1" diff --git a/examples/Project.toml b/examples/Project.toml index f8109c92..a3538ac6 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -2,4 +2,4 @@ ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" [compat] -ITensorNetworks = "0.13.2" +ITensorNetworks = "0.14.0" diff --git a/test/Project.toml b/test/Project.toml index 8f074e6a..372104b8 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -44,7 +44,7 @@ Glob = "1.3.1" Graphs = "1.12.0" GraphsFlows = "0.1.1" ITensorMPS = "0.3.6" -ITensorNetworks = "0.13.0" +ITensorNetworks = "0.14.0" ITensors = "0.7, 0.8, 0.9" KrylovKit = "0.8, 0.9" LinearAlgebra = "1.10.0" From 80cb6964afc4c24ed8b14fd0bb5f33e5b50732b5 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 9 Jul 2025 11:19:07 -0400 Subject: [PATCH 26/28] Rename file --- test/solvers/{test_tree_applyexp.jl => test_applyexp.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/solvers/{test_tree_applyexp.jl => test_applyexp.jl} (100%) diff --git a/test/solvers/test_tree_applyexp.jl b/test/solvers/test_applyexp.jl similarity index 100% rename from test/solvers/test_tree_applyexp.jl rename to test/solvers/test_applyexp.jl From 558f326c8211833d396313ffc05259d30f64146a Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 9 Jul 2025 11:57:06 -0400 Subject: [PATCH 27/28] Fix error in determining exponents. Change many variable names to more generic and correct exponent terminology. --- src/solvers/applyexp.jl | 59 ++++++++++++------- src/solvers/region_plans/tdvp_region_plans.jl | 22 +++---- test/solvers/test_applyexp.jl | 50 +++++++++++++++- 3 files changed, 96 insertions(+), 35 deletions(-) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index fe0cf72a..8e0f4325 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -3,23 +3,29 @@ using Printf: @printf @kwdef mutable struct ApplyExpProblem{State} operator state::State - current_time::Number = 0.0 + current_exponent::Number = 0.0 end ITensorNetworks.state(A::ApplyExpProblem) = A.state operator(A::ApplyExpProblem) = A.operator -current_time(A::ApplyExpProblem) = A.current_time +current_exponent(A::ApplyExpProblem) = A.current_exponent +function current_time(A::ApplyExpProblem) + t = im*A.current_exponent + return iszero(imag(t)) ? real(t) : t +end function set_operator(A::ApplyExpProblem, operator) - ApplyExpProblem(operator, A.state, A.current_time) + ApplyExpProblem(operator, A.state, A.current_exponent) +end +function set_state(A::ApplyExpProblem, state) + ApplyExpProblem(A.operator, state, A.current_exponent) end -set_state(A::ApplyExpProblem, state) = ApplyExpProblem(A.operator, state, A.current_time) -function set_current_time(A::ApplyExpProblem, current_time) - ApplyExpProblem(A.operator, A.state, current_time) +function set_current_exponent(A::ApplyExpProblem, current_exponent) + ApplyExpProblem(A.operator, A.state, current_exponent) end function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...) - return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) + return applyexp_regions(state(tdvp), time_step; nsites, sweep_kwargs...) end function update( @@ -27,15 +33,16 @@ function update( local_state, region_iterator; nsites, - time_step, + exponent_step, solver=runge_kutta_solver, outputlevel, kws..., ) + iszero(abs(exponent_step)) && return prob, local_state + local_state, info = solver( - x->optimal_map(operator(prob), x), time_step, local_state; kws... + x->optimal_map(operator(prob), x), exponent_step, local_state; kws... ) - if nsites==1 curr_reg = current_region(region_iterator) next_reg = next_region(region_iterator) @@ -45,21 +52,29 @@ function update( psi = copy(state(prob)) psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2])) shifted_operator = position(operator(prob), psi, NamedEdge(v1=>v2)) - R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...) + R_t, _ = solver(x->optimal_map(shifted_operator, x), -exponent_step, R; kws...) local_state = psi[v1]*R_t end end - prob = set_current_time(prob, current_time(prob)+time_step) + prob = set_current_exponent(prob, current_exponent(prob)+exponent_step) return prob, local_state end function sweep_printer( - problem::ApplyExpProblem; outputlevel, sweep, nsweeps, process_time=identity, kws... + problem::ApplyExpProblem; + exponent_description="exponent", + outputlevel, + sweep, + nsweeps, + process_time=identity, + kws..., ) if outputlevel >= 1 - @printf(" Current time = %s, ", process_time(current_time(problem))) + @printf( + " Current %s = %s, ", exponent_description, process_time(current_exponent(problem)) + ) @printf("maxlinkdim=%d", maxlinkdim(state(problem))) println() flush(stdout) @@ -75,30 +90,30 @@ function applyexp( insert_kwargs=(;), outputlevel=0, nsites=1, - tdvp_order=4, + order=4, kws..., ) init_prob = ApplyExpProblem(; state=align_indices(init_state), operator=ProjTTN(align_indices(operator)) ) - time_steps = diff([zero(eltype(exponents)); exponents])[2:end] - sweep_kws = (; - outputlevel, extract_kwargs, insert_kwargs, nsites, tdvp_order, update_kwargs - ) - kws_array = [(; sweep_kws..., time_step=t) for t in time_steps] + exponent_steps = diff([zero(eltype(exponents)); exponents]) + sweep_kws = (; outputlevel, extract_kwargs, insert_kwargs, nsites, order, update_kwargs) + kws_array = [(; sweep_kws..., time_step=t) for t in exponent_steps] sweep_iter = sweep_iterator(init_prob, kws_array) converged_prob = sweep_solve(sweep_iter; outputlevel, kws...) return state(converged_prob) end -process_real_times(z) = round(-imag(z); digits=10) +process_real_times(z) = iszero(abs(z)) ? 0.0 : round(-imag(z); digits=10) function time_evolve( operator, time_points, init_state; process_time=process_real_times, - sweep_printer=(a...; k...)->sweep_printer(a...; process_time, k...), + sweep_printer=( + a...; k... + )->sweep_printer(a...; exponent_description="time", process_time, k...), kws..., ) exponents = [-im*t for t in time_points] diff --git a/src/solvers/region_plans/tdvp_region_plans.jl b/src/solvers/region_plans/tdvp_region_plans.jl index cfe31555..c03ad4eb 100644 --- a/src/solvers/region_plans/tdvp_region_plans.jl +++ b/src/solvers/region_plans/tdvp_region_plans.jl @@ -1,26 +1,26 @@ -function tdvp_sub_time_steps(tdvp_order) - if tdvp_order == 1 +function applyexp_sub_steps(order) + if order == 1 return [1.0] - elseif tdvp_order == 2 + elseif order == 2 return [1 / 2, 1 / 2] - elseif tdvp_order == 4 + elseif order == 4 s = (2 - 2^(1 / 3))^(-1) return [s/2, s/2, 1/2 - s, 1/2 - s, s/2, s/2] else - error("TDVP order of $tdvp_order not supported") + error("Applyexp order of $order not supported") end end function first_order_sweep( - graph, time_step, dir=Base.Forward; update_kwargs, nsites, kws... + graph, exponent_step, dir=Base.Forward; update_kwargs, nsites, kws... ) basic_fwd_sweep = post_order_dfs_plan(graph; nsites, kws...) - update_kwargs = (; nsites, time_step, update_kwargs...) + update_kwargs = (; nsites, exponent_step, update_kwargs...) sweep = [] for (j, (region, region_kws)) in enumerate(basic_fwd_sweep) push!(sweep, (region, (; nsites, update_kwargs, region_kws...))) if length(region) == 2 && j < length(basic_fwd_sweep) - rev_kwargs = (; update_kwargs..., time_step=(-update_kwargs.time_step)) + rev_kwargs = (; update_kwargs..., exponent_step=(-update_kwargs.exponent_step)) push!(sweep, ([last(region)], (; update_kwargs=rev_kwargs, region_kws...))) end end @@ -31,13 +31,13 @@ function first_order_sweep( return sweep end -function tdvp_regions(graph, time_step; update_kwargs, tdvp_order, nsites, kws...) +function applyexp_regions(graph, exponent_step; update_kwargs, order, nsites, kws...) sweep_plan = [] - for (step, weight) in enumerate(tdvp_sub_time_steps(tdvp_order)) + for (step, weight) in enumerate(applyexp_sub_steps(order)) dir = isodd(step) ? Base.Forward : Base.Reverse append!( sweep_plan, - first_order_sweep(graph, weight*time_step, dir; update_kwargs, nsites, kws...), + first_order_sweep(graph, weight*exponent_step, dir; update_kwargs, nsites, kws...), ) end return sweep_plan diff --git a/test/solvers/test_applyexp.jl b/test/solvers/test_applyexp.jl index 431d5a6f..44e250ae 100644 --- a/test/solvers/test_applyexp.jl +++ b/test/solvers/test_applyexp.jl @@ -1,8 +1,9 @@ using Test: @test, @testset using ITensors -using ITensorNetworks: dmrg, maxlinkdim, siteinds, time_evolve, ttn +using ITensorNetworks: ITensorNetworks, applyexp, dmrg, maxlinkdim, siteinds, time_evolve, ttn using Graphs: add_vertex!, add_edge!, vertices using NamedGraphs: NamedGraph +using NamedGraphs.NamedGraphGenerators: named_path_graph using ITensorMPS: OpSum using TensorOperations: TensorOperations #For contraction order finding @@ -20,7 +21,7 @@ function chain_plus_ancilla(; nchain) return g end -@testset "Test Applyexp" begin +@testset "Test Tree Time Evolution" begin outputlevel = 0 N = 10 @@ -74,3 +75,48 @@ end z = inner(psi1_t, gs_psi) @test abs(atan(imag(z)/real(z)) - E*tmax) < 1E-4 end + +@testset "Applyexp Time Point Handling" begin + N = 10 + g = named_path_graph(N) + sites = siteinds("S=1/2", g) + + # Make Heisenberg model Hamiltonian + h = OpSum() + for j in 1:(N - 1) + h += "Sz", j, "Sz", j+1 + h += 1/2, "S+", j, "S-", j+1 + h += 1/2, "S-", j, "S+", j+1 + end + H = ttn(h, sites) + + # Initial product state + state = Dict{Int,String}() + for (j, v) in enumerate(vertices(sites)) + state[v] = iseven(j) ? "Up" : "Dn" + end + psi0 = ttn(state, sites) + + nsites = 2 + trunc = (; cutoff=1E-8, maxdim=100) + insert_kwargs=(; trunc) + + # Test that all time points are reached and reported correctly + time_points = [0.0,0.1,0.25,0.32,0.4] + times = Real[] + function collect_times(problem; kws...) + push!(times, ITensorNetworks.current_time(problem)) + end + time_evolve(H, time_points, psi0; insert_kwargs, nsites, sweep_callback=collect_times,outputlevel=1) + @test norm(times - time_points) < 10*eps(Float64) + + # Test that all exponents are reached and reported correctly + exponent_points = [-0.0,-0.1,-0.25,-0.32,-0.4] + exponents = Real[] + function collect_exponents(problem; kws...) + push!(exponents, ITensorNetworks.current_exponent(problem)) + end + applyexp(H, exponent_points, psi0; insert_kwargs, nsites, sweep_callback=collect_exponents,outputlevel=1) + @test norm(exponents - exponent_points) < 10*eps(Float64) +end + From fc5c4e83d654c2572e6538dae4fc7e5d127f1ed8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 10 Jul 2025 02:04:24 +0000 Subject: [PATCH 28/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- test/solvers/test_applyexp.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/solvers/test_applyexp.jl b/test/solvers/test_applyexp.jl index 44e250ae..69b83b7c 100644 --- a/test/solvers/test_applyexp.jl +++ b/test/solvers/test_applyexp.jl @@ -119,4 +119,3 @@ end applyexp(H, exponent_points, psi0; insert_kwargs, nsites, sweep_callback=collect_exponents,outputlevel=1) @test norm(exponents - exponent_points) < 10*eps(Float64) end -