Use NDTensors.with_auto_fermion in ITensorsExtensions and tests (#340)

mtfishman · claude · web-flow · commit 629b79cfbaaa · 2026-04-26T16:07:30.000Z
## Summary

Adopt the new exception-safe scope wrapper `NDTensors.with_auto_fermion`
(added in `NDTensors v0.4.25`), replacing the manual
`ITensors.enable_auto_fermion()` / `disable_auto_fermion()` save-restore
patterns in `ITensorsExtensions.map_eigvals` and across four test files.

The wrapper restores the auto-fermion state via `try` / `finally`,
fixing exception-safety bugs at callsites where a thrown exception would
have left the flag flipped, and also fixing one logic bug in a test that
always re-enabled the flag instead of restoring its prior state.

## Changes

### Source

- `src/lib/ITensorsExtensions/src/itensor.jl::map_eigvals`: replace the
manual disable/enable around the eigendecomposition with a single
`NDTensors.with_auto_fermion(!fermionic_itensor) do … end` do-block.
Numerical output is bit-identical across all four (UAF ×
has-fermionic-subspaces) cases on the existing test suite plus a Case-3
differential test.

### Tests

- `test/test_opsum_to_ttn.jl` ("Multiple onsite terms"): wrap with
`with_auto_fermion()`. Also fixes a real bug in the previous
save-restore: the conditional only re-enabled the flag if it was already
enabled at entry, so a run with the flag originally off would leave it
flipped on.
- `test/test_itensorsextensions.jl` ("Fermionic eigendecomp", "Fermionic
map eigvals tests"): wrap each fermionic testset with
`with_auto_fermion()`.
- `test/test_opsum_to_ttn_mpo_cross_check.jl`: wrap "OpSum to TTN
Fermions" with `with_auto_fermion()`. Remove a vestigial save/restore
from "OpSum to TTN" (the testset never actually toggled the flag, so the
captured `auto_fermion_enabled` and the conditional restore were dead
code).
- `test/test_ttn_position.jl` ("ProjTTN position"): fold the `if
use_qns; enable; else; disable; end` conditional into
`with_auto_fermion(use_qns) do … end`.

### Version

- Patch bump `0.19.0` → `0.19.1`.

## Out of scope

The four read-only `using_auto_fermion()` queries that gate logic in
package source (`ITensorsExtensions/itensor.jl:63`,
`opsum_to_ttn.jl:225/382/581`) are state queries, not state changes;
converting them would require larger algorithm-level changes (threading
a fermion-mode parameter through call signatures). Left unchanged here.

---------

Co-authored-by: Claude Opus 4.7 (1M context) &lt;noreply@anthropic.com&gt;
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
-version = "0.19.0"
+version = "0.19.1"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>, Joseph Tindall <jtindall@flatironinstitute.org> and contributors"]
 
 [workspace]
diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl
@@ -64,21 +64,18 @@ function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...)
     if fermionic_itensor
         A = make_bosonic(A::ITensor, Linds, Rinds)
         ordered_inds = inds(A)
-        ITensors.disable_auto_fermion()
     end
 
-    if isdiag(A)
-        mapped_A = map_diag(f, A)
-    else
+    mapped_A = NDTensors.with_auto_fermion(!fermionic_itensor) do
+        isdiag(A) && return map_diag(f, A)
         Ul, D, Ur = eigendecomp(A, Linds, Rinds; kws...)
-        mapped_A = Ul * map_diag(f, D) * Ur
+        return Ul * map_diag(f, D) * Ur
     end
 
     # <fermions>
     if fermionic_itensor
         # Ensure indices in "matrix" form before re-enabling fermion system
         mapped_A = permute(mapped_A, ordered_inds)
-        ITensors.enable_auto_fermion()
     end
 
     return mapped_A
diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl
@@ -1,7 +1,8 @@
 @eval module $(gensym())
 using ITensorNetworks.ITensorsExtensions: eigendecomp, map_eigvals
-using ITensors: ITensors, ITensor, Index, QN, apply, dag, delta, inds, mapprime, noprime,
-    norm, op, permute, prime, random_itensor, replaceind, replaceinds, sim, swapprime
+using ITensors.NDTensors: with_auto_fermion
+using ITensors: ITensor, Index, QN, apply, dag, delta, inds, mapprime, noprime, norm, op,
+    permute, prime, random_itensor, replaceind, replaceinds, sim, swapprime
 using StableRNGs: StableRNG
 using Test: @test, @testset
 @testset "ITensorsExtensions" begin
@@ -61,64 +62,75 @@ 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")
-
-        # Make a random Hermitian matrix-like 4th order ITensor
-        T = random_itensor(s1', s2', dag(s2), dag(s1))
-        T = apply(T, swapprime(dag(T), 0 => 1))
-        @test T ≈ swapprime(dag(T), 0 => 1) # check Hermitian
-
-        Ul, D, Ur = eigendecomp(T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian = true)
-
-        @test Ul * D * Ur ≈ T
-        ITensors.disable_auto_fermion()
+        with_auto_fermion() do
+            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"
+            )
+
+            # Make a random Hermitian matrix-like 4th order ITensor
+            T = random_itensor(s1', s2', dag(s2), dag(s1))
+            T = apply(T, swapprime(dag(T), 0 => 1))
+            @test T ≈ swapprime(dag(T), 0 => 1) # check Hermitian
+
+            Ul, D, Ur = eigendecomp(T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian = true)
+
+            @test Ul * D * Ur ≈ T
+        end
     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")
-
-        # Make a random Hermitian matrix ITensor
-        M = random_itensor(s1', dag(s1))
-        #M = mapprime(prime(M)*swapprime(dag(M),0=>1),2=>1)
-        M = apply(M, swapprime(dag(M), 0 => 1))
-
-        # Make a random Hermitian matrix-like 4th order ITensor
-        T = random_itensor(s1', s2', dag(s2), dag(s1))
-        T = apply(T, swapprime(dag(T), 0 => 1))
-
-        # Matrix test
-        sqrtM = map_eigvals(sqrt, M, [s1'], [dag(s1)]; ishermitian = true)
-        @test M ≈ apply(sqrtM, sqrtM)
-
-        ## Tensor test
-        sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian = true)
-        @test T ≈ apply(sqrtT, sqrtT)
-
-        # Permute and test again
-        T = permute(T, dag(s2), s2', dag(s1), s1')
-        sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian = true)
-        @test T ≈ apply(sqrtT, sqrtT)
-
-        ## Explicitly passing indices in different, valid orders
-        sqrtT = map_eigvals(sqrt, T, [s2', s1'], [dag(s2), dag(s1)]; ishermitian = true)
-        @test T ≈ apply(sqrtT, sqrtT)
-        sqrtT = map_eigvals(sqrt, T, [dag(s2), dag(s1)], [s2', s1'], ; ishermitian = true)
-        @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()
+        with_auto_fermion() do
+            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"
+            )
+
+            # Make a random Hermitian matrix ITensor
+            M = random_itensor(s1', dag(s1))
+            #M = mapprime(prime(M)*swapprime(dag(M),0=>1),2=>1)
+            M = apply(M, swapprime(dag(M), 0 => 1))
+
+            # Make a random Hermitian matrix-like 4th order ITensor
+            T = random_itensor(s1', s2', dag(s2), dag(s1))
+            T = apply(T, swapprime(dag(T), 0 => 1))
+
+            # Matrix test
+            sqrtM = map_eigvals(sqrt, M, [s1'], [dag(s1)]; ishermitian = true)
+            @test M ≈ apply(sqrtM, sqrtM)
+
+            ## Tensor test
+            sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian = true)
+            @test T ≈ apply(sqrtT, sqrtT)
+
+            # Permute and test again
+            T = permute(T, dag(s2), s2', dag(s1), s1')
+            sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian = true)
+            @test T ≈ apply(sqrtT, sqrtT)
+
+            ## Explicitly passing indices in different, valid orders
+            sqrtT = map_eigvals(sqrt, T, [s2', s1'], [dag(s2), dag(s1)]; ishermitian = true)
+            @test T ≈ apply(sqrtT, sqrtT)
+            sqrtT = map_eigvals(
+                sqrt, T, [dag(s2), dag(s1)], [s2', s1'], ; ishermitian = true
+            )
+            @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)
+        end
     end
 end
 end
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
@@ -1,31 +1,26 @@
 @eval module $(gensym())
 using Graphs: vertices
 using ITensorNetworks: ITensorNetworks, OpSum, siteinds, ttn
-using ITensors: ITensors
+using ITensors.NDTensors: with_auto_fermion
 using NamedGraphs.NamedGraphGenerators: named_grid
 using Test: @test, @testset
 
 @testset "OpSum to TTN converter" begin
     @testset "Multiple onsite terms (regression test for issue #62)" begin
-        auto_fermion_enabled = ITensors.using_auto_fermion()
-        if !auto_fermion_enabled
-            ITensors.enable_auto_fermion()
-        end
-        grid_dims = (2, 1)
-        g = named_grid(grid_dims)
-        s = siteinds("S=1/2", g)
+        with_auto_fermion() do
+            grid_dims = (2, 1)
+            g = named_grid(grid_dims)
+            s = siteinds("S=1/2", g)
 
-        os1 = OpSum()
-        os1 += 1.0, "Sx", (1, 1)
-        os2 = OpSum()
-        os2 += 1.0, "Sy", (1, 1)
-        H1 = ttn(os1, s)
-        H2 = ttn(os2, s)
-        H3 = ttn(os1 + os2, s)
+            os1 = OpSum()
+            os1 += 1.0, "Sx", (1, 1)
+            os2 = OpSum()
+            os2 += 1.0, "Sy", (1, 1)
+            H1 = ttn(os1, s)
+            H2 = ttn(os2, s)
+            H3 = ttn(os1 + os2, s)
 
-        @test H1 + H2 ≈ H3 rtol = 1.0e-6
-        if auto_fermion_enabled
-            ITensors.enable_auto_fermion()
+            @test H1 + H2 ≈ H3 rtol = 1.0e-6
         end
     end
 end
diff --git a/test/test_opsum_to_ttn_mpo_cross_check.jl b/test/test_opsum_to_ttn_mpo_cross_check.jl
@@ -5,9 +5,9 @@ using Graphs: add_edge!, add_vertex!, rem_edge!, vertices
 using ITensorMPS: ITensorMPS
 using ITensorNetworks.ITensorsExtensions: replace_vertices
 using ITensorNetworks: ITensorNetworks, siteinds, ttn
-using ITensors.NDTensors: matrix
-using ITensors: ITensors, @disable_warn_order, ITensor, Index, combinedind, combiner,
-    contract, dag, inds, removeqns
+using ITensors.NDTensors: matrix, with_auto_fermion
+using ITensors: @disable_warn_order, ITensor, Index, combinedind, combiner, contract, dag,
+    inds, removeqns
 using LinearAlgebra: norm
 using NamedGraphs.GraphsExtensions: leaf_vertices, post_order_dfs_vertices
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
@@ -30,7 +30,6 @@ end
 @testset "OpSum to TTN vs ITensorMPS.MPO" begin
     @testset "OpSum to TTN" begin
         # small comb tree
-        auto_fermion_enabled = ITensors.using_auto_fermion()
         tooth_lengths = fill(2, 3)
         c = named_comb_tree(tooth_lengths)
 
@@ -73,9 +72,6 @@ end
             end
             @test Tttno_lr ≈ Tmpo_lr rtol = 1.0e-6
         end
-        if auto_fermion_enabled
-            ITensors.enable_auto_fermion()
-        end
     end
 
     @testset "OpSum to TTN QN" begin
@@ -129,54 +125,50 @@ end
     end
 
     @testset "OpSum to TTN Fermions" begin
-        # small comb tree
-        auto_fermion_enabled = ITensors.using_auto_fermion()
-        if !auto_fermion_enabled
-            ITensors.enable_auto_fermion()
-        end
-        tooth_lengths = fill(2, 3)
-        c = named_comb_tree(tooth_lengths)
-        is = siteinds("Fermion", c; conserve_nf = true)
-
-        # test with next-nearest neighbor tight-binding model
-        t = 1.0
-        tp = 0.4
-        U = 0.0
-        h = 0.5
-        H = ModelHamiltonians.tight_binding(c; t, tp, h)
-
-        # add combination of longer range interactions
-        Hlr = copy(H)
-
-        @testset "Svd approach" for root_vertex in leaf_vertices(is)
-            # get TTN Hamiltonian directly
-            Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10)
-            # get corresponding MPO Hamiltonian
-            sites = [only(is[v]) for v in reverse(post_order_dfs_vertices(c, root_vertex))]
-            vmap = Dictionary(
-                reverse(post_order_dfs_vertices(c, root_vertex)),
-                1:length(sites)
-            )
-            Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites)
-            @disable_warn_order begin
-                Tmpo = prod(Hline)
-                Tttno = contract(Hsvd)
+        with_auto_fermion() do
+            # small comb tree
+            tooth_lengths = fill(2, 3)
+            c = named_comb_tree(tooth_lengths)
+            is = siteinds("Fermion", c; conserve_nf = true)
+
+            # test with next-nearest neighbor tight-binding model
+            t = 1.0
+            tp = 0.4
+            U = 0.0
+            h = 0.5
+            H = ModelHamiltonians.tight_binding(c; t, tp, h)
+
+            # add combination of longer range interactions
+            Hlr = copy(H)
+
+            @testset "Svd approach" for root_vertex in leaf_vertices(is)
+                # get TTN Hamiltonian directly
+                Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10)
+                # get corresponding MPO Hamiltonian
+                sites =
+                    [only(is[v]) for v in reverse(post_order_dfs_vertices(c, root_vertex))]
+                vmap = Dictionary(
+                    reverse(post_order_dfs_vertices(c, root_vertex)),
+                    1:length(sites)
+                )
+                Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites)
+                @disable_warn_order begin
+                    Tmpo = prod(Hline)
+                    Tttno = contract(Hsvd)
+                end
+
+                # verify that the norm isn't 0 and thus the same (which would indicate a problem with the autofermion system
+                @test norm(Tmpo) > 0
+                @test norm(Tttno) > 0
+                @test norm(Tmpo) ≈ norm(Tttno) rtol = 1.0e-6
+
+                # TODO: fix comparison for fermionic tensors
+                @test_broken Tmpo ≈ Tttno
+                # In the meantime: matricize tensors and convert to dense Matrix to compare element by element
+                dTmm = to_matrix(Tmpo)
+                dTtm = to_matrix(Tttno)
+                @test any(>(1.0e-14), dTmm - dTtm)
             end
-
-            # verify that the norm isn't 0 and thus the same (which would indicate a problem with the autofermion system
-            @test norm(Tmpo) > 0
-            @test norm(Tttno) > 0
-            @test norm(Tmpo) ≈ norm(Tttno) rtol = 1.0e-6
-
-            # TODO: fix comparison for fermionic tensors
-            @test_broken Tmpo ≈ Tttno
-            # In the meantime: matricize tensors and convert to dense Matrix to compare element by element
-            dTmm = to_matrix(Tmpo)
-            dTtm = to_matrix(Tttno)
-            @test any(>(1.0e-14), dTmm - dTtm)
-        end
-        if !auto_fermion_enabled
-            ITensors.disable_auto_fermion()
         end
     end
 
diff --git a/test/test_ttn_position.jl b/test/test_ttn_position.jl
