Network solvers redesign #242
Merged
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
Fyi, At least one test was failing on my machine, but it didn't seem to be related at all to the code in this PR. |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #242 +/- ##
===========================================
- Coverage 77.04% 58.75% -18.29%
===========================================
Files 76 82 +6
Lines 3868 4044 +176
===========================================
- Hits 2980 2376 -604
- Misses 888 1668 +780
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
Is the version check meaning I should just bump the version number? Which version number should I bump to? |
The version check is checking whether or not you've bumped the package version here: ITensorNetworks.jl/Project.toml Line 4 in 2d7db49 |
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
src/solvers/eigsolve.jl
Outdated
| end | ||
|
|
||
| eigenvalue(E::EigsolveProblem) = E.eigenvalue | ||
| ITensorNetworks.state(E::EigsolveProblem) = E.state |
There was a problem hiding this comment.
Suggested change
| ITensorNetworks.state(E::EigsolveProblem) = E.state | |
| state(E::EigsolveProblem) = E.state |
There was a problem hiding this comment.
I was getting an error if I didn't prepend ITensorNetworks here, which is why I put it. I'll find out the exact error and maybe we can fix the issue more at the root of what's causing it.
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
src/solvers/adapters.jl
Outdated
| # generates each tuple? | ||
| # | ||
|
|
||
| mutable struct TupleRegionIterator{RegionIter} |
There was a problem hiding this comment.
Maybe we can come up with a name that is more descriptive, like RegionIteratorWithKwargs. Tuple is a bit vague.
mtfishman
reviewed
Jul 2, 2025
src/solvers/applyexp.jl
Outdated
| current_time::Number = 0.0 | ||
| end | ||
|
|
||
| ITensorNetworks.state(A::ApplyExpProblem) = A.state |
There was a problem hiding this comment.
Suggested change
| ITensorNetworks.state(A::ApplyExpProblem) = A.state | |
| state(A::ApplyExpProblem) = A.state |
since we are in the ITensorNetworks module.
There was a problem hiding this comment.
I'll look into what error I was getting when I didn't prepend ITensorNetworks. I agree it shouldn't be needed in principle.
mtfishman
reviewed
Jul 2, 2025
src/solvers/iterators.jl
Outdated
|
|
||
| if !isnothing(which) | ||
| S.region_iter = region_iterator( | ||
| problem(S.region_iter); sweep=S.which_sweep, current_sweep_kws... |
There was a problem hiding this comment.
Suggested change
| problem(S.region_iter); sweep=S.which_sweep, current_sweep_kws... | |
| problem(S); sweep=S.which_sweep, current_sweep_kws... |
I think, based on the definition of problem(::SweepIterator) above.
There was a problem hiding this comment.
Also, this code pattern confused me a bit, what do you think of writing it like this:
update_region_iterator!(S; current_sweep_kws...)and hiding the implementation in update_region_iterator!?
mtfishman
reviewed
Jul 2, 2025
mtfishman
reviewed
Jul 2, 2025
src/solvers/iterators.jl
Outdated
Comment on lines
5
to
9
| mutable struct SweepIterator | ||
| sweep_kws | ||
| region_iter | ||
| which_sweep::Int | ||
| end |
There was a problem hiding this comment.
What about parameterizing this struct by the types of the fields sweep_kws and region_iter? Is it meant to be dynamic, i.e. is it expected that those types might change?
There was a problem hiding this comment.
It seems best to have a goal of having those be concrete and not changing type, but if they do you can still parameterize the type and then set the type parameter to an abstract type when it is being constructed (as needed), for example you can explicitly construct it as SweepIterator{Any,Any}(sweep_kws, region_iter, which_sweep).
mtfishman
reviewed
Jul 2, 2025
src/solvers/iterators.jl
Outdated
|
|
||
| problem(R::RegionIterator) = R.problem | ||
| current_region_plan(R::RegionIterator) = R.region_plan[R.which_region] | ||
| current_region(R::RegionIterator) = current_region_plan(R)[1] |
There was a problem hiding this comment.
What's the data structure for the output of current_region_plan(R)? Something doesn't feel right to me that the region is accessed with indexing by 1, maybe it should be a NamedTuple and the region could be accessed as current_region_plan(R).region or a struct where it is accessed with a function get_region(current_region_plan(R))?
mtfishman
reviewed
Jul 2, 2025
src/solvers/iterators.jl
Outdated
| 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] |
There was a problem hiding this comment.
It seems like maybe we should be using function accessors rather than field accessors, i.e. which_region(R) instead of R.which_region.
mtfishman
reviewed
Jul 15, 2025
src/solvers/applyexp.jl
Outdated
| operator(A::ApplyExpProblem) = A.operator | ||
| current_exponent(A::ApplyExpProblem) = A.current_exponent | ||
| function current_time(A::ApplyExpProblem) | ||
| t = im*A.current_exponent |
There was a problem hiding this comment.
Suggested change
| t = im*A.current_exponent | |
| t = im * current_exponent(A) |
Just a style change.
mtfishman
reviewed
Jul 15, 2025
| current_exponent(A::ApplyExpProblem) = A.current_exponent | ||
| function current_time(A::ApplyExpProblem) | ||
| t = im*A.current_exponent | ||
| return iszero(imag(t)) ? real(t) : t |
There was a problem hiding this comment.
I'm not sure how I feel about introducing this type instability, let's discuss this. I would lean towards a design where this code is not so clever and just returns im * current_exponent(A), and we provide a different interface for users to be more explicit about the types.
There was a problem hiding this comment.
Different interfaces for each case is a good suggestion. How about current_time which returns just the real part of im*A.current_exponent (so usual "physics time") and also current_complex_time which returns a time (im*exponent) but as a complex value, in the sense of "complex time evolution" methods.
There was a problem hiding this comment.
That sounds good. Maybe current_time should check that iszero(imag(t)) and if not error, silently ignoring the imaginary part may be confusing.
mtfishman
reviewed
Jul 15, 2025
|
|
||
| function sweep_printer( | ||
| problem::ApplyExpProblem; | ||
| exponent_description="exponent", |
There was a problem hiding this comment.
This input is a bit odd to me, let's discuss this.
mtfishman
reviewed
Jul 15, 2025
| @@ -0,0 +1,16 @@ | |||
|
|
|||
| function align_indices(tn) | |||
There was a problem hiding this comment.
Suggested change
| function align_indices(tn) | |
| function align_inds(tn) |
since in general ITensor indices are referred to as inds in function names.
JoeyT1994
reviewed
Jul 16, 2025
There was a problem hiding this comment.
A similar change to the function function ITensors.apply( o, ψ::VidalITensorNetwork; ...) lower down (lime 322) will be necessary to avoid ambiguity and the current failing test.
Co-authored-by: Jack Dunham <jdunham@flatironinstitute.org> Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
|
Your PR requires formatting changes to meet the project's style guidelines. Click here to view the suggested changes.diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl
index 57af3ed..5643fce 100644
--- a/src/itensornetwork.jl
+++ b/src/itensornetwork.jl
@@ -180,7 +180,7 @@ function generic_state(x::Op, inds::NamedTuple)
end
function generic_state(s::AbstractString, inds::NamedTuple)
# TODO: Figure out what to do if there is more than one site.
- site_tensors = [ITensors.state(s, only(inds.siteinds))]
+ site_tensors = [ITensors.state(s, only(inds.siteinds))]
link_tensors = [[onehot(i => 1) for i in inds.linkinds[e]] for e in keys(inds.linkinds)]
return contract(reduce(vcat, link_tensors; init=site_tensors))
end
diff --git a/src/solvers/adapters.jl b/src/solvers/adapters.jl
index e1c58d1..a9ec55b 100644
--- a/src/solvers/adapters.jl
+++ b/src/solvers/adapters.jl
@@ -5,8 +5,8 @@ Iterator wrapper whos `compute!` function simply returns itself, doing nothing i
process. This allows one to manually call a custom `compute!` or insert their own code it in
the loop body in place of `compute!`.
"""
-struct IncrementOnly{S<:AbstractNetworkIterator} <: AbstractNetworkIterator
- parent::S
+struct IncrementOnly{S <: AbstractNetworkIterator} <: AbstractNetworkIterator
+ parent::S
end
islaststep(adapter::IncrementOnly) = islaststep(adapter.parent)
@@ -22,8 +22,8 @@ IncrementOnly(adapter::IncrementOnly) = adapter
Adapter that flattens each region iterator in the parent sweep iterator into a single
iterator.
"""
-struct EachRegion{SI<:SweepIterator} <: AbstractNetworkIterator
- parent::SI
+struct EachRegion{SI <: SweepIterator} <: AbstractNetworkIterator
+ parent::SI
end
# In keeping with Julia convention.
@@ -31,16 +31,16 @@ eachregion(iter::SweepIterator) = EachRegion(iter)
# Essential definitions
function islaststep(adapter::EachRegion)
- region_iter = region_iterator(adapter.parent)
- return islaststep(adapter.parent) && islaststep(region_iter)
+ region_iter = region_iterator(adapter.parent)
+ return islaststep(adapter.parent) && islaststep(region_iter)
end
function increment!(adapter::EachRegion)
- region_iter = region_iterator(adapter.parent)
- islaststep(region_iter) ? increment!(adapter.parent) : increment!(region_iter)
- return adapter
+ region_iter = region_iterator(adapter.parent)
+ islaststep(region_iter) ? increment!(adapter.parent) : increment!(region_iter)
+ return adapter
end
function compute!(adapter::EachRegion)
- region_iter = region_iterator(adapter.parent)
- compute!(region_iter)
- return adapter
+ region_iter = region_iterator(adapter.parent)
+ compute!(region_iter)
+ return adapter
end
diff --git a/src/solvers/align_indices.jl b/src/solvers/align_indices.jl
index 9a8e33c..8612114 100644
--- a/src/solvers/align_indices.jl
+++ b/src/solvers/align_indices.jl
@@ -1,16 +1,16 @@
function align_indices(tn)
- si = siteinds(tn)
- ptn = copy(tn)
- for v in 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])
+ si = siteinds(tn)
+ ptn = copy(tn)
+ for v in 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
- ptn[v] = permute(tn[v], new_is)
- end
- return ptn
+ return ptn
end
diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl
index 3b91ff8..fc78037 100644
--- a/src/solvers/applyexp.jl
+++ b/src/solvers/applyexp.jl
@@ -1,123 +1,123 @@
using Printf: @printf
@kwdef mutable struct ApplyExpProblem{State} <: AbstractProblem
- operator
- state::State
- current_exponent::Number = 0.0
+ operator
+ state::State
+ current_exponent::Number = 0.0
end
operator(A::ApplyExpProblem) = A.operator
state(A::ApplyExpProblem) = A.state
current_exponent(A::ApplyExpProblem) = A.current_exponent
function current_time(A::ApplyExpProblem)
- t = im * A.current_exponent
- return iszero(imag(t)) ? real(t) : t
+ t = im * A.current_exponent
+ return iszero(imag(t)) ? real(t) : t
end
# Rename region_plan
function region_plan(A::ApplyExpProblem; nsites, exponent_step, sweep_kwargs...)
- # The `exponent_step` kwarg for the `update!` function needs some pre-processing.
- return applyexp_regions(state(A), exponent_step; nsites, sweep_kwargs...)
+ # The `exponent_step` kwarg for the `update!` function needs some pre-processing.
+ return applyexp_regions(state(A), exponent_step; nsites, sweep_kwargs...)
end
function update!(
- region_iter::RegionIterator{<:ApplyExpProblem},
- local_state;
- nsites,
- exponent_step,
- solver=runge_kutta_solver,
-)
- prob = problem(region_iter)
-
- if iszero(abs(exponent_step))
- return region_iter, local_state
- end
-
- solver_kwargs = region_kwargs(solver, region_iter)
-
- local_state, _ = solver(
- x -> optimal_map(operator(prob), x), exponent_step, local_state; solver_kwargs...
- )
- if nsites == 1
- curr_reg = current_region(region_iter)
- next_reg = next_region(region_iter)
- if !isnothing(next_reg) && next_reg != curr_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(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), -exponent_step, R; solver_kwargs...
- )
- local_state = psi[v1] * R_t
+ region_iter::RegionIterator{<:ApplyExpProblem},
+ local_state;
+ nsites,
+ exponent_step,
+ solver = runge_kutta_solver,
+ )
+ prob = problem(region_iter)
+
+ if iszero(abs(exponent_step))
+ return region_iter, local_state
+ end
+
+ solver_kwargs = region_kwargs(solver, region_iter)
+
+ local_state, _ = solver(
+ x -> optimal_map(operator(prob), x), exponent_step, local_state; solver_kwargs...
+ )
+ if nsites == 1
+ curr_reg = current_region(region_iter)
+ next_reg = next_region(region_iter)
+ if !isnothing(next_reg) && next_reg != curr_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(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), -exponent_step, R; solver_kwargs...
+ )
+ local_state = psi[v1] * R_t
+ end
end
- end
- prob.current_exponent += exponent_step
+ prob.current_exponent += exponent_step
- return region_iter, local_state
+ return region_iter, local_state
end
function default_sweep_callback(
- sweep_iterator::SweepIterator{<:ApplyExpProblem};
- exponent_description="exponent",
- outputlevel=0,
- process_time=identity,
-)
- if outputlevel >= 1
- the_problem = problem(sweep_iterator)
- @printf(
- " Current %s = %s, ",
- exponent_description,
- process_time(current_exponent(the_problem))
+ sweep_iterator::SweepIterator{<:ApplyExpProblem};
+ exponent_description = "exponent",
+ outputlevel = 0,
+ process_time = identity,
)
- @printf("maxlinkdim=%d", maxlinkdim(state(the_problem)))
- println()
- flush(stdout)
- end
+ return if outputlevel >= 1
+ the_problem = problem(sweep_iterator)
+ @printf(
+ " Current %s = %s, ",
+ exponent_description,
+ process_time(current_exponent(the_problem))
+ )
+ @printf("maxlinkdim=%d", maxlinkdim(state(the_problem)))
+ println()
+ flush(stdout)
+ end
end
function applyexp(
- init_prob::AbstractProblem,
- exponents;
- sweep_callback=default_sweep_callback,
- order=4,
- nsites=2,
- sweep_kwargs...,
-)
- exponent_steps = diff([zero(eltype(exponents)); exponents])
-
- kws_array = [
- (; order, nsites, sweep_kwargs..., exponent_step) for exponent_step in exponent_steps
- ]
- sweep_iter = SweepIterator(init_prob, kws_array)
-
- converged_prob = problem(sweep_solve!(sweep_callback, sweep_iter))
-
- return state(converged_prob)
+ init_prob::AbstractProblem,
+ exponents;
+ sweep_callback = default_sweep_callback,
+ order = 4,
+ nsites = 2,
+ sweep_kwargs...,
+ )
+ exponent_steps = diff([zero(eltype(exponents)); exponents])
+
+ kws_array = [
+ (; order, nsites, sweep_kwargs..., exponent_step) for exponent_step in exponent_steps
+ ]
+ sweep_iter = SweepIterator(init_prob, kws_array)
+
+ converged_prob = problem(sweep_solve!(sweep_callback, sweep_iter))
+
+ return state(converged_prob)
end
function applyexp(operator, exponents, init_state; kws...)
- init_prob = ApplyExpProblem(;
- state=align_indices(init_state),
- operator=ProjTTN(align_indices(operator)),
- current_exponent=first(exponents),
- )
- return applyexp(init_prob, exponents; kws...)
+ init_prob = ApplyExpProblem(;
+ state = align_indices(init_state),
+ operator = ProjTTN(align_indices(operator)),
+ current_exponent = first(exponents),
+ )
+ return applyexp(init_prob, exponents; kws...)
end
-process_real_times(z) = iszero(abs(z)) ? 0.0 : 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_callback=iter ->
- default_sweep_callback(iter; exponent_description="time", process_time),
- sweep_kwargs...,
-)
- exponents = [-im * t for t in time_points]
- return applyexp(operator, exponents, init_state; sweep_callback, sweep_kwargs...)
+ operator,
+ time_points,
+ init_state;
+ process_time = process_real_times,
+ sweep_callback = iter ->
+ default_sweep_callback(iter; exponent_description = "time", process_time),
+ sweep_kwargs...,
+ )
+ exponents = [-im * t for t in time_points]
+ return applyexp(operator, exponents, init_state; sweep_callback, sweep_kwargs...)
end
diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl
index 907e71f..8fbed9b 100644
--- a/src/solvers/eigsolve.jl
+++ b/src/solvers/eigsolve.jl
@@ -1,11 +1,11 @@
using Printf: @printf
using ITensors: truncerror
-@kwdef mutable struct EigsolveProblem{State,Operator} <: AbstractProblem
- operator::Operator
- state::State
- eigenvalue::Number = Inf
- max_truncerror::Real = 0.0
+@kwdef mutable struct EigsolveProblem{State, Operator} <: AbstractProblem
+ operator::Operator
+ state::State
+ eigenvalue::Number = Inf
+ max_truncerror::Real = 0.0
end
eigenvalue(E::EigsolveProblem) = E.eigenvalue
@@ -13,69 +13,69 @@ state(E::EigsolveProblem) = E.state
operator(E::EigsolveProblem) = E.operator
max_truncerror(E::EigsolveProblem) = E.max_truncerror
-function set_truncation_info!(E::EigsolveProblem; spectrum=nothing)
- if !isnothing(spectrum)
- E.max_truncerror = max(max_truncerror(E), truncerror(spectrum))
- end
- return E
+function set_truncation_info!(E::EigsolveProblem; spectrum = nothing)
+ if !isnothing(spectrum)
+ E.max_truncerror = max(max_truncerror(E), truncerror(spectrum))
+ end
+ return E
end
function update!(
- region_iter::RegionIterator{<:EigsolveProblem},
- local_state;
- outputlevel=0,
- solver=eigsolve_solver,
-)
- prob = problem(region_iter)
+ region_iter::RegionIterator{<:EigsolveProblem},
+ local_state;
+ outputlevel = 0,
+ solver = eigsolve_solver,
+ )
+ prob = problem(region_iter)
- eigval, local_state = solver(
- ψ -> optimal_map(operator(prob), ψ), local_state; region_kwargs(solver, region_iter)...
- )
+ eigval, local_state = solver(
+ ψ -> optimal_map(operator(prob), ψ), local_state; region_kwargs(solver, region_iter)...
+ )
- prob.eigenvalue = eigval
+ prob.eigenvalue = eigval
- if outputlevel >= 2
- @printf(" Region %s: energy = %.12f\n", current_region(region_iter), eigenvalue(prob))
- end
- return region_iter, local_state
+ if outputlevel >= 2
+ @printf(" Region %s: energy = %.12f\n", current_region(region_iter), eigenvalue(prob))
+ end
+ return region_iter, local_state
end
function default_sweep_callback(
- sweep_iterator::SweepIterator{<:EigsolveProblem}; outputlevel=0
-)
- if outputlevel >= 1
- nsweeps = length(sweep_iterator)
- current_sweep = sweep_iterator.which_sweep
- if length(sweep_iterator) >= 10
- @printf("After sweep %02d/%d ", current_sweep, nsweeps)
- else
- @printf("After sweep %d/%d ", current_sweep, nsweeps)
+ sweep_iterator::SweepIterator{<:EigsolveProblem}; outputlevel = 0
+ )
+ return if outputlevel >= 1
+ nsweeps = length(sweep_iterator)
+ current_sweep = sweep_iterator.which_sweep
+ if length(sweep_iterator) >= 10
+ @printf("After sweep %02d/%d ", current_sweep, nsweeps)
+ else
+ @printf("After sweep %d/%d ", current_sweep, nsweeps)
+ end
+ @printf("eigenvalue=%.12f", eigenvalue(problem))
+ @printf(" maxlinkdim=%d", maxlinkdim(state(problem)))
+ @printf(" max truncerror=%d", max_truncerror(problem))
+ println()
+ flush(stdout)
end
- @printf("eigenvalue=%.12f", eigenvalue(problem))
- @printf(" maxlinkdim=%d", maxlinkdim(state(problem)))
- @printf(" max truncerror=%d", max_truncerror(problem))
- println()
- flush(stdout)
- end
end
function eigsolve(
- operator, init_state; nsweeps, nsites=1, outputlevel=0, factorize_kwargs, sweep_kwargs...
-)
- init_prob = EigsolveProblem(;
- state=align_indices(init_state), operator=ProjTTN(align_indices(operator))
- )
- sweep_iter = SweepIterator(
- init_prob,
- nsweeps;
- nsites,
- outputlevel,
- factorize_kwargs,
- subspace_expand!_kwargs=(; eigen_kwargs=factorize_kwargs),
- sweep_kwargs...,
- )
- prob = problem(sweep_solve!(sweep_iter))
- return eigenvalue(prob), state(prob)
+ operator, init_state; nsweeps, nsites = 1, outputlevel = 0, factorize_kwargs, sweep_kwargs...
+ )
+ init_prob = EigsolveProblem(;
+ state = align_indices(init_state), operator = ProjTTN(align_indices(operator))
+ )
+ sweep_iter = SweepIterator(
+ init_prob,
+ nsweeps;
+ nsites,
+ outputlevel,
+ factorize_kwargs,
+ subspace_expand!_kwargs = (; eigen_kwargs = factorize_kwargs),
+ sweep_kwargs...,
+ )
+ prob = problem(sweep_solve!(sweep_iter))
+ return eigenvalue(prob), state(prob)
end
dmrg(operator, init_state; kwargs...) = eigsolve(operator, init_state; kwargs...)
diff --git a/src/solvers/extract.jl b/src/solvers/extract.jl
index 629b70f..28ebddd 100644
--- a/src/solvers/extract.jl
+++ b/src/solvers/extract.jl
@@ -1,16 +1,16 @@
-function extract!(region_iter::RegionIterator; subspace_algorithm="nothing")
- prob = problem(region_iter)
- region = current_region(region_iter)
+function extract!(region_iter::RegionIterator; subspace_algorithm = "nothing")
+ prob = problem(region_iter)
+ region = current_region(region_iter)
- psi = orthogonalize(state(prob), region)
- local_state = prod(psi[v] for v in region)
+ psi = orthogonalize(state(prob), region)
+ local_state = prod(psi[v] for v in region)
- prob.state = psi
+ prob.state = psi
- _, local_state = subspace_expand!(region_iter, local_state; subspace_algorithm)
- shifted_operator = position(operator(prob), state(prob), region)
+ _, local_state = subspace_expand!(region_iter, local_state; subspace_algorithm)
+ shifted_operator = position(operator(prob), state(prob), region)
- prob.operator = shifted_operator
+ prob.operator = shifted_operator
- return region_iter, local_state
+ return region_iter, local_state
end
diff --git a/src/solvers/fitting.jl b/src/solvers/fitting.jl
index c03852e..ffd127e 100644
--- a/src/solvers/fitting.jl
+++ b/src/solvers/fitting.jl
@@ -3,12 +3,12 @@ using NamedGraphs: AbstractNamedGraph, NamedEdge
using NamedGraphs.PartitionedGraphs: partitionedges
using Printf: @printf
-@kwdef mutable struct FittingProblem{State<:AbstractBeliefPropagationCache} <:
- AbstractProblem
- state::State
- ket_graph::AbstractNamedGraph
- overlap::Number = 0
- gauge_region
+@kwdef mutable struct FittingProblem{State <: AbstractBeliefPropagationCache} <:
+ AbstractProblem
+ state::State
+ ket_graph::AbstractNamedGraph
+ overlap::Number = 0
+ gauge_region
end
state(F::FittingProblem) = F.state
@@ -17,86 +17,86 @@ overlap(F::FittingProblem) = F.overlap
gauge_region(F::FittingProblem) = F.gauge_region
function ket(F::FittingProblem)
- ket_vertices = vertices(ket_graph(F))
- return first(induced_subgraph(tensornetwork(state(F)), ket_vertices))
+ ket_vertices = vertices(ket_graph(F))
+ return first(induced_subgraph(tensornetwork(state(F)), ket_vertices))
end
function extract!(region_iter::RegionIterator{<:FittingProblem})
- prob = problem(region_iter)
-
- region = current_region(region_iter)
- prev_region = gauge_region(prob)
- tn = state(prob)
- path = edge_sequence_between_regions(ket_graph(prob), 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...)
-
- prob.state = tn
- prob.gauge_region = region
-
- return region_iter, local_tensor
+ prob = problem(region_iter)
+
+ region = current_region(region_iter)
+ prev_region = gauge_region(prob)
+ tn = state(prob)
+ path = edge_sequence_between_regions(ket_graph(prob), 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...)
+
+ prob.state = tn
+ prob.gauge_region = region
+
+ return region_iter, local_tensor
end
function update!(
- region_iter::RegionIterator{<:FittingProblem}, local_tensor; outputlevel=0
-)
- F = problem(region_iter)
+ region_iter::RegionIterator{<:FittingProblem}, local_tensor; outputlevel = 0
+ )
+ F = problem(region_iter)
- region = current_region(region_iter)
+ region = current_region(region_iter)
- n = (local_tensor * dag(local_tensor))[]
- F.overlap = n / sqrt(n)
+ n = (local_tensor * dag(local_tensor))[]
+ F.overlap = n / sqrt(n)
- if outputlevel >= 2
- @printf(" Region %s: squared overlap = %.12f\n", region, overlap(F))
- end
+ if outputlevel >= 2
+ @printf(" Region %s: squared overlap = %.12f\n", region, overlap(F))
+ end
- return region_iter, local_tensor
+ return region_iter, local_tensor
end
function region_plan(F::FittingProblem; nsites, sweep_kwargs...)
- return euler_sweep(ket_graph(F); 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,
- normalize=true,
- factorize_kwargs,
- extra_sweep_kwargs...,
-)
- 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 = (; normalize, set_orthogonal_region=false)
- update!_kwargs = (; outputlevel)
-
- sweep_kwargs = (; nsites, outputlevel, update!_kwargs, insert!_kwargs, factorize_kwargs)
- kwargs_array = [(; sweep_kwargs..., extra_sweep_kwargs..., sweep) for sweep in 1:nsweeps]
-
- sweep_iter = SweepIterator(init_prob, kwargs_array)
- converged_prob = problem(sweep_solve!(sweep_iter))
-
- return rename_vertices(inv_vertex_map(overlap_network), ket(converged_prob))
+ overlap_network,
+ args...;
+ nsweeps = 25,
+ nsites = 1,
+ outputlevel = 0,
+ normalize = true,
+ factorize_kwargs,
+ extra_sweep_kwargs...,
+ )
+ 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 = (; normalize, set_orthogonal_region = false)
+ update!_kwargs = (; outputlevel)
+
+ sweep_kwargs = (; nsites, outputlevel, update!_kwargs, insert!_kwargs, factorize_kwargs)
+ kwargs_array = [(; sweep_kwargs..., extra_sweep_kwargs..., sweep) for sweep in 1:nsweeps]
+
+ sweep_iter = SweepIterator(init_prob, kwargs_array)
+ converged_prob = problem(sweep_solve!(sweep_iter))
+
+ 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...)
+ return fit_tensornetwork(inner_network(tn, init_state), args; kwargs...)
end
#function truncate(tn; maxdim=default_maxdim(), cutoff=default_cutoff(), kwargs...)
@@ -109,15 +109,15 @@ end
#end
function ITensors.apply(
- A::AbstractITensorNetwork,
- x::AbstractITensorNetwork;
- maxdim=typemax(Int),
- cutoff=0.0,
- sweep_kwargs...,
-)
- init_state = ITensorNetwork(v -> inds -> delta(inds), siteinds(x); link_space=maxdim)
- overlap_network = inner_network(x, A, init_state)
- return fit_tensornetwork(
- overlap_network; factorize_kwargs=(; maxdim, cutoff), sweep_kwargs...
- )
+ A::AbstractITensorNetwork,
+ x::AbstractITensorNetwork;
+ maxdim = typemax(Int),
+ cutoff = 0.0,
+ sweep_kwargs...,
+ )
+ init_state = ITensorNetwork(v -> inds -> delta(inds), siteinds(x); link_space = maxdim)
+ overlap_network = inner_network(x, A, init_state)
+ return fit_tensornetwork(
+ overlap_network; factorize_kwargs = (; maxdim, cutoff), sweep_kwargs...
+ )
end
diff --git a/src/solvers/insert.jl b/src/solvers/insert.jl
index 87ffaf6..3896395 100644
--- a/src/solvers/insert.jl
+++ b/src/solvers/insert.jl
@@ -1,32 +1,32 @@
using NamedGraphs: edgetype
-function insert!(region_iter, local_tensor; normalize=false, set_orthogonal_region=true)
- prob = problem(region_iter)
+function insert!(region_iter, local_tensor; normalize = false, set_orthogonal_region = true)
+ prob = problem(region_iter)
- region = current_region(region_iter)
- psi = copy(state(prob))
- if length(region) == 1
- C = local_tensor
- elseif length(region) == 2
- e = edgetype(psi)(first(region), last(region))
- indsTe = inds(psi[first(region)])
- tags = ITensors.tags(psi, e)
+ region = current_region(region_iter)
+ psi = copy(state(prob))
+ if length(region) == 1
+ C = local_tensor
+ elseif length(region) == 2
+ e = edgetype(psi)(first(region), last(region))
+ indsTe = inds(psi[first(region)])
+ tags = ITensors.tags(psi, e)
- U, C, spectrum = factorize(
- local_tensor, indsTe; tags, region_kwargs(factorize, region_iter)...
- )
+ U, C, spectrum = factorize(
+ local_tensor, indsTe; tags, region_kwargs(factorize, region_iter)...
+ )
- @preserve_graph psi[first(region)] = U
- prob = set_truncation_info!(prob; spectrum)
- else
- error("Region of length $(length(region)) not currently supported")
- end
- v = last(region)
- @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])
+ @preserve_graph psi[first(region)] = U
+ prob = set_truncation_info!(prob; spectrum)
+ else
+ error("Region of length $(length(region)) not currently supported")
+ end
+ v = last(region)
+ @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])
- prob.state = psi
+ prob.state = psi
- return region_iter
+ return region_iter
end
diff --git a/src/solvers/iterators.jl b/src/solvers/iterators.jl
index 16497f0..990d86d 100644
--- a/src/solvers/iterators.jl
+++ b/src/solvers/iterators.jl
@@ -11,14 +11,14 @@ abstract type AbstractNetworkIterator end
# We use greater than or equals here as we increment the state at the start of the iteration
islaststep(iterator::AbstractNetworkIterator) = state(iterator) >= length(iterator)
-function Base.iterate(iterator::AbstractNetworkIterator, init=true)
- islaststep(iterator) && return nothing
- # We seperate increment! from step! and demand that any AbstractNetworkIterator *must*
- # define a method for increment! This way we avoid cases where one may wish to nest
- # calls to different step! methods accidentaly incrementing multiple times.
- init || increment!(iterator)
- rv = compute!(iterator)
- return rv, false
+function Base.iterate(iterator::AbstractNetworkIterator, init = true)
+ islaststep(iterator) && return nothing
+ # We seperate increment! from step! and demand that any AbstractNetworkIterator *must*
+ # define a method for increment! This way we avoid cases where one may wish to nest
+ # calls to different step! methods accidentaly incrementing multiple times.
+ init || increment!(iterator)
+ rv = compute!(iterator)
+ return rv, false
end
function increment! end
@@ -26,10 +26,10 @@ compute!(iterator::AbstractNetworkIterator) = iterator
step!(iterator::AbstractNetworkIterator) = step!(identity, iterator)
function step!(f, iterator::AbstractNetworkIterator)
- compute!(iterator)
- f(iterator)
- increment!(iterator)
- return iterator
+ compute!(iterator)
+ f(iterator)
+ increment!(iterator)
+ return iterator
end
#
@@ -38,23 +38,23 @@ end
"""
struct RegionIterator{Problem, RegionPlan} <: AbstractNetworkIterator
"""
-mutable struct RegionIterator{Problem,RegionPlan} <: AbstractNetworkIterator
- problem::Problem
- region_plan::RegionPlan
- which_region::Int
- const which_sweep::Int
- function RegionIterator(problem::P, region_plan::R, sweep::Int) where {P,R}
- return new{P,R}(problem, region_plan, 1, sweep)
- end
+mutable struct RegionIterator{Problem, RegionPlan} <: AbstractNetworkIterator
+ problem::Problem
+ region_plan::RegionPlan
+ which_region::Int
+ const which_sweep::Int
+ function RegionIterator(problem::P, region_plan::R, sweep::Int) where {P, R}
+ return new{P, R}(problem, region_plan, 1, sweep)
+ end
end
function RegionIterator(problem; sweep, sweep_kwargs...)
- plan = region_plan(problem; sweep_kwargs...)
- return RegionIterator(problem, plan, sweep)
+ plan = region_plan(problem; sweep_kwargs...)
+ return RegionIterator(problem, plan, sweep)
end
function new_region_iterator(iterator::RegionIterator; sweep_kwargs...)
- return RegionIterator(iterator.problem; sweep_kwargs...)
+ return RegionIterator(iterator.problem; sweep_kwargs...)
end
state(region_iter::RegionIterator) = region_iter.which_region
@@ -63,48 +63,48 @@ Base.length(region_iter::RegionIterator) = length(region_iter.region_plan)
problem(region_iter::RegionIterator) = region_iter.problem
function current_region_plan(region_iter::RegionIterator)
- return region_iter.region_plan[region_iter.which_region]
+ return region_iter.region_plan[region_iter.which_region]
end
function current_region(region_iter::RegionIterator)
- region, _ = current_region_plan(region_iter)
- return region
+ region, _ = current_region_plan(region_iter)
+ return region
end
function region_kwargs(region_iter::RegionIterator)
- _, kwargs = current_region_plan(region_iter)
- return kwargs
+ _, kwargs = current_region_plan(region_iter)
+ return kwargs
end
function region_kwargs(f::Function, iter::RegionIterator)
- return get(region_kwargs(iter), Symbol(f, :_kwargs), (;))
+ return get(region_kwargs(iter), Symbol(f, :_kwargs), (;))
end
function prev_region(region_iter::RegionIterator)
- state(region_iter) <= 1 && return nothing
- prev, _ = region_iter.region_plan[region_iter.which_region - 1]
- return prev
+ state(region_iter) <= 1 && return nothing
+ prev, _ = region_iter.region_plan[region_iter.which_region - 1]
+ return prev
end
function next_region(region_iter::RegionIterator)
- islaststep(region_iter) && return nothing
- next, _ = region_iter.region_plan[region_iter.which_region + 1]
- return next
+ islaststep(region_iter) && return nothing
+ next, _ = region_iter.region_plan[region_iter.which_region + 1]
+ return next
end
#
# Functions associated with RegionIterator
#
function increment!(region_iter::RegionIterator)
- region_iter.which_region += 1
- return region_iter
+ region_iter.which_region += 1
+ return region_iter
end
function compute!(iter::RegionIterator)
- _, local_state = extract!(iter; region_kwargs(extract!, iter)...)
- _, local_state = update!(iter, local_state; region_kwargs(update!, iter)...)
- insert!(iter, local_state; region_kwargs(insert!, iter)...)
+ _, local_state = extract!(iter; region_kwargs(extract!, iter)...)
+ _, local_state = update!(iter, local_state; region_kwargs(update!, iter)...)
+ insert!(iter, local_state; region_kwargs(insert!, iter)...)
- return iter
+ return iter
end
region_plan(problem; sweep_kwargs...) = euler_sweep(state(problem); sweep_kwargs...)
@@ -113,16 +113,16 @@ region_plan(problem; sweep_kwargs...) = euler_sweep(state(problem); sweep_kwargs
# SweepIterator
#
-mutable struct SweepIterator{Problem,Iter} <: AbstractNetworkIterator
- region_iter::RegionIterator{Problem}
- sweep_kwargs::Iterators.Stateful{Iter}
- which_sweep::Int
- function SweepIterator(problem::Prob, sweep_kwargs::Iter) where {Prob,Iter}
- stateful_sweep_kwargs = Iterators.Stateful(sweep_kwargs)
- first_kwargs, _ = Iterators.peel(stateful_sweep_kwargs)
- region_iter = RegionIterator(problem; sweep=1, first_kwargs...)
- return new{Prob,Iter}(region_iter, stateful_sweep_kwargs, 1)
- end
+mutable struct SweepIterator{Problem, Iter} <: AbstractNetworkIterator
+ region_iter::RegionIterator{Problem}
+ sweep_kwargs::Iterators.Stateful{Iter}
+ which_sweep::Int
+ function SweepIterator(problem::Prob, sweep_kwargs::Iter) where {Prob, Iter}
+ stateful_sweep_kwargs = Iterators.Stateful(sweep_kwargs)
+ first_kwargs, _ = Iterators.peel(stateful_sweep_kwargs)
+ region_iter = RegionIterator(problem; sweep = 1, first_kwargs...)
+ return new{Prob, Iter}(region_iter, stateful_sweep_kwargs, 1)
+ end
end
islaststep(sweep_iter::SweepIterator) = isnothing(peek(sweep_iter.sweep_kwargs))
@@ -133,27 +133,28 @@ problem(sweep_iter::SweepIterator) = problem(region_iterator(sweep_iter))
state(sweep_iter::SweepIterator) = sweep_iter.which_sweep
Base.length(sweep_iter::SweepIterator) = length(sweep_iter.sweep_kwargs)
function increment!(sweep_iter::SweepIterator)
- sweep_iter.which_sweep += 1
- sweep_kwargs, _ = Iterators.peel(sweep_iter.sweep_kwargs)
- update_region_iterator!(sweep_iter; sweep_kwargs...)
- return sweep_iter
+ sweep_iter.which_sweep += 1
+ sweep_kwargs, _ = Iterators.peel(sweep_iter.sweep_kwargs)
+ update_region_iterator!(sweep_iter; sweep_kwargs...)
+ return sweep_iter
end
function update_region_iterator!(iterator::SweepIterator; kwargs...)
- sweep = state(iterator)
- iterator.region_iter = new_region_iterator(iterator.region_iter; sweep, kwargs...)
- return iterator
+ sweep = state(iterator)
+ iterator.region_iter = new_region_iterator(iterator.region_iter; sweep, kwargs...)
+ return iterator
end
function compute!(sweep_iter::SweepIterator)
- for _ in sweep_iter.region_iter
- # TODO: Is it sensible to execute the default region callback function?
- end
+ for _ in sweep_iter.region_iter
+ # TODO: Is it sensible to execute the default region callback function?
+ end
+ return
end
# More basic constructor where sweep_kwargs are constant throughout sweeps
function SweepIterator(problem, nsweeps::Int; sweep_kwargs...)
- # Initialize this to an empty RegionIterator
- sweep_kwargs_iter = Iterators.repeated(sweep_kwargs, nsweeps)
- return SweepIterator(problem, sweep_kwargs_iter)
+ # Initialize this to an empty RegionIterator
+ sweep_kwargs_iter = Iterators.repeated(sweep_kwargs, nsweeps)
+ return SweepIterator(problem, sweep_kwargs_iter)
end
diff --git a/src/solvers/local_solvers/eigsolve.jl b/src/solvers/local_solvers/eigsolve.jl
index 96f1c9f..e86fd85 100644
--- a/src/solvers/local_solvers/eigsolve.jl
+++ b/src/solvers/local_solvers/eigsolve.jl
@@ -1,20 +1,20 @@
using KrylovKit: KrylovKit
function eigsolve_solver(
- operator,
- init,
- howmany=1;
+ operator,
+ init,
+ howmany = 1;
which_eigval=:SR,
ishermitian=true,
- tol=1E-14,
+ tol = 1.0e-14,
krylovdim=3,
maxiter=1,
verbosity=0,
eager=false,
- kws...,
+ kws...,
)
- vals, vecs, info = KrylovKit.eigsolve(
- operator,
+ vals, vecs, info = KrylovKit.eigsolve(
+ operator,
init,
howmany,
which_eigval;
@@ -25,5 +25,5 @@ function eigsolve_solver(
verbosity,
eager,
)
- return vals[1], vecs[1]
+ return vals[1], vecs[1]
end
diff --git a/src/solvers/local_solvers/exponentiate.jl b/src/solvers/local_solvers/exponentiate.jl
index b092991..d167f40 100644
--- a/src/solvers/local_solvers/exponentiate.jl
+++ b/src/solvers/local_solvers/exponentiate.jl
@@ -1,9 +1,9 @@
using KrylovKit: KrylovKit
function exponentiate_solver(
- operator,
- time,
- init;
+ operator,
+ time,
+ init;
krylovdim=30,
maxiter=100,
verbosity=0,
@@ -11,11 +11,11 @@ function exponentiate_solver(
ishermitian=true,
issymmetric=true,
eager=true,
- kws...,
+ kws...,
)
- result, exp_info = KrylovKit.exponentiate(
- operator,
- time,
+ result, exp_info = KrylovKit.exponentiate(
+ operator,
+ time,
init;
eager,
krylovdim,
@@ -25,5 +25,5 @@ function exponentiate_solver(
ishermitian,
issymmetric,
)
- return result, exp_info
+ return result, exp_info
end
diff --git a/src/solvers/local_solvers/runge_kutta.jl b/src/solvers/local_solvers/runge_kutta.jl
index 12ecd2e..2998bcf 100644
--- a/src/solvers/local_solvers/runge_kutta.jl
+++ b/src/solvers/local_solvers/runge_kutta.jl
@@ -1,12 +1,12 @@
-function runge_kutta_solver(H, t, ψ; order=4, kws...)
- # For linear ODE, Runge-Kutta is a Taylor series.
- # Pattern below derived as:
- # exp(tH)ψ = ψ + tHψ + (tH)^2(ψ)/2! + (tH)^3(ψ)/3! + ...
- # = ψ + (tH) * (ψ + (tH)/2 * (ψ + (tH)/3 * (ψ + ...)))
- eHψ = copy(ψ)
- for ord in reverse(1:order)
- eHψ = (t/ord)*H(eHψ) + ψ
- end
- return eHψ, (;)
+function runge_kutta_solver(H, t, ψ; order = 4, kws...)
+ # For linear ODE, Runge-Kutta is a Taylor series.
+ # Pattern below derived as:
+ # exp(tH)ψ = ψ + tHψ + (tH)^2(ψ)/2! + (tH)^3(ψ)/3! + ...
+ # = ψ + (tH) * (ψ + (tH)/2 * (ψ + (tH)/3 * (ψ + ...)))
+ eHψ = copy(ψ)
+ for ord in reverse(1:order)
+ eHψ = (t / ord) * H(eHψ) + ψ
+ end
+ return eHψ, (;)
end
diff --git a/src/solvers/operator_map.jl b/src/solvers/operator_map.jl
index 78a6f26..2f52f22 100644
--- a/src/solvers/operator_map.jl
+++ b/src/solvers/operator_map.jl
@@ -1,9 +1,9 @@
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 = contraction_sequence(contract_list; alg="optimal")
- Pψ = contract(contract_list; sequence)
- return noprime(Pψ)
+ 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 = contraction_sequence(contract_list; alg = "optimal")
+ Pψ = contract(contract_list; sequence)
+ return noprime(Pψ)
end
diff --git a/src/solvers/region_plans/dfs_plans.jl b/src/solvers/region_plans/dfs_plans.jl
index 9b44b98..7a5c32b 100644
--- a/src/solvers/region_plans/dfs_plans.jl
+++ b/src/solvers/region_plans/dfs_plans.jl
@@ -1,22 +1,22 @@
using Graphs: dst, src
using NamedGraphs.GraphsExtensions:
- default_root_vertex, post_order_dfs_edges, post_order_dfs_vertices
+ default_root_vertex, post_order_dfs_edges, post_order_dfs_vertices
function post_order_dfs_plan(
- graph, sweep_kwargs; nsites, root_vertex=default_root_vertex(graph)
-)
- 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
+ graph, sweep_kwargs; nsites, root_vertex = default_root_vertex(graph)
+ )
+ 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...]
+ 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
index 68548fd..aff6e3c 100644
--- a/src/solvers/region_plans/euler_plans.jl
+++ b/src/solvers/region_plans/euler_plans.jl
@@ -1,15 +1,15 @@
using Graphs: dst, src
using NamedGraphs.GraphsExtensions: default_root_vertex
-function euler_sweep(graph; nsites, root_vertex=default_root_vertex(graph), sweep_kwargs...)
- sweep_kwargs = (; nsites, root_vertex, sweep_kwargs...)
+function euler_sweep(graph; nsites, root_vertex = default_root_vertex(graph), sweep_kwargs...)
+ sweep_kwargs = (; nsites, root_vertex, 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
+ 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
index e5010d8..6ced448 100644
--- a/src/solvers/region_plans/euler_tour.jl
+++ b/src/solvers/region_plans/euler_tour.jl
@@ -2,47 +2,47 @@ using Graphs: dst, edges, src, vertices
using NamedGraphs: edgetype, vertextype
function compute_adjacencies(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))
- end
- return adj
+ 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))
+ end
+ return adj
end
function euler_tour_edges(G, start_vertex)
- adj = compute_adjacencies(G)
- etype = edgetype(G)
- vtype = 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
+ adj = compute_adjacencies(G)
+ etype = edgetype(G)
+ vtype = 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
- end
- return tour
+ 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)...]
+ 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
index 8ee9008..829fecd 100644
--- a/src/solvers/region_plans/tdvp_region_plans.jl
+++ b/src/solvers/region_plans/tdvp_region_plans.jl
@@ -1,60 +1,60 @@
using Accessors: @modify
function applyexp_sub_steps(order)
- if order == 1
- return [1.0]
- elseif order == 2
- return [1 / 2, 1 / 2]
- 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("Applyexp order of $order not supported")
- end
+ if order == 1
+ return [1.0]
+ elseif order == 2
+ return [1 / 2, 1 / 2]
+ 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("Applyexp order of $order not supported")
+ end
end
function first_order_sweep(graph, sweep_kwargs; nsites)
- basic_fwd_sweep = post_order_dfs_plan(graph, sweep_kwargs; nsites)
- region_plan = []
+ basic_fwd_sweep = post_order_dfs_plan(graph, sweep_kwargs; nsites)
+ region_plan = []
- for (j, (region, region_kwargs)) in enumerate(basic_fwd_sweep)
- push!(region_plan, region => region_kwargs)
+ for (j, (region, region_kwargs)) in enumerate(basic_fwd_sweep)
+ push!(region_plan, region => region_kwargs)
- if length(region) == 2 && j < length(basic_fwd_sweep)
- region_kwargs = @modify(-, region_kwargs.update!_kwargs.exponent_step)
- push!(region_plan, [last(region)] => region_kwargs)
+ if length(region) == 2 && j < length(basic_fwd_sweep)
+ region_kwargs = @modify(-, region_kwargs.update!_kwargs.exponent_step)
+ push!(region_plan, [last(region)] => region_kwargs)
+ end
end
- end
- return region_plan
+ return region_plan
end
function reverse_regions(region_plan)
- return map(reverse(region_plan)) do (region, kwargs)
- return reverse(region) => kwargs
- end
+ return map(reverse(region_plan)) do (region, kwargs)
+ return reverse(region) => kwargs
+ end
end
# Generate the kwargs for each region.
function applyexp_regions(
- graph, raw_exponent_step; order, nsites, update!_kwargs=(; nsites), remaining_kwargs...
-)
- sweep_plan = []
+ graph, raw_exponent_step; order, nsites, update!_kwargs = (; nsites), remaining_kwargs...
+ )
+ sweep_plan = []
- for (step, weight) in enumerate(applyexp_sub_steps(order))
- # Use this exponent step only if none provided
- new_update!_kwargs = (; exponent_step=weight * raw_exponent_step, update!_kwargs...)
+ for (step, weight) in enumerate(applyexp_sub_steps(order))
+ # Use this exponent step only if none provided
+ new_update!_kwargs = (; exponent_step = weight * raw_exponent_step, update!_kwargs...)
- sweep_kwargs = (; remaining_kwargs..., update!_kwargs=new_update!_kwargs)
+ sweep_kwargs = (; remaining_kwargs..., update!_kwargs = new_update!_kwargs)
- region_plan = first_order_sweep(graph, sweep_kwargs; nsites)
+ region_plan = first_order_sweep(graph, sweep_kwargs; nsites)
- if iseven(step)
- region_plan = reverse_regions(region_plan)
- end
+ if iseven(step)
+ region_plan = reverse_regions(region_plan)
+ end
- append!(sweep_plan, region_plan)
- end
+ append!(sweep_plan, region_plan)
+ end
- return sweep_plan
+ return sweep_plan
end
diff --git a/src/solvers/subspace/densitymatrix.jl b/src/solvers/subspace/densitymatrix.jl
index ae2ff50..d868cce 100644
--- a/src/solvers/subspace/densitymatrix.jl
+++ b/src/solvers/subspace/densitymatrix.jl
@@ -2,70 +2,70 @@ using NamedGraphs.GraphsExtensions: incident_edges
using Printf: @printf
function subspace_expand!(
- ::Backend"densitymatrix",
- region_iter,
- local_state;
- expansion_factor=1.5,
- maxexpand=typemax(Int),
- north_pass=1,
- eigen_kwargs=(;),
-)
- prob = problem(region_iter)
+ ::Backend"densitymatrix",
+ region_iter,
+ local_state;
+ expansion_factor = 1.5,
+ maxexpand = typemax(Int),
+ north_pass = 1,
+ eigen_kwargs = (;),
+ )
+ prob = problem(region_iter)
- region = current_region(region_iter)
- psi = copy(state(prob))
+ region = current_region(region_iter)
+ psi = copy(state(prob))
- prev_vertex_set = setdiff(pos(operator(prob)), region)
- (length(prev_vertex_set) != 1) && return region_iter, local_state
- prev_vertex = only(prev_vertex_set)
- A = psi[prev_vertex]
+ prev_vertex_set = setdiff(pos(operator(prob)), region)
+ (length(prev_vertex_set) != 1) && return region_iter, local_state
+ prev_vertex = only(prev_vertex_set)
+ A = psi[prev_vertex]
- next_vertices = filter(v -> (hascommoninds(psi[v], A)), region)
- isempty(next_vertices) && return region_iter, local_state
- next_vertex = only(next_vertices)
- C = psi[next_vertex]
+ next_vertices = filter(v -> (hascommoninds(psi[v], A)), region)
+ isempty(next_vertices) && return region_iter, local_state
+ next_vertex = only(next_vertices)
+ C = psi[next_vertex]
- a = commonind(A, C)
- isnothing(a) && return region_iter, local_state
- basis_size = prod(dim.(uniqueinds(A, C)))
+ a = commonind(A, C)
+ isnothing(a) && return region_iter, local_state
+ basis_size = prod(dim.(uniqueinds(A, C)))
- expanded_maxdim = compute_expansion(
- dim(a), basis_size; expansion_factor, maxexpand, eigen_kwargs.maxdim
- )
- expanded_maxdim <= 0 && return region_iter, local_state
+ expanded_maxdim = compute_expansion(
+ dim(a), basis_size; expansion_factor, maxexpand, eigen_kwargs.maxdim
+ )
+ expanded_maxdim <= 0 && return region_iter, local_state
- envs = environments(operator(prob))
- H = operator(operator(prob))
- sqrt_rho = A
- for e in incident_edges(operator(prob))
- (src(e) ∈ region || dst(e) ∈ region) && continue
- sqrt_rho *= envs[e]
- end
- sqrt_rho *= H[prev_vertex]
+ envs = environments(operator(prob))
+ H = operator(operator(prob))
+ sqrt_rho = A
+ for e in incident_edges(operator(prob))
+ (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 _ in 1:north_pass
- sqrt_rho = conj_proj_A(sqrt_rho)
- end
- rho = sqrt_rho * dag(noprime(sqrt_rho))
- D, U = eigen(rho; eigen_kwargs..., ishermitian=true)
+ conj_proj_A(T) = (T - prime(A) * (dag(prime(A)) * T))
+ for _ in 1:north_pass
+ sqrt_rho = conj_proj_A(sqrt_rho)
+ end
+ rho = sqrt_rho * dag(noprime(sqrt_rho))
+ D, U = eigen(rho; eigen_kwargs..., ishermitian = true)
- Uproj(T) = (T - prime(A, a) * (dag(prime(A, a)) * T))
- for _ 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 region_iter, local_state
- end
+ Uproj(T) = (T - prime(A, a) * (dag(prime(A, a)) * T))
+ for _ in 1:north_pass
+ U = Uproj(U)
+ end
+ if norm(dag(U) * A) > 1.0e-10
+ @printf("Warning: |U*A| = %.3E in subspace expansion\n", norm(dag(U) * A))
+ return region_iter, local_state
+ end
- Ax, ax = directsum(A => a, U => commonind(U, D))
- expander = dag(Ax) * A
- psi[prev_vertex] = Ax
- psi[next_vertex] = expander * C
- local_state = expander * local_state
+ Ax, ax = directsum(A => a, U => commonind(U, D))
+ expander = dag(Ax) * A
+ psi[prev_vertex] = Ax
+ psi[next_vertex] = expander * C
+ local_state = expander * local_state
- prob.state = psi
+ prob.state = psi
- return region_iter, local_state
+ return region_iter, local_state
end
diff --git a/src/solvers/subspace/ortho_subspace.jl b/src/solvers/subspace/ortho_subspace.jl
index 2646530..61b410f 100644
--- a/src/solvers/subspace/ortho_subspace.jl
+++ b/src/solvers/subspace/ortho_subspace.jl
@@ -4,7 +4,7 @@ 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]
+ return [q => expand_space(d, expansion_factor) for (q, d) in χs]
end
#
@@ -17,61 +17,61 @@ end
#
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 = prev_region(region_iterator)
- region = current_region(region_iterator)
- if isnothing(prev_region) || isa(region, AbstractEdge)
- return local_tensor
- end
+ ::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 = prev_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)
+ 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]
+ 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]
+ 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))
+ # 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")
+ 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")
+ 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) <= 1.0e-15 || expand_maxdim <= 0) && return local_tensor
+ Ux, S, V = svd(Y, basis_inds; cutoff = 1.0e-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
+ 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
+ return local_tensor
end
diff --git a/src/solvers/subspace/subspace.jl b/src/solvers/subspace/subspace.jl
index d538824..99ea3f9 100644
--- a/src/solvers/subspace/subspace.jl
+++ b/src/solvers/subspace/subspace.jl
@@ -1,41 +1,41 @@
using NDTensors: NDTensors
using NDTensors.BackendSelection: Backend, @Backend_str
-function subspace_expand!(region_iter, local_state; subspace_algorithm="nothing")
- backend = Backend(subspace_algorithm)
+function subspace_expand!(region_iter, local_state; subspace_algorithm = "nothing")
+ backend = Backend(subspace_algorithm)
- if backend isa Backend"nothing"
- return region_iter, local_state
- end
+ if backend isa Backend"nothing"
+ return region_iter, local_state
+ end
- _, local_state = subspace_expand!(
- backend, region_iter, local_state; region_kwargs(subspace_expand!, region_iter)...
- )
+ _, local_state = subspace_expand!(
+ backend, region_iter, local_state; region_kwargs(subspace_expand!, region_iter)...
+ )
- return region_iter, local_state
+ return region_iter, local_state
end
function subspace_expand!(backend, region_iterator, local_state; kwargs...)
- # We allow passing of any kwargs here is this method throws an error anyway
- return error(
- "Subspace expansion (subspace_expand!) not defined for requested combination of subspace_algorithm and problem types",
- )
+ # We allow passing of any kwargs here is this method throws an error anyway
+ return error(
+ "Subspace expansion (subspace_expand!) not defined for requested combination of subspace_algorithm and problem types",
+ )
end
# Have these defaults set per backend in `subspace_expand!`
function compute_expansion(current_dim, basis_size; expansion_factor, maxexpand, 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(maxexpand, 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 insert step)
- expand_maxdim = min(maxdim - current_dim, expand_maxdim)
- # Ensure expand_maxdim is non-negative
- expand_maxdim = max(0, expand_maxdim)
-
- return expand_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(maxexpand, 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 insert 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_solve.jl b/src/solvers/sweep_solve.jl
index f67ce4d..8084ee0 100644
--- a/src/solvers/sweep_solve.jl
+++ b/src/solvers/sweep_solve.jl
@@ -2,42 +2,42 @@
default_region_callback(sweep_iterator) = sweep_iterator
default_sweep_callback(sweep_iterator) = sweep_iterator
-# In this implementation the function `sweep_solve` is essentially just a wrapper around
+# In this implementation the function `sweep_solve` is essentially just a wrapper around
# the iterate interface that allows one to pass callbacks.
function sweep_solve!(
- sweep_iterator;
- sweep_callback=default_sweep_callback,
- region_callback=default_region_callback,
-)
- # Don't compute the region iteration automatically as we wish to insert a callback.
- for _ in IncrementOnly(sweep_iterator)
- for _ in region_iterator(sweep_iterator)
- region_callback(sweep_iterator)
+ sweep_iterator;
+ sweep_callback = default_sweep_callback,
+ region_callback = default_region_callback,
+ )
+ # Don't compute the region iteration automatically as we wish to insert a callback.
+ for _ in IncrementOnly(sweep_iterator)
+ for _ in region_iterator(sweep_iterator)
+ region_callback(sweep_iterator)
+ end
+ sweep_callback(sweep_iterator)
end
- sweep_callback(sweep_iterator)
- end
- return sweep_iterator
+ return sweep_iterator
end
# I suspect that `sweep_callback` is the more commonly used callback, so allow this to
# be set using the `do` syntax.
function sweep_solve!(
- sweep_callback, sweep_iterator; region_callback=default_region_callback
-)
- return sweep_solve!(sweep_iterator; sweep_callback, region_callback)
+ sweep_callback, sweep_iterator; region_callback = default_region_callback
+ )
+ return sweep_solve!(sweep_iterator; sweep_callback, region_callback)
end
function sweep_solve!(
- each_region_iterator::EachRegion; region_callback=default_region_callback
-)
- return sweep_solve!(region_callback, each_region_iterator)
+ each_region_iterator::EachRegion; region_callback = default_region_callback
+ )
+ return sweep_solve!(region_callback, each_region_iterator)
end
function sweep_solve!(region_callback, each_region_iterator::EachRegion)
- for _ in each_region_iterator
- # I don't think it is obvious what object this particular callback should take,
- # but for now be consistant and pass the parent sweep iterator.
- sweep_iterator = each_region_iterator.parent
- region_callback(sweep_iterator)
- end
- return each_region_iterator
+ for _ in each_region_iterator
+ # I don't think it is obvious what object this particular callback should take,
+ # but for now be consistant and pass the parent sweep iterator.
+ sweep_iterator = each_region_iterator.parent
+ region_callback(sweep_iterator)
+ end
+ return each_region_iterator
end
diff --git a/src/solvers/truncation_parameters.jl b/src/solvers/truncation_parameters.jl
index bcd7d94..6fb979b 100644
--- a/src/solvers/truncation_parameters.jl
+++ b/src/solvers/truncation_parameters.jl
@@ -5,10 +5,10 @@ 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)
+ 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
diff --git a/test/runtests.jl b/test/runtests.jl
index fb2673d..fbde7af 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,7 +24,7 @@ end
@time begin
# tests in groups based on folder structure
- for testgroup in filter(f -> isdir(joinpath(@__DIR__, f)), readdir(@__DIR__))
+ for testgroup in filter(f -> isdir(joinpath(@__DIR__, f)), readdir(@__DIR__))
if GROUP == "ALL" || GROUP == uppercase(testgroup)
groupdir = joinpath(@__DIR__, testgroup)
for file in filter(istestfile, readdir(groupdir))
diff --git a/test/solvers/test_applyexp.jl b/test/solvers/test_applyexp.jl
index c3464c8..79b3700 100644
--- a/test/solvers/test_applyexp.jl
+++ b/test/solvers/test_applyexp.jl
@@ -8,115 +8,115 @@ using ITensorMPS: OpSum
using TensorOperations: TensorOperations #For contraction order finding
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
+ 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 "Time Evolution" begin
- @testset "Test Tree Time Evolution" begin
- outputlevel = 0
-
- 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
- factorize_kwargs = (; cutoff, maxdim)
- E, gs_psi = dmrg(H, psi0; factorize_kwargs, nsites, nsweeps, outputlevel)
- (outputlevel >= 1) && println("2-site DMRG energy = ", E)
-
- nsites = 1
- tmax = 0.10
- time_range = 0.0:0.02:tmax
- psi1_t = time_evolve(H, time_range, gs_psi; factorize_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; factorize_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 atan(imag(z) / real(z)) ≈ E * tmax atol = 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(...*[Comment body truncated]* |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
4 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR is a rewrite of codes for "sweep solver" algorithms such as DMRG, TDVP, etc.
It introduces a simplified designs for the whole system, especially regarding the creation of "region plans", how the double loop comprising what is now called
sweep_solveis coded, and the handling of keyword arguments.The PR also adds a subspace expansion system with one currently implementation (based on a "projected density matrix perturbation" idea, basically McCulloch+White hybrid method). Having a subspace expansion is crucial to converge DMRG on certain tree lattices, not to mention other cases like 2D DMRG with QN conservation.
Internally, the core design revolves around two iterators: a sweep iterator and a region iterator. Each of these is in principle wrappable by "iteration adapters" (see for example the tuple adapter in adapters.jl). At each iteration of the RegionIterator type, it calls a
region_iterator_actionfunction which can be overloaded, but by default calls three "subactions":extracter,updater, andinserter. Currentlyextracteralso calls down into asubspace_expandfunction which can be customized through multiple backends. These "action" functions all dispatch on a "problem" type which can hold arbitrary data, making these codes more flexible and future-proof toward cases like optimizing two tensor network states at once, working with sets of operators, etc.Other improvements that may fit better into future PRs:
truncateandapplyalready in the NetworkSolvers repo (there were some issues bringing them over)dmrgwhich automatically propagate truncation arguments into more "expert" keyword argument packs