Skip to content

Commit

Permalink
Merge pull request #22 from JoeyT1994/fix
Browse files Browse the repository at this point in the history
Small fixes and new features
  • Loading branch information
JoeyT1994 authored Apr 25, 2024
2 parents 4103432 + 735d3c9 commit bbdec04
Show file tree
Hide file tree
Showing 7 changed files with 167 additions and 114 deletions.
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
4 changes: 3 additions & 1 deletion src/ITensorNumericalAnalysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
17 changes: 17 additions & 0 deletions src/elementary_functions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
100 changes: 37 additions & 63 deletions src/elementary_operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
104 changes: 56 additions & 48 deletions src/indexmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand All @@ -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"
)
Expand All @@ -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
2 changes: 1 addition & 1 deletion src/itensornetworkfunction.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
using Base: Base
using ITensorNetworks:
ITensorNetworks, AbstractITensorNetwork, data_graph, data_graph_type, scalar
using ITensors: ITensor, dim, contract, siteinds, onehot
using ITensors: ITensor, dim, contract, siteinds, onehot, maxlinkdim
using Graphs: Graphs

default_contraction_alg() = "bp"
Expand Down
Loading

2 comments on commit bbdec04

@JoeyT1994
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request updated: JuliaRegistries/General/105552

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.0 -m "<description of version>" bbdec04d7f152a46e2bbada896220ddc95d4c703
git push origin v0.1.0

Please sign in to comment.