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

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)

src/indexmap.jl

@@ -8,8 +8,8 @@ struct IndexMap{VB,VD}
   index_dimension::VD
 end
 
-index_digit(im::IndexMap) = im.index_digit
-index_dimension(im::IndexMap) = im.index_dimension
+index_digit(imap::IndexMap) = imap.index_digit
+index_dimension(imap::IndexMap) = imap.index_dimension
 
 function default_digit_map(indices::Vector{Index}; map_dimension::Int=1)
   return Dictionary(
@@ -48,83 +48,83 @@ function IndexMap(dimension_indices::Vector{Vector{V}}) where {V<:Index}
   return IndexMap(index_digit, index_dimension)
 end
 
-function Base.copy(im::IndexMap)
-  return IndexMap(copy(index_digit(im)), copy(index_dimension(im)))
+function Base.copy(imap::IndexMap)
+  return IndexMap(copy(index_digit(imap)), copy(index_dimension(imap)))
 end
 
-dimension(im::IndexMap) = maximum(collect(values(index_dimension(im))))
-dimension(im::IndexMap, ind::Index) = index_dimension(im)[ind]
-dimensions(im::IndexMap, inds::Vector{Index}) = dimension.(inds)
-digit(im::IndexMap, ind::Index) = index_digit(im)[ind]
-digits(im::IndexMap, inds::Vector{Index}) = digit.(inds)
-function index_value_to_scalar(im::IndexMap, ind::Index, value::Int)
-  return (value) / (dim(ind)^digit(im, ind))
+dimension(imap::IndexMap) = maximum(collect(values(index_dimension(imap))))
+dimension(imap::IndexMap, ind::Index) = index_dimension(imap)[ind]
+dimensions(imap::IndexMap, inds::Vector{Index}) = dimension.(inds)
+digit(imap::IndexMap, ind::Index) = index_digit(imap)[ind]
+digits(imap::IndexMap, inds::Vector{Index}) = digit.(inds)
+function index_value_to_scalar(imap::IndexMap, ind::Index, value::Int)
+  return (value) / (dim(ind)^digit(imap, ind))
 end
-function index_values_to_scalars(im::IndexMap, ind::Index)
-  return [index_value_to_scalar(im, ind, i) for i in 0:(dim(ind) - 1)]
+function index_values_to_scalars(imap::IndexMap, ind::Index)
+  return [index_value_to_scalar(imap, ind, i) for i in 0:(dim(ind) - 1)]
 end
 
-function ITensors.inds(im::IndexMap)
-  @assert keys(index_dimension(im)) == keys(index_digit(im))
-  return collect(keys(index_dimension(im)))
+function ITensors.inds(imap::IndexMap)
+  @assert keys(index_dimension(imap)) == keys(index_digit(imap))
+  return collect(keys(index_dimension(imap)))
 end
-function dimension_inds(im::IndexMap, dimension::Int)
+function dimension_inds(imap::IndexMap, dimension::Int)
   return collect(
-    filter(i -> index_dimension(im)[i] == dimension, keys(index_dimension(im)))
+    filter(i -> index_dimension(imap)[i] == dimension, keys(index_dimension(imap)))
   )
 end
-function ind(im::IndexMap, dimension::Int, digit::Int)
+function ind(imap::IndexMap, dimension::Int, digit::Int)
   return only(
     filter(
-      i -> index_dimension(im)[i] == dimension && index_digit(im)[i] == digit,
-      keys(index_dimension(im)),
+      i -> index_dimension(imap)[i] == dimension && index_digit(imap)[i] == digit,
+      keys(index_dimension(imap)),
     ),
   )
 end
 
-function dimension_indices(im::IndexMap, dim::Int)
-  return filter(ind -> dimension(im, ind) == dim, inds(im))
+function dimension_indices(imap::IndexMap, dimap::Int)
+  return filter(ind -> dimension(imap, ind) == dim, inds(imap))
 end
 
-function calculate_xyz(im::IndexMap, ind_to_ind_value_map, dimensions::Vector{Int})
+function calculate_xyz(imap::IndexMap, ind_to_ind_value_map, dimensions::Vector{Int})
   out = Number[]
   for dimension in dimensions
-    indices = dimension_inds(im, dimension)
+    indices = dimension_inds(imap, dimension)
     push!(
       out,
-      sum([index_value_to_scalar(im, ind, ind_to_ind_value_map[ind]) for ind in indices]),
+      sum([index_value_to_scalar(imap, ind, ind_to_ind_value_map[ind]) for ind in indices]),
     )
   end
   return out
 end
 
-function calculate_xyz(im::IndexMap, ind_to_ind_value_map)
-  return calculate_xyz(im, ind_to_ind_value_map, [i for i in 1:dimension(im)])
+function calculate_xyz(imap::IndexMap, ind_to_ind_value_map)
+  return calculate_xyz(imap, ind_to_ind_value_map, [i for i in 1:dimension(imap)])
 end
-function calculate_x(im::IndexMap, ind_to_ind_value_map, dimension::Int)
-  return only(calculate_xyz(im, ind_to_ind_value_map, [dimension]))
+function calculate_x(imap::IndexMap, ind_to_ind_value_map, dimension::Int)
+  return only(calculate_xyz(imap, ind_to_ind_value_map, [dimension]))
 end
-function calculate_x(im::IndexMap, ind_to_ind_value_map)
-  return calculate_x(im, ind_to_ind_value_map, 1)
+function calculate_x(imap::IndexMap, ind_to_ind_value_map)
+  return calculate_x(imap, ind_to_ind_value_map, 1)
 end
 
 function calculate_ind_values(
-  im::IndexMap, xs::Vector, dimensions::Vector{Int}; print_x=false
+  imap::IndexMap, xs::Vector, dimensions::Vector{Int}; print_x=false
 )
   @assert length(xs) == length(dimensions)
   ind_to_ind_value_map = Dictionary()
   for (i, x) in enumerate(xs)
     dimension = dimensions[i]
     x_rn = x
-    indices = dimension_inds(im, dimension)
-    sorted_inds = sort(indices; by=indices -> digit(im, indices))
+    indices = dimension_inds(imap, dimension)
+    sorted_inds = sort(indices; by=indices -> digit(imap, indices))
     for ind in sorted_inds
       ind_val = dim(ind) - 1
       ind_set = false
       while (!ind_set)
-        if x_rn >= index_value_to_scalar(im, ind, ind_val)
+        if x_rn >= index_value_to_scalar(imap, ind, ind_val)
           set!(ind_to_ind_value_map, ind, ind_val)
-          x_rn -= index_value_to_scalar(im, ind, ind_val)
+          x_rn -= index_value_to_scalar(imap, ind, ind_val)
           ind_set = true
         else
           ind_val = ind_val - 1
@@ -133,7 +133,7 @@ function calculate_ind_values(
     end
 
     if print_x
-      x_bitstring = calculate_x(im, ind_to_ind_value_map, dimension)
+      x_bitstring = calculate_x(imap, ind_to_ind_value_map, dimension)
       println(
         "Dimension $dimension. Actual value of x is $x but bitstring rep. is $x_bitstring"
       )
@@ -142,23 +142,31 @@ function calculate_ind_values(
   return ind_to_ind_value_map
 end
 
-function calculate_ind_values(im::IndexMap, x::Number, dimension::Int; kwargs...)
-  return calculate_ind_values(im, [x], [dimension]; kwargs...)
+function calculate_ind_values(imap::IndexMap, x::Number, dimension::Int; kwargs...)
+  return calculate_ind_values(imap, [x], [dimension]; kwargs...)
 end
-function calculate_ind_values(im::IndexMap, xs::Vector; kwargs...)
-  return calculate_ind_values(im, xs, [i for i in 1:length(xs)]; kwargs...)
+function calculate_ind_values(imap::IndexMap, xs::Vector; kwargs...)
+  return calculate_ind_values(imap, xs, [i for i in 1:length(xs)]; kwargs...)
 end
-function calculate_ind_values(im::IndexMap, x::Number; kwargs...)
-  return calculate_ind_values(im, [x], [1]; kwargs...)
+function calculate_ind_values(imap::IndexMap, x::Number; kwargs...)
+  return calculate_ind_values(imap, [x], [1]; kwargs...)
 end
 
-function grid_points(im::IndexMap, N::Int, dimension::Int)
-  dims = dim.(dimension_inds(im, dimension))
+function grid_points(imap::IndexMap, N::Int, dimension::Int)
+  dims = dim.(dimension_inds(imap, dimension))
   @assert all(y -> y == first(dims), dims)
   base = first(dims)
-  vals = Vector
-  L = length(dimension_inds(im, dimension))
+  L = length(dimension_inds(imap, dimension))
   a = round(base^L / N)
   grid_points = [i * (a / base^L) for i in 0:(N + 1)]
   return filter(x -> x <= 1, grid_points)
 end
+
+" Obtain all grid points of a given dimension"
+function grid_points(imap::IndexMap, dimension::Int)
+  dims = dim.(dimension_inds(imap, dimension))
+  @assert all(y -> y == first(dims), dims)
+  base = dims[dimension]
+  L = length(dimension_inds(imap, dimension))
+  return grid_points(imap, base^L, dimension)
+end