Account for BP edge case where network evaluates to 0 (#178)

Author: JoeyT1994

src/caches/beliefpropagationcache.jl

@@ -19,7 +19,10 @@ default_message_norm(m::ITensor) = norm(m)
 function default_message_update(contract_list::Vector{ITensor}; kwargs...)
   sequence = optimal_contraction_sequence(contract_list)
   updated_messages = contract(contract_list; sequence, kwargs...)
-  updated_messages /= norm(updated_messages)
+  message_norm = norm(updated_messages)
+  if !iszero(message_norm)
+    updated_messages /= message_norm
+  end
   return ITensor[updated_messages]
 end
 @traitfn default_bp_maxiter(g::::(!IsDirected)) = is_tree(g) ? 1 : nothing

src/contract.jl

@@ -73,6 +73,7 @@ function logscalar(
     denominator_terms
   end
 
+  any(iszero, denominator_terms) && return -Inf
   return sum(log.(numerator_terms)) - sum(log.((denominator_terms)))
 end
 

test/test_belief_propagation.jl

@@ -2,7 +2,6 @@
 using Compat: Compat
 using Graphs: vertices
 # Trigger package extension.
-using GraphsFlows: GraphsFlows
 using ITensorNetworks:
   ITensorNetworks,
   BeliefPropagationCache,
@@ -18,6 +17,7 @@ using ITensorNetworks:
   message,
   partitioned_tensornetwork,
   random_tensornetwork,
+  scalar,
   siteinds,
   split_index,
   tensornetwork,
@@ -28,7 +28,7 @@ using ITensors: ITensors, ITensor, combiner, dag, inds, inner, op, prime, random
 using ITensorNetworks.ModelNetworks: ModelNetworks
 using ITensors.NDTensors: array
 using LinearAlgebra: eigvals, tr
-using NamedGraphs: NamedEdge
+using NamedGraphs: NamedEdge, NamedGraph
 using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid
 using NamedGraphs.PartitionedGraphs: PartitionVertex, partitionedges
 using Random: Random
@@ -75,5 +75,12 @@ using Test: @test, @testset
 
   @test all(eig -> imag(eig) ≈ 0, eigs)
   @test all(eig -> real(eig) >= -eps(eltype(eig)), eigs)
+
+  #Test edge case of network which evalutes to 0
+  χ = 2
+  g = named_grid((3, 1))
+  ψ = random_tensornetwork(ComplexF64, g; link_space=χ)
+  ψ[(1, 1)] = 0.0 * ψ[(1, 1)]
+  @test iszero(scalar(ψ; alg="bp"))
 end
 end