Backward_shift_op definition. Double shift tests

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 =

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,15 +170,15 @@ 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
 
   for i in [4, 5]
     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
 

test/test_shift_operators.jl

@@ -8,41 +8,42 @@ 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)
   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])
 
-  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,63 +71,65 @@ 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)
   bit_map = BitMap(g)
   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