Zeroing out operator

examples/2d_laplace_solver.jl

@@ -3,6 +3,7 @@ using ITensorNumericalAnalysis
 
 using Graphs: SimpleGraph, uniform_tree
 using NamedGraphs: NamedGraph
+using NamedGraphs.GraphsExtensions: rename_vertices
 using ITensors: ITensors, Index, siteinds, dim, tags, replaceprime!, MPO, MPS, inner
 using ITensorNetworks: ITensorNetwork, dmrg, ttn, maxlinkdim
 using Dictionaries: Dictionary

examples/boundary_conditions.jl

@@ -0,0 +1,66 @@
+using Test
+using ITensorNumericalAnalysis
+
+using NamedGraphs.NamedGraphGenerators: named_comb_tree
+using ITensors: ITensors, inner
+using ITensorNetworks: ITensorNetwork, ttn, maxlinkdim
+using Random: seed!
+using ITensors: Ops
+
+using UnicodePlots
+
+seed!(1234)
+L = 12
+g = named_comb_tree((2, L ÷ 2))
+lastDigit = 1 - 1 / 2^(L ÷ 2)
+
+s = continuous_siteinds(g; map_dimension=2)
+
+ψ_fxy = cos_itn(s; k=π)
+#ψ_fxy = const_itn(s; c=3, linkdim=3) # note if you use const, need big linkdim
+@show maxlinkdim(ψ_fxy)
+
+maxdim = 10
+cutoff = 1e-6
+@show cutoff
+ϕ_fxy = map_to_zeros(ψ_fxy, [0, lastDigit, 0, lastDigit], [1, 1, 2, 2]; cutoff, maxdim)
+@show maxlinkdim(ϕ_fxy)
+
+n_grid = 100
+x_vals, y_vals = grid_points(s, n_grid, 1)[1:2:(end - 1)],
+grid_points(s, n_grid, 2)[1:2:(end - 1)]
+# fix for if we don't include all 1s
+if x_vals[end] != lastDigit
+  push!(x_vals, lastDigit)
+end
+if y_vals[end] != lastDigit
+  push!(y_vals, lastDigit)
+end
+vals = zeros((length(x_vals), length(y_vals)))
+for (i, x) in enumerate(x_vals)
+  for (j, y) in enumerate(y_vals)
+    vals[i, j] = real(calculate_fxyz(ϕ_fxy, [x, y]; alg="exact"))
+  end
+end
+
+println("Here is the heatmap of the 2D function")
+display(heatmap(vals; xfact=1 / 32, yfact=1 / 32, xoffset=0, yoffset=0, colormap=:inferno))
+
+y = 0.5
+vals2 = zeros(length(x_vals))
+for (i, x) in enumerate(x_vals)
+  vals2[i] = real(calculate_fxyz(ϕ_fxy, [x, y]; alg="exact"))
+end
+
+lp = lineplot(x_vals, vals2; name="cut y=$y")
+
+x = 0.5
+vals3 = zeros(length(y_vals))
+for (i, y) in enumerate(y_vals)
+  vals3[i] = real(calculate_fxyz(ϕ_fxy, [x, y]; alg="exact"))
+end
+
+println("Here is a cut of the function at x = $x or y = $y")
+display(lineplot!(lp, y_vals, vals3; name="cut x=$x"))
+
+@show vals2[1], vals2[end], vals3[1], vals3[end]

src/ITensorNumericalAnalysis.jl

@@ -45,7 +45,10 @@ export const_itensornetwork,
   fourth_derivative_operator,
   identity_operator,
   delta_x,
-  delta_xyz
+  delta_xyz,
+  map_to_zero_operator,
+  map_to_zeros,
+  const_plane_op
 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

@@ -195,23 +195,103 @@ 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))
+"Create a product state of a given bit configuration. Will make planes if all dims not specificed"
+function delta_xyz(
+  s::IndsNetworkMap,
+  xs::Vector{<:Number},
+  dims::Vector{Int}=[i for i in 1:length(xs)];
+  kwargs...,
+)
+  ivmap = calculate_ind_values(s, xs, dims)
+  tn = ITensorNetwork(
+    v -> only(s[v]) in keys(ivmap) ? string(ivmap[only(s[v])]) : ones(dim(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
 
+function delta_xyz(
+  s::IndsNetworkMap,
+  points::Vector{<:Vector},
+  points_dims::Vector{<:Vector}=[[i for i in 1:length(xs)] for xs in points];
+  kwargs...,
+)
+  @assert length(points) != 0
+  @assert length(points) == length(points_dims)
+  ψ = reduce(
+    +, [delta_xyz(s, xs, dims; kwargs...) for (xs, dims) in zip(points, points_dims)]
+  )
+  return ψ
+end
+
+" Function to manipulate delta functions. Defaults to map_to_zero behavior"
+function delta_kernel(
+  s::IndsNetworkMap,
+  points::Vector{<:Vector},
+  points_dims::Vector{<:Vector}=[[i for i in 1:length(xs)] for xs in points];
+  remove_overlap=true,
+  coeff::Number=-1,
+  include_identity=true,
+  truncate_kwargs...,
+)
+  ψ = coeff * delta_xyz(s, points, points_dims; truncate_kwargs...)
+
+  if include_identity
+    ψ = const_itn(s) + ψ
+  end
+
+  if remove_overlap && length(points) > 1
+    overlap_points, overlap_dims = Vector{Vector}(), Vector{Vector}()
+    # determine intersection of any points, 
+    # and remove them with the opposite sign
+    for i in 1:length(points)
+      p1, d1 = points[i], points_dims[i]
+      for j in (i + 1):length(points)
+        p2, d2 = points[j], points_dims[j]
+
+        # same dimensions, and no point overlap,
+        # can safely ignore
+        (all(d1 .≈ d2) && !all(p1 .≈ p2)) && continue
+
+        # check if dims are the same. 
+        # If they are, check the corresponding dim
+        ps_ = [p1; p2]
+        ds_ = [d1; d2]
+        order = sortperm(ds_)
+        ps, ds = [ps_[order[1]]], [ds_[order[1]]]
+        for k in 2:length(ds_)
+          if (ds_[order[k]] != ds_[order[k - 1]])
+            push!(ps, ps_[order[k]])
+            push!(ds, ds_[order[k]])
+            continue
+          end
+          # found two matching elements
+          if ps_[k] ≈ ps_[k - 1]
+            continue # added previously
+          else # there's no overap here, continue
+            ps, ds = [], []
+            break
+          end
+        end
+        #(length(Set(ds)) != length(ds)) && continue
+        (length(ds) == 0) && continue
+        push!(overlap_points, Vector(ps))
+        push!(overlap_dims, Vector(ds))
+      end
+    end
+    if length(overlap_points) != 0
+      ψ = ψ + -coeff * delta_xyz(s, overlap_points, overlap_dims; truncate_kwargs...)
+    end
+  end
+
+  return ψ
+end
 const const_itn = const_itensornetwork
 const poly_itn = polynomial_itensornetwork
 const cosh_itn = cosh_itensornetwork

src/elementary_operators.jl

@@ -7,6 +7,7 @@ using ITensors:
   noprime,
   op,
   Op,
+  Ops,
   truncate,
   replaceinds,
   delta,
@@ -191,6 +192,56 @@ function identity_operator(s::IndsNetworkMap; kwargs...)
   return ITensorNetwork(Op("I"), operator_inds)
 end
 
+" Create an operator bitstring corresponding to the number x"
+function point_to_opsum(s::IndsNetworkMap, x::Number, dim::Int)
+  ttn_op = OpSum()
+  ind_to_ind_value_map = calculate_ind_values(s, x, dim)
+  string_sites = [
+    (ind_to_ind_value_map[only(s[v])] == 1) ? ("Dup", v) : ("Ddn", v) for
+    v in dimension_vertices(s, dim)
+  ]
+  add!(ttn_op, 1.0, (string_sites...)...)
+  return ttn_op
+end
+
+" Create an operator which maps a function to 0 at all points in xs"
+function map_to_zero_operator(
+  s::IndsNetworkMap, xs::Vector, dims::Vector=[1 for _ in xs]; truncate_kwargs...
+)
+  return operator_proj(
+    delta_kernel(
+      s,
+      [[x] for x in xs],
+      [[dim] for dim in dims];
+      remove_overlap=true,
+      coeff=-1,
+      include_identity=true,
+      truncate_kwargs...,
+    ),
+  )
+end
+
+function map_to_zero_operator(s::IndsNetworkMap, x::Number, dim::Int=1; truncate_kwargs...)
+  return map_to_zero_operator(s, [x], [dim]; truncate_kwargs...)
+end
+
+" Map the points xs in dimension dims of the function f to 0"
+function map_to_zeros(
+  f::ITensorNetworkFunction,
+  xs::Vector,
+  dims::Vector=[1 for _ in xs];
+  truncate_kwargs=(;), # for map_operator
+  kwargs..., # for operate
+)
+  s = indsnetworkmap(f)
+  zero_op = map_to_zero_operator(s, xs, dims; truncate_kwargs...)
+  return operate(zero_op, f; kwargs...)
+end
+
+function map_to_zeros(f::ITensorNetworkFunction, x::Number, dim::Int=1; kwargs...)
+  return map_to_zeros(f, [x], [dim]; kwargs...)
+end
+
 " Take |f> and create an operator |f><δ| "
 function operator_proj(fx::ITensorNetworkFunction)
   fx = copy(fx)

src/indsnetworkmap.jl

@@ -94,8 +94,12 @@ function vertex_digit(inm::IndsNetworkMap, v)
   return digit(inm, only(inds(inm, v)))
 end
 
-function dimension_vertices(inm::IndsNetworkMap, dimension::Int)
-  return filter(v -> vertex_dimension(inm, v) == dimension, vertices(inm))
+function dimension_vertices(inm::IndsNetworkMap, dim::Int)
+  return dimension_vertices(inm, [dim])
+end
+
+function dimension_vertices(inm::IndsNetworkMap, dims::Vector{Int})
+  return filter(v -> vertex_dimension(inm, v) in dims, vertices(inm))
 end
 
 function vertex(inm::IndsNetworkMap, dimension::Int, digit::Int)

test/test_examples.jl

@@ -3,7 +3,9 @@ using ITensorNumericalAnalysis: ITensorNumericalAnalysis
 using Test: @testset
 
 @testset "Test examples" begin
-  example_files = ["2d_laplace_solver.jl", "construct_multi_dimensional_function.jl"]
+  example_files = [
+    "2d_laplace_solver.jl", "construct_multi_dimensional_function.jl", "fredholm_solver.jl"
+  ]
   @testset "Test $example_file" for example_file in example_files
     include(joinpath(pkgdir(ITensorNumericalAnalysis), "examples", example_file))
   end

test/test_itensorfunction.jl

@@ -206,4 +206,54 @@ Random.seed!(1234)
       @test c * tanh(k * x + a) + c * tanh(k * y + a) ≈ fxy_xy
     end
   end
+  @testset "test delta_xyz" begin
+    L = 10
+    g = named_grid((L, 1))
+    s = continuous_siteinds(g; map_dimension=2)
+    x0, y0 = 0.625, 0.25
+    delta = 2.0^(-1.0 * L)
+    lastDigit = 1 - delta
+    xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, lastDigit]
+    @testset "test single point" begin
+      ψ = delta_xyz(s, [x0, y0])
+      @test calculate_fxyz(ψ, [x0, y0], [1, 2]) ≈ 1
+      # test another point
+      @test calculate_fxyz(ψ, [y0, x0], [1, 2]) ≈ 0
+    end
+    @testset "test plane" begin
+      ψ = delta_xyz(s, [y0], [2])
+
+      # should be 1 in the plane
+      for x in xs
+        @test calculate_fxyz(ψ, [x, y0], [1, 2]) ≈ 1
+      end
+      # test random points
+      for x in xs
+        @test calculate_fxyz(ψ, [x, 0.5], [1, 2]) ≈ 0
+      end
+    end
+    @testset "test sums of points" begin
+      points = [[x0, y0], [y0, x0]]
+      ψ = delta_xyz(s, points)
+      @test calculate_fxyz(ψ, [x0, y0], [1, 2]) ≈ 1
+      @test calculate_fxyz(ψ, [y0, x0], [1, 2]) ≈ 1
+      # test other points
+      @test calculate_fxyz(ψ, [0, 0], [1, 2]) ≈ 0
+      @test calculate_fxyz(ψ, [0, y0], [1, 2]) ≈ 0
+    end
+
+    @testset "test sums of points and plane" begin
+      p0 = 0.5
+      points = [[x0, y0], [p0]]
+      dims = [[1, 2], [2]]
+      ψ = delta_xyz(s, points, dims)
+      @test calculate_fxyz(ψ, [x0, y0], [1, 2]) ≈ 1
+      for x in xs
+        @test calculate_fxyz(ψ, [x, p0], [1, 2]) ≈ 1
+      end
+      ## test other points
+      @test calculate_fxyz(ψ, [0, 0], [1, 2]) ≈ 0
+      @test calculate_fxyz(ψ, [0, y0], [1, 2]) ≈ 0
+    end
+  end
 end