From bc155dee0f91054cca25d44ed4fbed20399bd528 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 16 Apr 2024 09:28:37 -0400 Subject: [PATCH] Backward_shift_op definition. Double shift tests --- examples/2d_laplace_solver.jl | 5 +- src/elementary_operators.jl | 18 +++- test/test_shift_operators.jl | 189 ++++++++++++++-------------------- 3 files changed, 94 insertions(+), 118 deletions(-) diff --git a/examples/2d_laplace_solver.jl b/examples/2d_laplace_solver.jl index 7f4838c..8fdb01f 100644 --- a/examples/2d_laplace_solver.jl +++ b/examples/2d_laplace_solver.jl @@ -12,15 +12,14 @@ using UnicodePlots #Solve the 2D Laplace equation on a random tree seed!(1234) -L = 14 +L = 12 g = NamedGraph(SimpleGraph(uniform_tree(L))) bit_map = BitMap(g; map_dimension=2) s = siteinds(g, bit_map) ψ_fxy = 0.1 * rand_itn(s, bit_map; link_space=2) -∇ = laplacian_operator(s, bit_map; scale=false) -∇ = truncate(∇; cutoff=1e-8) +∇ = laplacian_operator(s, bit_map; scale=false, cutoff=1e-8) println("2D Laplacian constructed for this tree, bond dimension is $(maxlinkdim(∇))") init_energy = diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index b9461ae..7522480 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -45,7 +45,7 @@ function ITensors.op(::OpName"Dup", ::SiteType"Digit", s::Index) return ITensor(o, s, s') end -function plus_shift_opsum( +function forward_shift_opsum( s::IndsNetwork, bit_map; dimension=default_dimension(), @@ -90,7 +90,7 @@ function plus_shift_opsum( return ttn_op end -function minus_shift_opsum( +function backward_shift_opsum( s::IndsNetwork, bit_map; dimension=default_dimension(), @@ -142,6 +142,16 @@ function no_shift_opsum(s::IndsNetwork) return ttn_op end +function backward_shift_op(s::IndsNetwork, bit_map::BitMap; truncate_kwargs=(;), kwargs...) + ttn_opsum = backward_shift_opsum(s, bit_map; kwargs...) + return ttn(ttn_opsum, s; algorithm="svd", truncate_kwargs...) +end + +function forward_shift_op(s::IndsNetwork, bit_map::BitMap; truncate_kwargs=(;), kwargs...) + ttn_opsum = forward_shift_opsum(s, bit_map; kwargs...) + return ttn(ttn_opsum, s; algorithm="svd", truncate_kwargs...) +end + function stencil( s::IndsNetwork, bit_map, @@ -160,7 +170,7 @@ function stencil( n = i == 1 ? 1 : 0 if !iszero(shifts[i]) stencil_opsum += - shifts[i] * plus_shift_opsum(s, bit_map; dimension, boundary=right_boundary, n) + shifts[i] * forward_shift_opsum(s, bit_map; dimension, boundary=right_boundary, n) end end @@ -168,7 +178,7 @@ function stencil( n = i == 5 ? 1 : 0 if !iszero(shifts[i]) stencil_opsum += - shifts[i] * minus_shift_opsum(s, bit_map; dimension, boundary=left_boundary, n) + shifts[i] * backward_shift_opsum(s, bit_map; dimension, boundary=left_boundary, n) end end diff --git a/test/test_shift_operators.jl b/test/test_shift_operators.jl index 07aa9fa..c48b5c9 100644 --- a/test/test_shift_operators.jl +++ b/test/test_shift_operators.jl @@ -8,10 +8,10 @@ using NamedGraphs: named_grid, named_comb_tree, NamedGraph, nv, vertices using ITensorNumericalAnalysis: itensornetwork using Dictionaries: Dictionary -using ITensorNumericalAnalysis: stencil +using ITensorNumericalAnalysis: backward_shift_op, forward_shift_op -@testset "test shift operators in 1D on MPS" begin - g = named_grid((6, 1)) +@testset "test shift operators in 1D on Tree" begin + g = named_comb_tree((2, 3)) L = nv(g) delta = 2.0^(-1.0 * L) bit_map = BitMap(g) @@ -19,30 +19,31 @@ using ITensorNumericalAnalysis: stencil xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] ψ_fx = poly_itn(s, bit_map, [1.0, 0.5, 0.25]) - plus_shift_dirichlet = stencil( - s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Dirichlet" + forward_shift_dirichlet = forward_shift_op( + s, bit_map; boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) ) - minus_shift_dirichlet = stencil( - s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Dirichlet" + backward_shift_dirichlet = backward_shift_op( + s, bit_map; boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) ) - plus_shift_pbc = stencil( - s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Periodic" + forward_shift_pbc = forward_shift_op( + s, bit_map; boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) ) - minus_shift_pbc = stencil( - s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Periodic" + backward_shift_pbc = backward_shift_op( + s, bit_map; boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) ) - plus_shift_neumann = stencil( - s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Neumann" + forward_shift_neumann = forward_shift_op( + s, bit_map; boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) ) - minus_shift_neumann = stencil( - s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Neumann" + backward_shift_neumann = backward_shift_op( + s, bit_map; boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) ) - ψ_fx_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_neumann = operate(minus_shift_neumann, ψ_fx; cutoff=1e-12) + + ψ_fx_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_pbc = operate(forward_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_pbc = operate(backward_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_neumann = operate(forward_shift_neumann, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_neumann = operate(backward_shift_neumann, ψ_fx; cutoff=1e-12) for x in xs if x + delta < 1 @@ -70,7 +71,7 @@ using ITensorNumericalAnalysis: stencil end end -@testset "test shift operators in 1D on Tree" begin +@testset "test double shift operators in 1D on Tree" begin g = named_comb_tree((2, 3)) L = nv(g) delta = 2.0^(-1.0 * L) @@ -78,55 +79,57 @@ end s = siteinds(g, bit_map) xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta] ψ_fx = poly_itn(s, bit_map, [1.0, 0.5, 0.25]) + n = 1 - plus_shift_dirichlet = stencil( - s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Dirichlet" + forward_shift_dirichlet = forward_shift_op( + s, bit_map; n, boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) ) - minus_shift_dirichlet = stencil( - s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Dirichlet" + backward_shift_dirichlet = backward_shift_op( + s, bit_map; n, boundary="Dirichlet", truncate_kwargs=(; cutoff=1e-10) ) - plus_shift_pbc = stencil( - s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Periodic" + forward_shift_pbc = forward_shift_op( + s, bit_map; n, boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) ) - minus_shift_pbc = stencil( - s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Periodic" + backward_shift_pbc = backward_shift_op( + s, bit_map; n, boundary="Periodic", truncate_kwargs=(; cutoff=1e-10) ) - plus_shift_neumann = stencil( - s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Neumann" + forward_shift_neumann = forward_shift_op( + s, bit_map; n, boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) ) - minus_shift_neumann = stencil( - s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Neumann" + backward_shift_neumann = backward_shift_op( + s, bit_map; n, boundary="Neumann", truncate_kwargs=(; cutoff=1e-10) ) - ψ_fx_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; cutoff=1e-12) - ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; cutoff=1e-12) - ψ_fx_mshift_neumann = operate(minus_shift_neumann, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_pbc = operate(forward_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_pbc = operate(backward_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_neumann = operate(forward_shift_neumann, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_neumann = operate(backward_shift_neumann, ψ_fx; cutoff=1e-12) for x in xs - if x + delta < 1 - fx_xplus = calculate_fx(ψ_fx, x + delta) + if x + 2.0 * delta < 1 + fx_xplus = calculate_fx(ψ_fx, x + 2.0 * delta) @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_dirichlet, x) atol = 1e-8 @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_pbc, x) atol = 1e-8 @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_neumann, x) atol = 1e-8 - elseif x == 1.0 - delta + elseif x == 1.0 - 2.0 * delta || x == 1.0 - delta @test calculate_fx(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1e-8 - @test calculate_fx(ψ_fx_pshift_pbc, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8 - @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = + @test calculate_fx(ψ_fx_pshift_pbc, x) ≈ calculate_fx(ψ_fx, x + 2.0 * delta - 1.0) atol = 1e-8 + @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, x) atol = 1e-8 end - if x - delta >= 0.0 - fx_xminus = calculate_fx(ψ_fx, x - delta) + if x - 2.0 * delta >= 0.0 + fx_xminus = calculate_fx(ψ_fx, x - 2.0 * delta) @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_dirichlet, x) atol = 1e-8 @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_pbc, x) atol = 1e-8 @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_neumann, x) atol = 1e-8 - elseif x == 0.0 + else @test calculate_fx(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8 - @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8 - @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8 + @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, x - 2.0 * delta + 1.0) atol = + 1e-8 + @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, x) atol = 1e-8 end end end @@ -143,67 +146,31 @@ end ψ_fy = cos_itn(s, bit_map; dimension=2) ψ_fxy = ψ_fx + ψ_fx - plus_shift_dirichlet = stencil( - s, - bit_map, - [0.0, 1.0, 0.0, 0.0, 0.0], - 0; - scale=false, - dimension=2, - right_boundary="Dirichlet", - ) - minus_shift_dirichlet = stencil( - s, - bit_map, - [0.0, 0.0, 0.0, 1.0, 0.0], - 0; - scale=false, - dimension=2, - left_boundary="Dirichlet", - ) - plus_shift_pbc = stencil( - s, - bit_map, - [0.0, 1.0, 0.0, 0.0, 0.0], - 0; - scale=false, - dimension=2, - right_boundary="Periodic", - ) - minus_shift_pbc = stencil( - s, - bit_map, - [0.0, 0.0, 0.0, 1.0, 0.0], - 0; - scale=false, - dimension=2, - left_boundary="Periodic", - ) - plus_shift_neumann = stencil( - s, - bit_map, - [0.0, 1.0, 0.0, 0.0, 0.0], - 0; - scale=false, - dimension=2, - right_boundary="Neumann", - ) - minus_shift_neumann = stencil( - s, - bit_map, - [0.0, 0.0, 0.0, 1.0, 0.0], - 0; - scale=false, - dimension=2, - left_boundary="Neumann", - ) - - ψ_fxy_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fxy; cutoff=1e-12) - ψ_fxy_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fxy; cutoff=1e-12) - ψ_fxy_pshift_pbc = operate(plus_shift_pbc, ψ_fxy; cutoff=1e-12) - ψ_fxy_mshift_pbc = operate(minus_shift_pbc, ψ_fxy; cutoff=1e-12) - ψ_fxy_pshift_neumann = operate(plus_shift_neumann, ψ_fxy; cutoff=1e-12) - ψ_fxy_mshift_neumann = operate(minus_shift_neumann, ψ_fxy; cutoff=1e-12) + forward_shift_dirichlet = forward_shift_op( + s, bit_map; boundary="Dirichlet", dimension=2, truncate_kwargs=(; cutoff=1e-10) + ) + backward_shift_dirichlet = backward_shift_op( + s, bit_map; boundary="Dirichlet", dimension=2, truncate_kwargs=(; cutoff=1e-10) + ) + forward_shift_pbc = forward_shift_op( + s, bit_map; boundary="Periodic", dimension=2, truncate_kwargs=(; cutoff=1e-10) + ) + backward_shift_pbc = backward_shift_op( + s, bit_map; boundary="Periodic", dimension=2, truncate_kwargs=(; cutoff=1e-10) + ) + forward_shift_neumann = forward_shift_op( + s, bit_map; boundary="Neumann", dimension=2, truncate_kwargs=(; cutoff=1e-10) + ) + backward_shift_neumann = backward_shift_op( + s, bit_map; boundary="Neumann", dimension=2, truncate_kwargs=(; cutoff=1e-10) + ) + + ψ_fxy_pshift_dirichlet = operate(forward_shift_dirichlet, ψ_fxy; cutoff=1e-12) + ψ_fxy_mshift_dirichlet = operate(backward_shift_dirichlet, ψ_fxy; cutoff=1e-12) + ψ_fxy_pshift_pbc = operate(forward_shift_pbc, ψ_fxy; cutoff=1e-12) + ψ_fxy_mshift_pbc = operate(backward_shift_pbc, ψ_fxy; cutoff=1e-12) + ψ_fxy_pshift_neumann = operate(forward_shift_neumann, ψ_fxy; cutoff=1e-12) + ψ_fxy_mshift_neumann = operate(backward_shift_neumann, ψ_fxy; cutoff=1e-12) for y in ys if y + delta < 1