Fix scalar contract bug and other bug fixes and debugging tools (#547)

* Bump NDTensors to v0.1.21 and ITensors to v0.1.32.

* Add regression test for scalar contract bug with mixed element types.

* Add `@debug_check` macro, turn on and off with `enable_debug_checks()` and `disable_debug_checks()`.

* Add `index_id_rng()`, an RNG specifically for generating Index IDs.

* Add check for arrow directions in QN IndexSet contraction.

NEWS.md

@@ -1,3 +1,11 @@
+ITensors v0.1.32 Release Notes
+==============================
+* Update to NDTensors v0.1.21, which includes a bug fix for scalar-like tensor contractions involving mixed element types (NDTensors PR #58).
+* Docs for observer system and DMRGObserver (PR #546).
+* Add `ITensors.@debug_check`, `ITensors.enable_debug_checks()`, and `ITensors.disable_debug_checks()` for denoting that a block of code is a debug check, and turning on and off debug checks (off by default). Use to check for repeated indices in IndexSet construction and other checks (PR #547).
+* Add `index_id_rng()`, an RNG specifically for generating Index IDs. Set the seed with `Random.seed!(index_id_rng(), 1234)`. This makes the random stream of number seperate for Index IDs and random elements, and helps avoid Index ID clashes with reading and writing (PR #547).
+* Add back checking for proper QN Index directions in contraction (PR #547).
+
 ITensors v0.1.31 Release Notes
 ==============================
 * Update to NDTensors v0.1.20, which includes some more general block sparse slicing operations as well as optimizations for contracting scalar-like (length 1) tensors (NDTensors PR #57).

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensors"
 uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 authors = ["mfishman <mfishman@caltech.edu>"]
-version = "0.1.31"
+version = "0.1.32"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -19,7 +19,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Compat = "2.1, 3"
 HDF5 = "0.13.1, 0.14"
 KrylovKit = "0.4.2, 0.5"
-NDTensors = "= 0.1.20"
+NDTensors = "= 0.1.21"
 PackageCompiler = "1.0.0"
 StaticArrays = "0.12.3"
 TimerOutputs = "0.5.5"

src/ITensors.jl

@@ -44,6 +44,11 @@ include("imports.jl")
 #
 include("global_variables.jl")
 
+#####################################
+# Debug checking
+#
+include("debug_checks.jl")
+
 #####################################
 # Index and IndexSet
 #
@@ -120,6 +125,10 @@ include("packagecompile/compile.jl")
 #
 include("developer_tools.jl")
 
+function __init__()
+  resize!(empty!(INDEX_ID_RNGs), Threads.nthreads()) # ensures that we didn't save a bad object
+end
+
 #####################################
 # Precompile certain functions
 # (generated from precompile/make_precompile.jl

src/debug_checks.jl

@@ -0,0 +1,24 @@
+
+# XXX: rename use_debug_checks
+use_debug_checks() = false
+
+macro debug_check(ex)
+  quote
+    if use_debug_checks()
+      $(esc(ex))
+    end
+  end
+end
+
+function enable_debug_checks()
+  if !getfield(@__MODULE__, :use_debug_checks)()
+    Core.eval(@__MODULE__, :(use_debug_checks() = true))
+  end
+end
+
+function disable_debug_checks()
+  if getfield(@__MODULE__, :use_debug_checks)()
+    Core.eval(@__MODULE__, :(use_debug_checks() = false))
+  end
+end
+

src/decomp.jl

@@ -149,7 +149,7 @@ iterate(E::TruncEigen, ::Val{:r}) = (E.r, Val(:done))
 iterate(E::TruncEigen, ::Val{:done}) = nothing
 
 function eigen(A::ITensor{N}, Linds, Rinds; kwargs...) where {N}
-  @debug begin
+  @debug_check begin
     if hasqns(A)
       @assert flux(A) == QN()
     end
@@ -240,7 +240,7 @@ function eigen(A::ITensor{N}, Linds, Rinds; kwargs...) where {N}
   # The right eigenvectors, after being applied to A
   Vt = replaceinds(V, (Ris..., r), (Lis..., l))
 
-  @debug begin
+  @debug_check begin
     if hasqns(A)
       @assert flux(D) == QN()
       @assert flux(V) == QN()

src/exports.jl

@@ -20,6 +20,9 @@ export
   # Deprecated
   addblock!,
 
+# ITensors.jl
+  index_id_rng,
+
 # argsdict/argsdict.jl
   argsdict,
 

src/index.jl

@@ -1,6 +1,24 @@
 
+#const IDType = UInt128
 const IDType = UInt64
 
+# Custom RNG for Index id
+# Vector of RNGs, one for each thread
+const INDEX_ID_RNGs = MersenneTwister[]
+@inline index_id_rng() = index_id_rng(Threads.threadid())
+@noinline function index_id_rng(tid::Int)
+  0 < tid <= length(INDEX_ID_RNGs) || _index_id_rng_length_assert()
+  if @inbounds isassigned(INDEX_ID_RNGs, tid)
+    @inbounds MT = INDEX_ID_RNGs[tid]
+  else
+    MT = MersenneTwister()
+    @inbounds INDEX_ID_RNGs[tid] = MT
+  end
+  return MT
+end
+@noinline _index_id_rng_length_assert() =  @assert false "0 < tid <= length(INDEX_ID_RNGs)"
+
+
 """
 An `Index` represents a single tensor index with fixed dimension `dim`. Copies of an Index compare equal unless their 
 `tags` are different.
@@ -54,7 +72,7 @@ julia> tags(i)
 ```
 """
 function Index(dim::Int; tags="", plev=0)
-  return Index(rand(IDType), dim, Neither, tags, plev)
+  return Index(rand(index_id_rng(), IDType), dim, Neither, tags, plev)
 end
 
 """
@@ -182,7 +200,7 @@ Produces an `Index` with the same properties (dimension or QN structure)
 but with a new `id`.
 """
 sim(i::Index; tags = copy(tags(i)), plev = plev(i), dir = dir(i)) =
-  Index(rand(IDType), copy(space(i)), dir, tags, plev)
+  Index(rand(index_id_rng(), IDType), copy(space(i)), dir, tags, plev)
 
 """
     dag(i::Index)

src/indexset.jl

@@ -4,6 +4,11 @@ struct IndexSet{N, IndexT <: Index, DataT <: Tuple}
   data::DataT
 
   function IndexSet{N, IndexT, DataT}(data) where {N, IndexT, DataT <: NTuple{N, IndexT}}
+    @debug_check begin
+      if !allunique(data)
+        error("Trying to create IndexSet with collection of indices $data. Indices must be unique.")
+      end
+    end
     return new{N, IndexT, DataT}(data)
   end
 
@@ -973,15 +978,16 @@ hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
 function replaceinds(is::IndexSet, inds1, inds2)
   is1 = IndexSet(inds1)
   poss = indexin(is1, is)
+  is_tuple = Tuple(is)
   for (j, pos) in enumerate(poss)
     isnothing(pos) && continue
-    i1 = is[pos]
+    i1 = is_tuple[pos]
     i2 = inds2[j]
     i2 = setdir(i2, dir(i1))
     space(i1) ≠ space(i2) && error("Indices must have the same spaces to be replaced")
-    is = setindex(is, i2, pos)
+    is_tuple = setindex(is_tuple, i2, pos)
   end
-  return is
+  return IndexSet(is_tuple)
 end
 
 replaceind(is::IndexSet, i1::Index, i2::Index) =
@@ -1008,12 +1014,18 @@ end
 
 function compute_contraction_labels(Ais::IndexSet{NA},
                                     Bis::IndexSet{NB}) where {NA,NB}
+  have_qns = hasqns(Ais) && hasqns(Bis)
   Alabels = MVector{NA,Int}(ntuple(_->0,Val(NA)))
   Blabels = MVector{NB,Int}(ntuple(_->0,Val(NB)))
 
   ncont = 0
   for i = 1:NA, j = 1:NB
-    if Ais[i] == Bis[j]
+    Ais_i = @inbounds Ais[i]
+    Bis_j = @inbounds Bis[j]
+    if Ais_i == Bis_j
+      if have_qns && (dir(Ais_i) ≠ -dir(Bis_j))
+        error("Attempting to contract IndexSet:\n$(Ais)with IndexSet:\n$(Bis)QN indices must have opposite direction to contract.")
+      end
       Alabels[i] = Blabels[j] = -(1+ncont)
       ncont += 1
     end

src/itensor.jl

@@ -1448,7 +1448,7 @@ computed internally.
 function exp(A::ITensor{N}, Linds, Rinds; kwargs...) where {N}
   ishermitian=get(kwargs,:ishermitian,false)
 
-  @debug begin
+  @debug_check begin
     if hasqns(A)
       @assert flux(A) == QN()
     end
@@ -1791,7 +1791,7 @@ If the ITensor is empty or it has no QNs, returns `nothing`.
 """
 function flux(T::ITensor)
   (!hasqns(T) || isempty(T)) && return nothing
-  @debug checkflux(T)
+  @debug_check checkflux(T)
   block1 = first(eachnzblock(T))
   return flux(T, block1)
 end

src/mps/abstractmps.jl

@@ -351,6 +351,7 @@ end
 
 """
     firstsiteind(M::Union{MPS,MPO}, j::Integer; kwargs...)
+    siteind(::typeof(first), M::Union{MPS,MPO}, j::Integer; kwargs...)
 
 Return the first site Index found on the MPS or MPO
 (the first Index unique to the `j`th MPS/MPO tensor).
@@ -372,6 +373,9 @@ function firstsiteind(M::AbstractMPS, j::Integer;
   return si
 end
 
+siteind(::typeof(first), M::AbstractMPS, j::Integer; kwargs...) =
+  firstsiteind(M, j; kwargs...)
+
 """
     siteinds(M::Union{MPS, MPO}}, j::Integer; kwargs...)
 

src/mps/dmrg.jl

@@ -94,8 +94,9 @@ end
 
 
 function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
-  # Debug level checks
-  @debug begin
+  @debug_check begin
+    # Debug level checks
+    # Enable with ITensors.enable_debug_checks()
     checkflux(psi0)
     checkflux(PH)
   end
@@ -153,7 +154,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
 
     for (b, ha) in sweepnext(N)
 
-      @debug begin
+      @debug_check begin
         checkflux(psi)
         checkflux(PH)
       end
@@ -162,7 +163,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
       position!(PH, psi, b)
       end
 
-      @debug begin
+      @debug_check begin
         checkflux(psi)
         checkflux(PH)
       end
@@ -188,7 +189,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
         drho = noise(sweeps, sw) * noiseterm(PH,phi,ortho)
       end
 
-      @debug begin
+      @debug_check begin
         checkflux(phi)
       end
 
@@ -203,7 +204,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
                                        svd_alg = svd_alg)
       end
 
-      @debug begin
+      @debug_check begin
         checkflux(psi)
         checkflux(PH)
       end

src/mps/mps.jl

@@ -329,9 +329,22 @@ productMPS(sites::Vector{ <: Index},
 """
     siteind(M::MPS, j::Int; kwargs...)
 
-Get the site Index of the MPS.
+Get the first site Index of the MPS. Return `nothing` if none is found.
 """
-siteind(M::MPS, j::Int; kwargs...) = firstsiteind(M, j; kwargs...)
+siteind(M::MPS, j::Int; kwargs...) = siteind(first, M, j; kwargs...)
+
+"""
+    siteind(::typeof(only), M::MPS, j::Int; kwargs...)
+
+Get the only site Index of the MPS. Return `nothing` if none is found.
+"""
+function siteind(::typeof(only), M::MPS, j::Int; kwargs...)
+  is = siteinds(M, j; kwargs...)
+  if isempty(is)
+    return nothing
+  end
+  return only(siteinds(M, j; kwargs...))
+end
 
 """
     siteinds(M::MPS)
@@ -342,7 +355,7 @@ siteinds(M::MPS) = [siteind(M, j) for j in 1:length(M)]
 
 function replace_siteinds!(M::MPS, sites)
   for j in eachindex(M)
-    sj = siteind(M, j)
+    sj = siteind(only, M, j)
     replaceind!(M[j], sj, sites[j])
   end
   return M

src/qn/qnindex.jl

@@ -93,7 +93,7 @@ function Index(qnblocks::QNBlocks; dir::Arrow = Out, tags = "", plev = 0)
   # TODO: make this a debug check?
   #have_same_qns(qnblocks) || error("When creating a QN Index, the QN blocks must have the same QNs")
   #have_same_mods(qnblocks) || error("When creating a QN Index, the QN blocks must have the same mods")
-  return Index(rand(IDType), qnblocks, dir, tags, plev)
+  return Index(rand(index_id_rng(), IDType), qnblocks, dir, tags, plev)
 end
 
 """

test/debug_checks.jl

@@ -0,0 +1,26 @@
+using ITensors
+using Test
+
+@testset "Test debug checks on IndexSet construction" begin
+  i = Index(2, "i")
+
+  @test !ITensors.use_debug_checks()
+  # Test that no error is thrown in constructor
+  @test IndexSet(i, i) isa IndexSet{2}
+  @test IndexSet(i, i') isa IndexSet{2}
+
+  # Turn on debug checks
+  ITensors.enable_debug_checks()
+  @test ITensors.use_debug_checks()
+  @test_throws ErrorException IndexSet(i, i)
+  # Test that no error is thrown in constructor
+  @test IndexSet(i, i') isa IndexSet{2}
+
+  # Turn off debug checks
+  ITensors.disable_debug_checks()
+  @test !ITensors.use_debug_checks()
+  # Test that no error is thrown in constructor
+  @test IndexSet(i, i) isa IndexSet{2}
+  @test IndexSet(i, i') isa IndexSet{2}
+end
+