Merge pull request #37 from JoeyT1994/fix_multiply

Fix bugs in multiply and identity_operator

src/elementary_operators.jl

@@ -15,9 +15,16 @@ using ITensors:
   prime,
   sim,
   noprime!,
-  contract
+  contract,
+  replaceinds
 using ITensorNetworks:
-  IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn, union_all_inds
+  IndsNetwork,
+  ITensorNetwork,
+  TreeTensorNetwork,
+  combine_linkinds,
+  ttn,
+  union_all_inds,
+  map_inds
 default_boundary() = "Dirichlet"
 
 ## TODO: turn this into a proper system ala sites which can be externally overloaded
@@ -183,7 +190,7 @@ end
 
 function identity_operator(s::IndsNetworkMap; kwargs...)
   operator_inds = ITensorNetworks.union_all_inds(indsnetwork(s), prime(indsnetwork(s)))
-  return ITensorNetwork(Op("I"), operator_inds)
+  return ITensorNetwork(Op("I"), operator_inds; kwargs...)
 end
 
 " Create an operator bitstring corresponding to the number x"
@@ -250,14 +257,17 @@ function operator_proj(fx::ITensorNetworkFunction)
   return operator
 end
 
+#gx and fx need to live in different virtual spaces (can't share linkinds)
 function multiply(gx::ITensorNetworkFunction, fx::ITensorNetworkFunction)
-  @assert vertices(gx) == vertices(fx)
   fx, fxgx = copy(fx), copy(gx)
+  fx = map_inds(prime, fx; sites=[])
+  @assert vertices(gx) == vertices(fx)
   s = siteinds(fxgx)
   for v in vertices(fxgx)
     ssim = sim(s[v])
     temp_tensor = replaceinds(fx[v], s[v], ssim)
-    fxgx[v] = noprime!(fxgx[v] * delta(s[v], s[v]', ssim) * temp_tensor)
+    fxgx[v] = fxgx[v] * delta(s[v], s[v]', ssim) * temp_tensor
+    fxgx[v] = replaceinds(fxgx[v], s[v]', s[v])
   end
 
   return combine_linkinds(fxgx)

test/test_operators.jl

@@ -104,10 +104,13 @@ using Dictionaries: Dictionary
     ψ_fx = cos_itn(s; k=0.25 * Number(pi))
 
     ψ_fxgx = ψ_gx * ψ_fx
+    ψ_sq = ψ_fx * ψ_fx
     xs = [0.025, 0.1, 0.25, 0.625, 0.875]
     for x in xs
       ψ_fxgx_x = real(evaluate(ψ_fxgx, x))
       @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-3
+      ψ_sq_x = real(evaluate(ψ_sq, x))
+      @test ψ_sq_x ≈ cos(0.25 * pi * x) * cos(0.25 * pi * x) atol = 1e-3
     end
   end
 
@@ -123,6 +126,7 @@ using Dictionaries: Dictionary
     @assert dimension(ψ_fx) == dimension(ψ_gy) == 2
 
     ψ_fxgy = ψ_fx * ψ_gy
+    ψ_sq = ψ_gy * ψ_gy
 
     xs = [0.125, 0.25, 0.625, 0.875]
     ys = [0.125, 0.25, 0.625, 0.875]
@@ -134,10 +138,32 @@ using Dictionaries: Dictionary
         @test ψ_gy_y ≈ sin(0.5 * pi * y)
         ψ_fxgy_xy = real(evaluate(ψ_fxgy, [x, y]))
         @test ψ_fxgy_xy ≈ cos(0.25 * pi * x) * sin(0.5 * pi * y) atol = 1e-3
+        ψ_sq_xy = real(evaluate(ψ_sq, [x, y]))
+        @test ψ_sq_xy ≈ (sin(0.5 * pi * y))^2 atol = 1e-3
       end
     end
   end
 
+  @testset "test operator_proj in 1D" begin
+    g = named_comb_tree((4, 3))
+    L = nv(g)
+    s = continuous_siteinds(g)
+
+    ψ_gx = sin_itn(s; k=0.5 * Number(pi))
+    ψ_fx = cos_itn(s; k=0.25 * Number(pi))
+    O = operator_proj(ψ_fx)
+
+    ψ_fxgx = operate(O, ψ_gx; cutoff=1e-14)
+    ψ_sq = operate(O, ψ_fx; cutoff=1e-14)
+    xs = [0.025, 0.1, 0.25, 0.625, 0.875]
+    for x in xs
+      ψ_fxgx_x = real(evaluate(ψ_fxgx, x))
+      @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-3
+      ψ_sq_x = real(evaluate(ψ_sq, x))
+      @test ψ_sq_x ≈ cos(0.25 * pi * x) * cos(0.25 * pi * x) atol = 1e-3
+    end
+  end
+
   @testset "test shift operators in 1D on Tree" begin
     g = named_comb_tree((2, 3))
     L = nv(g)