Zeroing out operator

[ryanlevy]

This PR will hopefully begin to build custom boundary conditions operator.

Currently have added:

1. zero_point_op(s::IndsNetworkMap, xs::Vector, dimension::Int) (feel free to suggest a better name) - This operator generates (1-\sum_xp delta(x-xp))f(x) which will zero out all points xp in xs

[JoeyT1994]

@ryanlevy thanks for this this will be super helpful.

I think we should change the example to a test which constructs some functions and checks these zero operators correctly map them? It could either be added to test_operators or a separate test?

[ryanlevy]

I like where you're headed with the delta functions but I still think the interface is a bit off. So I imagine the following scenarios:

• I want to make a function that is a single point (currently delta_xyz)
• I want to make a function that is a plane (currently ~ const*const_plane_op)
• I want to make a plane operator from a delta functions (O = delta(x)* 1_y * 1_z)
• I want to sum up a series of planes (O = delta(x) + delta(y)) or series of points
• I want to sum up a series of points (O= delta(r - p)+delta(r-p))

Then we have as input Vector[Vector] and the corresponding dimensions Vector[Vector]. I think that should include all of these, where the first two would be Vector[Array{Number,1}], The second Vector[Array{Number,dim}]. And could write functions that conveniently slice the above into nicer interfaces if needed. Particularly in the last scenario, we want to include if there should be an identity and a sign/coeff (i.e. 1-operators or operators). I am inclined to leave the BC function out of all this and just build in a more robust delta setup also.

Then move delta_xyz to something like running over dimension_vertices(s, dim) instead of everything

What do you think?

[JoeyT1994]

Yeah I totally agree that we should do it all by interfacing off of the delta_xyz function. Then all of those scenarios should be possible.

I think all we actually need to do is "elevate" the delta function a bit to make this work. At the moment it relies on you specifying a single point xs in the full d dimensional space where d = dim(s::IndsNetworkMap)

We could define the following

"Create prod_{x \in xs}\delta(x)"
function delta_xyz(s::IndsNetworkMap, xs::Vector, dimensions::Vector{Int}; kwargs...)
  ind_map = calculate_ind_values(s, xs, dimensions)
  tn = ITensorNetwork(v -> only(s[v]) in keys(ind_map) ? string(ind_to_ind_value_map[only(s[v])]) : [1.0 for i in 1:dim(only(s[v]))], indsnetwork(s))
  return ITensorNetworkFunction(tn, s)
end

this would create the same delta function as before if you specify a single point in the full d dimensional space but if you specify, say, xs = [x], dims = [1] and s is a 2D IndsNetworkMap then it will gives you a function which is non-zero only along a single line of the two dimensional plane.

With this function and operator_proj all of the funcionality you describe in the bullet points above can be built with very minimal code and without ever using OpSum

[ryanlevy]

Thanks to our discussion, I've implemented a generalized delta_xyz which can arbitrarily make points or lower dimensional structures (lines, planes). Then there's a new function delta_kernel, which map_to_zero/const_plane_op is now a wrapper for, which does the logic of adding/subtracting the correct delta functions to create what you want. So for example, in 3D:

Zo = ITensorNumericalAnalysis.delta_kernel(
  s,
  [[0], [0,0 ]],
  [[1], [1, 2]];
  remove_overlap=true,
  coeff=-1,
  include_identity=true,
  promote_operator=true,
)

The two delta functions are a plane at x=0 (running over y,z) and a line at (0,0) (running over z). They intersect at a point which is properly removed. Likewise for [[0],[0]]], [[3],[1]], planes at z=0 and x=0, the intersecting line is taken care of. TODO, maybe in the future is to make the coeff more general for heterogeneous BC conditions