Update map_eigvals to use newer pattern and make more robust

Author: emstoudenmire

src/lib/ITensorsExtensions/src/itensor.jl

@@ -5,13 +5,15 @@ using ITensors:
   Index,
   commonind,
   dag,
+  dir,
   hasqns,
   inds,
   isdiag,
   itensor,
   map_diag,
   noncommonind,
   noprime,
+  permute,
   replaceind,
   replaceinds,
   sim,
@@ -63,17 +65,20 @@ function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...)
 end
 
 function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...)
+  Linds = isa(Linds, Index) ? [Linds] : collect(Linds)
+  Rinds = isa(Rinds, Index) ? [Rinds] : collect(Rinds)
 
   # <fermions>
-  auto_fermion_enabled = ITensors.using_auto_fermion()
-  if auto_fermion_enabled
+  fermionic_itensor =
+    ITensors.using_auto_fermion() && ITensors.has_fermionic_subspaces(inds(A))
+  if fermionic_itensor
     # If fermionic, bring indices into i',j',..,dag(j),dag(i)
     # ordering with Out indices coming before In indices
     # Resulting tensor acts like a normal matrix (no extra signs
     # when taking powers A^n)
-    if all(j->dir(j)==Out, Linds) && all(j->dir(j)==In, Rinds)
+    if all(j->dir(j)==ITensors.Out, Linds) && all(j->dir(j)==ITensors.In, Rinds)
       ordered_inds = [Linds..., reverse(Rinds)...]
-    elseif all(j->dir(j)==Out, Rinds) && all(j->dir(j)==In, Linds)
+    elseif all(j->dir(j)==ITensors.Out, Rinds) && all(j->dir(j)==ITensors.In, Linds)
       ordered_inds = [Rinds..., reverse(Linds)...]
     else
       error(
@@ -96,7 +101,7 @@ function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...)
   end
 
   # <fermions>
-  if auto_fermion_enabled
+  if fermionic_itensor
     # Ensure expA indices in "matrix" form before re-enabling fermion system
     mapped_A = permute(mapped_A, ordered_inds)
     ITensors.enable_auto_fermion()

test/test_itensorsextensions.jl

@@ -79,6 +79,7 @@ using Test: @test, @testset
   end
 
   @testset "Fermionic eigendecomp" begin
+    ITensors.enable_auto_fermion()
     s1 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=1")
     s2 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=2")
 
@@ -90,9 +91,11 @@ using Test: @test, @testset
     Ul, D, Ur = eigendecomp(T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true)
 
     @test Ul*D*Ur ≈ T
+    ITensors.disable_auto_fermion()
   end
 
   @testset "Fermionic map eigvals tests" begin
+    ITensors.enable_auto_fermion()
     s1 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=1")
     s2 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=2")
 
@@ -125,6 +128,15 @@ using Test: @test, @testset
     @test T ≈ apply(sqrtT, sqrtT)
     sqrtT = map_eigvals(sqrt, T, [dag(s1), dag(s2)], [s1', s2'], ; ishermitian=true)
     @test T ≈ apply(sqrtT, sqrtT)
+
+    # Test bosonic index case while fermion system is enabled
+    b = Index([QN("Nb", 0)=>2, QN("Nb", 1)=>2])
+    T = random_itensor(b', dag(b))
+    T = apply(T, swapprime(dag(T), 0=>1))
+    sqrtT = map_eigvals(sqrt, T, b', dag(b); ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+
+    ITensors.disable_auto_fermion()
   end
 end
 end