Small fixes and new features

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2024 Joseph Tindall <[email protected]> and contributors
+Copyright (c) 2024 Joseph Tindall <[email protected]>, Ryan Levy <[email protected]> and contributors
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

Project.toml

@@ -16,6 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 ITensorNetworks = "0.8.1"
 Glob = "1.3.1"
+Graphs = "1.8"
 Dictionaries = "0.4.2"
 julia = "1.8"
 ITensors = "0.3.58"

examples/2d_laplace_solver.jl

@@ -17,7 +17,7 @@ g = NamedGraph(SimpleGraph(uniform_tree(L)))
 
 s = continuous_siteinds(g; map_dimension=2)
 
-ψ_fxy = 0.1 * rand_itn(s; link_space=2)
+ψ_fxy = 0.1 * rand_itn(s; link_space=3)
 ∇ = laplacian_operator(s; scale=false, cutoff=1e-8)
 println("2D Laplacian constructed for this tree, bond dimension is $(maxlinkdim(∇))")
 

src/ITensorNumericalAnalysis.jl

@@ -42,7 +42,9 @@ export const_itensornetwork,
   second_derivative_operator,
   third_derivative_operator,
   fourth_derivative_operator,
-  identity_operator
+  identity_operator,
+  delta_x,
+  delta_xyz
 export const_itn,
   poly_itn, cosh_itn, sinh_itn, tanh_itn, exp_itn, sin_itn, cos_itn, rand_itn
 export calculate_fx, calculate_fxyz

src/elementary_functions.jl

@@ -201,6 +201,23 @@ function random_itensornetwork(s::IndsNetworkMap; kwargs...)
   return ITensorNetworkFunction(random_tensornetwork(indsnetwork(s); kwargs...), s)
 end
 
+"Create a product state of a given bit configuration"
+function delta_xyz(s::IndsNetworkMap, xs::Vector, dimensions::Vector{Int}; kwargs...)
+  ind_to_ind_value_map = calculate_ind_values(s, xs, dimensions)
+  tn = ITensorNetwork(v -> string(ind_to_ind_value_map[only(s[v])]), indsnetwork(s))
+  return ITensorNetworkFunction(tn, s)
+end
+
+function delta_xyz(s::IndsNetworkMap, xs::Vector; kwargs...)
+  return delta_xyz(s, xs, [i for i in 1:length(xs)]; kwargs...)
+end
+
+"Create a product state of a given bit configuration of a 1D function"
+function delta_x(s::IndsNetworkMap, x::Number, kwargs...)
+  @assert dimension(s) == 1
+  return delta_xyz(s, [x], [1]; kwargs...)
+end
+
 const const_itn = const_itensornetwork
 const poly_itn = polynomial_itensornetwork
 const cosh_itn = cosh_itensornetwork

src/elementary_operators.jl

@@ -2,12 +2,11 @@ using Graphs: is_tree
 using NamedGraphs: undirected_graph
 using ITensors:
   OpSum,
-  @OpName_str,
-  @SiteType_str,
   SiteType,
   siteinds,
   noprime,
   op,
+  Op,
   truncate,
   replaceinds,
   delta,
@@ -16,57 +15,23 @@ using ITensors:
   sim,
   noprime!,
   contract
-using ITensorNetworks: IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn
-
+using ITensorNetworks:
+  IndsNetwork, ITensorNetwork, TreeTensorNetwork, combine_linkinds, ttn, union_all_inds
 default_boundary() = "Dirichlet"
 
-function ITensors.op(::OpName"D+", ::SiteType"Digit", s::Index)
-  d = dim(s)
-  o = zeros(d, d)
-  o[2, 1] = 1
-  return ITensor(o, s, s')
-end
-function ITensors.op(::OpName"D-", ::SiteType"Digit", s::Index)
-  d = dim(s)
-  o = zeros(d, d)
-  o[1, 2] = 1
-  return ITensor(o, s, s')
-end
-function ITensors.op(::OpName"Ddn", ::SiteType"Digit", s::Index)
-  d = dim(s)
-  o = zeros(d, d)
-  o[1, 1] = 1
-  return ITensor(o, s, s')
-end
-function ITensors.op(::OpName"Dup", ::SiteType"Digit", s::Index)
-  d = dim(s)
-  o = zeros(d, d)
-  o[2, 2] = 1
-  return ITensor(o, s, s')
-end
+## TODO: turn this into a proper system ala sites which can be externally overloaded
 
-function forward_shift_opsum(
-  s::IndsNetworkMap; dimension=default_dimension(), boundary=default_boundary(), n::Int=0
+function boundary_term(
+  s::IndsNetworkMap, boundary::String, dimension, isfwd::Bool, n::Int=0
 )
-  @assert is_tree(s)
-  @assert base(s) == 2
   ttn_op = OpSum()
   dim_vertices = dimension_vertices(s, dimension)
   L = length(dim_vertices)
 
-  string_site = [("D+", vertex(s, dimension, L - n))]
-  add!(ttn_op, 1.0, "D+", vertex(s, dimension, L - n))
-  for i in (L - n):-1:2
-    pop!(string_site)
-    push!(string_site, ("D-", vertex(s, dimension, i)))
-    push!(string_site, ("D+", vertex(s, dimension, i - 1)))
-    add!(ttn_op, 1.0, (string_site...)...)
-  end
-
   if boundary == "Neumann"
     string_site = [
       if j <= (L - n)
-        ("Dup", vertex(s, dimension, j))
+        (isfwd ? "Dup" : "Ddn", vertex(s, dimension, j))
       else
         ("I", vertex(s, dimension, j))
       end for j in 1:L
@@ -75,13 +40,35 @@ function forward_shift_opsum(
   elseif boundary == "Periodic"
     string_site = [
       if j <= (L - n)
-        ("D-", vertex(s, dimension, j))
+        (isfwd ? "D-" : "D+", vertex(s, dimension, j))
       else
         ("I", vertex(s, dimension, j))
       end for j in 1:L
     ]
     add!(ttn_op, 1.0, (string_site...)...)
   end
+  return ttn_op
+end
+
+function forward_shift_opsum(
+  s::IndsNetworkMap; dimension=default_dimension(), boundary=default_boundary(), n::Int=0
+)
+  @assert is_tree(s)
+  @assert base(s) == 2
+  ttn_op = OpSum()
+  dim_vertices = dimension_vertices(s, dimension)
+  L = length(dim_vertices)
+
+  string_site = [("D+", vertex(s, dimension, L - n))]
+  add!(ttn_op, 1.0, "D+", vertex(s, dimension, L - n))
+  for i in (L - n):-1:2
+    pop!(string_site)
+    push!(string_site, ("D-", vertex(s, dimension, i)))
+    push!(string_site, ("D+", vertex(s, dimension, i - 1)))
+    add!(ttn_op, 1.0, (string_site...)...)
+  end
+
+  ttn_op += boundary_term(s, boundary, dimension, true, n)
 
   return ttn_op
 end
@@ -104,25 +91,7 @@ function backward_shift_opsum(
     add!(ttn_op, 1.0, (string_site...)...)
   end
 
-  if boundary == "Neumann"
-    string_site = [
-      if j <= (L - n)
-        ("Ddn", vertex(s, dimension, j))
-      else
-        ("I", vertex(s, dimension, j))
-      end for j in 1:L
-    ]
-    add!(ttn_op, 1.0, (string_site...)...)
-  elseif boundary == "Periodic"
-    string_site = [
-      if j <= (L - n)
-        ("D+", vertex(s, dimension, j))
-      else
-        ("I", vertex(s, dimension, j))
-      end for j in 1:L
-    ]
-    add!(ttn_op, 1.0, (string_site...)...)
-  end
+  ttn_op += boundary_term(s, boundary, dimension, false, n)
 
   return ttn_op
 end
@@ -155,6 +124,7 @@ function stencil(
   truncate_op=true,
   kwargs...,
 )
+  # shifts = [ x+2Δh, x+Δh, x, x-Δh, x-2Δh]
   @assert length(shifts) == 5
   b = base(s)
   stencil_opsum = shifts[3] * no_shift_opsum(s)
@@ -212,9 +182,13 @@ function laplacian_operator(
   end
   return ∇
 end
+function laplacian_operator(s::IndsNetworkMap, boundary::String; kwargs...)
+  return laplacian_operator(s; left_boundary=boundary, right_boundary=boundary, kwargs...)
+end
 
 function identity_operator(s::IndsNetworkMap; kwargs...)
-  return stencil(s, [0.0, 1.0, 0.0], 0; kwargs...)
+  operator_inds = ITensorNetworks.union_all_inds(indsnetwork(s), prime(indsnetwork(s)))
+  return ITensorNetwork(Op("I"), operator_inds)
 end
 
 function operator(fx::ITensorNetworkFunction)