From d6914af540bbe76d2e430921bc3210bf56a6e836 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 15 Jul 2026 08:34:47 -0600 Subject: [PATCH 01/16] Implement coefficient thresholding in Hamiltonian generation - Introduced `hamiltonian.coefficient_threshold` to control the absolute cutoff for spin-orbital and final Pauli coefficients. - Updated relevant functions to utilize the new threshold setting, allowing for disabling magnitude-based pruning. - Documented changes in README and added a detailed hotfix explanation in THRESHOLDING_HOTFIX.md. --- hamiltonian_generator/THRESHOLDING_HOTFIX.md | 72 ++++++++ hamiltonian_generator/thresholding.py | 181 +++++++++++++++++++ 2 files changed, 253 insertions(+) create mode 100644 hamiltonian_generator/THRESHOLDING_HOTFIX.md create mode 100644 hamiltonian_generator/thresholding.py diff --git a/hamiltonian_generator/THRESHOLDING_HOTFIX.md b/hamiltonian_generator/THRESHOLDING_HOTFIX.md new file mode 100644 index 00000000..318bd93b --- /dev/null +++ b/hamiltonian_generator/THRESHOLDING_HOTFIX.md @@ -0,0 +1,72 @@ +# Issue #31 thresholding hotfix + +This is a QHAT-only hotfix. It does not edit or monkey-patch OpenFermion. + +## What changed + +Add one setting to a Hamiltonian-generator configuration: + +```python +# Existing/default behavior +hamiltonian.coefficient_threshold = 1e-8 + +# Disable magnitude thresholding +hamiltonian.coefficient_threshold = 0.0 +``` + +QHAT now owns two small operations that previously had hidden OpenFermion cutoffs: + +1. Spatial integrals are expanded into spin-orbital tensors in + `hamiltonian_generator/thresholding.py`, using the configured threshold. +2. JW/BK mapping is done one unit-coefficient fermion term at a time. QHAT restores the actual + coefficients, sums matching Pauli strings, and applies the configured threshold only once at + the end. + +The second step matters because OpenFermion also removes small coefficients while incrementally +building a `QubitOperator`. Fixing only the tensor conversion does not disable thresholding in the +final Pauli file. + +## Cache handling + +The active-space pickle already contains thresholded tensors. QHAT now records the threshold on +that object and checks it when loading the cache. If the requested value differs, QHAT reloads the +Hartree-Fock pickle and rebuilds the active-space tensors. Old pickles are treated as `1e-8`, which +was the only behavior before this hotfix. + +The output filenames are unchanged. Running the same `general.file_stub` with a new threshold +overwrites the active-space and final files, so use another stub if both versions must be kept. + +This hotfix is intentionally limited to `hamiltonian_generator/hamgen.py`. Its generated `.dat` +Pauli file follows the configured threshold. The separate analysis path that remaps a +`.tensors.npz` file is unchanged and still uses OpenFermion's normal mapping behavior. + +## Li2 result + +Using the bundled five-spatial-orbital Li2 fixture: + +| Threshold | One-body tensor entries | Two-body tensor entries | JW/BK Pauli terms | Ground energy (Ha) | +| ---: | ---: | ---: | ---: | ---: | +| `1e-8` | 14 | 436 | 156 | about `-14.6426835699246` | +| `0.0` | 50 | 2500 | 876 | about `-14.6426835699246` | + +This reproduces the Issue #31 observation: `1e-8` removes 36 one-body and 2,064 two-body tensor +entries. The one-norm of the difference between the final Hamiltonians is only about `1.30e-14` +Ha, so their ground energies are numerically equivalent for this fixture. + +## Important zero-threshold caveat + +`0.0` keeps floating-point symmetry and cancellation residue as well as intentionally small terms. +The extra Li2 terms are mostly machine-scale noise, so zero does **not** guarantee identical term +counts at every bond length. A small positive threshold is normally more stable; zero is useful +for auditing and sensitivity checks. + +For `m` dense active spatial orbitals, the threshold-free spin tensors have `2*m^2` one-body and +`4*m^4` two-body entries. For this Li2 case (`m = 5`), that gives 50 and 2,500. + +## Files in the hotfix + +- `hamiltonian_generator/hamgen.py`: uses the configurable QHAT path and validates cached tensors. +- `hamiltonian_generator/hamgen_types.py`: exposes and validates the configuration setting. +- `hamiltonian_generator/thresholding.py`: contains the local tensor conversion and JW/BK mapper. +- `hamiltonian_generator/config.py`: shows the setting in the example configuration. +- `hamiltonian_generator/README.md`: documents normal usage. diff --git a/hamiltonian_generator/thresholding.py b/hamiltonian_generator/thresholding.py new file mode 100644 index 00000000..ef98bd84 --- /dev/null +++ b/hamiltonian_generator/thresholding.py @@ -0,0 +1,181 @@ +"""Small QHAT-owned hotfix for OpenFermion's fixed coefficient cutoff.""" + +from collections import defaultdict +import math +import numbers + +import numpy as np +from openfermion import ( + FermionOperator, + InteractionOperator, + QubitOperator, + bravyi_kitaev, + jordan_wigner, +) + + +DEFAULT_COEFFICIENT_THRESHOLD = 1.0e-8 + + +def validate_coefficient_threshold(value) -> float: + if isinstance(value, bool) or not isinstance(value, numbers.Real): + raise TypeError("coefficient threshold must be a real number") + threshold = float(value) + if not math.isfinite(threshold) or threshold < 0.0: + raise ValueError("coefficient threshold must be finite and non-negative") + return threshold + + +def _validate_spatial_integrals(one_body_integrals, two_body_integrals): + one_body = np.asarray(one_body_integrals) + two_body = np.asarray(two_body_integrals) + + if one_body.ndim != 2 or one_body.shape[0] != one_body.shape[1]: + raise ValueError("one-body spatial integrals must be a square rank-2 array") + expected_two_body_shape = (one_body.shape[0],) * 4 + if two_body.shape != expected_two_body_shape: + raise ValueError( + "two-body spatial integrals must have shape " + f"{expected_two_body_shape}, got {two_body.shape}" + ) + return one_body, two_body + + +def spinorb_from_spatial( + one_body_integrals, + two_body_integrals, + coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, +): + """Expand spatial integrals into spin-orbital tensors. + + This is the narrow piece of ``MolecularData.get_molecular_hamiltonian`` + QHAT needs to own in order to expose the otherwise fixed OpenFermion + cutoff. It preserves OpenFermion's strict ``abs(value) < threshold`` + boundary and supports ``0.0`` to disable magnitude truncation. + """ + + threshold = validate_coefficient_threshold(coefficient_threshold) + one_spatial, two_spatial = _validate_spatial_integrals( + one_body_integrals, two_body_integrals + ) + n_qubits = 2 * one_spatial.shape[0] + + one_spin = np.zeros((n_qubits, n_qubits), dtype=one_spatial.dtype) + two_spin = np.zeros((n_qubits,) * 4, dtype=two_spatial.dtype) + + one_spin[0::2, 0::2] = one_spatial + one_spin[1::2, 1::2] = one_spatial + + two_spin[0::2, 1::2, 1::2, 0::2] = two_spatial + two_spin[1::2, 0::2, 0::2, 1::2] = two_spatial + two_spin[0::2, 0::2, 0::2, 0::2] = two_spatial + two_spin[1::2, 1::2, 1::2, 1::2] = two_spatial + + if threshold > 0.0: + one_spin[np.abs(one_spin) < threshold] = 0.0 + two_spin[np.abs(two_spin) < threshold] = 0.0 + + return one_spin, two_spin + + +def molecular_hamiltonian_from_spatial( + constant, + one_body_integrals, + two_body_integrals, + coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, +): + """Construct an ``InteractionOperator`` and record cutoff provenance.""" + + threshold = validate_coefficient_threshold(coefficient_threshold) + one_spatial, two_spatial = _validate_spatial_integrals( + one_body_integrals, two_body_integrals + ) + one_before = 2 * int(np.count_nonzero(one_spatial)) + two_before = 4 * int(np.count_nonzero(two_spatial)) + one_spin, two_spin = spinorb_from_spatial(one_spatial, two_spatial, threshold) + + operator = InteractionOperator(constant, one_spin, 0.5 * two_spin) + operator.coefficient_threshold = threshold + operator.threshold_stats = { + "one_body_before": one_before, + "one_body_after": int(np.count_nonzero(one_spin)), + "one_body_removed": one_before - int(np.count_nonzero(one_spin)), + "two_body_before": two_before, + "two_body_after": int(np.count_nonzero(two_spin)), + "two_body_removed": two_before - int(np.count_nonzero(two_spin)), + } + return operator + + +def get_molecular_hamiltonian( + molecule, + occupied_indices=None, + active_indices=None, + coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, +): + """QHAT-local counterpart of ``MolecularData.get_molecular_hamiltonian``.""" + + if occupied_indices is None and active_indices is None: + one_body_integrals, two_body_integrals = molecule.get_integrals() + constant = molecule.nuclear_repulsion + else: + core_adjustment, one_body_integrals, two_body_integrals = ( + molecule.get_active_space_integrals(occupied_indices, active_indices) + ) + constant = molecule.nuclear_repulsion + core_adjustment + + return molecular_hamiltonian_from_spatial( + constant, + one_body_integrals, + two_body_integrals, + coefficient_threshold, + ) + + +def map_interaction_operator( + operator, + mapping_name, + coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, +): + """Map without changing OpenFermion's global ``EQ_TOLERANCE``. + + Each unit-coefficient monomial safely survives OpenFermion's internal + cutoff. QHAT restores the real coefficient, sums matching Pauli terms, + and applies the requested cutoff once at the end. + """ + + threshold = validate_coefficient_threshold(coefficient_threshold) + if mapping_name not in ("jordan-wigner", "bravyi-kitaev"): + raise ValueError(f"invalid fermion-to-qubit mapping {mapping_name!r}") + + contributions = defaultdict(list) + for fermion_term in operator: + coefficient = operator[fermion_term] + unit_term = FermionOperator(fermion_term, 1.0) + if mapping_name == "jordan-wigner": + mapped_unit = jordan_wigner(unit_term) + else: + mapped_unit = bravyi_kitaev(unit_term, n_qubits=operator.n_qubits) + for pauli_term, unit_coefficient in mapped_unit.terms.items(): + contributions[pauli_term].append(coefficient * unit_coefficient) + + terms = {} + for pauli_term, values in contributions.items(): + coefficient = complex( + math.fsum(complex(value).real for value in values), + math.fsum(complex(value).imag for value in values), + ) + # Molecular Hamiltonians are Hermitian. Permit only machine-scale + # imaginary residue from floating-point tensor symmetries. + scale = math.fsum(abs(value) for value in values) + imag_tolerance = 128.0 * np.finfo(float).eps * max(1.0, scale) + if abs(coefficient.imag) > imag_tolerance: + raise ValueError(f"non-Hermitian Pauli coefficient {coefficient!r}") + coefficient = coefficient.real + if coefficient != 0.0 and abs(coefficient) >= threshold: + terms[pauli_term] = coefficient + + result = QubitOperator() + # Do not use += here; it would reintroduce OpenFermion's fixed cutoff. + result.terms = terms + return result From 1429097d9f040e301149481dda68e081bdad2f36 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 15 Jul 2026 08:34:55 -0600 Subject: [PATCH 02/16] Add coefficient thresholding functionality to Hamiltonian generation - Introduced `coefficient_threshold` in configuration with a default of 1e-8. - Updated README to document the coefficient cutoff behavior. - Modified `hamgen.py` to apply the coefficient threshold during Hamiltonian construction. - Enhanced metadata logging to include threshold statistics. --- hamiltonian_generator/README.md | 19 +++- hamiltonian_generator/config.py | 5 +- hamiltonian_generator/hamgen.py | 130 +++++++++++++++++++------- hamiltonian_generator/hamgen_types.py | 13 ++- 4 files changed, 130 insertions(+), 37 deletions(-) diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index b56645e9..a53acb4e 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -74,6 +74,21 @@ Currently allowed values are - for Bravyi-Kitaev: "Bravyi-Kitaev", "Bravyi Kitaev", "BK", "bravyi-kitaev", "bravyi\_kitaev", "bravyi kitaev", or "bk" +The absolute coefficient cutoff is set with **`hamiltonian.coefficient_threshold`**. Its default +is `1e-8`, preserving the historical QHAT/OpenFermion scale. QHAT applies this cutoff when it +builds the active-space spin-orbital tensors and again after all contributions to each final Pauli +coefficient have been summed. Set it to `0.0` to disable magnitude-based pruning: + +```python +hamiltonian.coefficient_threshold = 0.0 +``` + +The value must be finite and non-negative. The comparison is strict, so a coefficient exactly +equal to the threshold is retained. Literal zero also retains floating-point symmetry noise, so +it does not guarantee identical term counts for molecular geometries that differ only in bond +length. See [THRESHOLDING_HOTFIX.md](THRESHOLDING_HOTFIX.md) for the Issue #31 diagnosis, +Li2 results, cache behavior, and verification instructions. + ## Generated Files This script will generate various files for intermediate and final stages of the process. These @@ -91,7 +106,9 @@ locking some high-energy orbitals as being vacant. The results of this step wil active occupied orbitals and the number of active vacant orbitals. The one-body and two-body tensors from this stage of the calculation will also be saved in a file called "[filestub]\_[astag].tensors.npz", which can be loaded into the resource estimation software using -the `load_numpy` function. +the `load_numpy` function. Both files record the effective coefficient threshold. If a cached +pickle was generated with a different threshold, QHAT rebuilds this stage from the Hartree-Fock +pickle instead of silently reusing incompatible tensors. The third stage is to transform the fermionic creation and annihilation operators to qubit operators, yielding some metadata and a set of Pauli strings. This information will be written to diff --git a/hamiltonian_generator/config.py b/hamiltonian_generator/config.py index ce166da0..679d7be2 100644 --- a/hamiltonian_generator/config.py +++ b/hamiltonian_generator/config.py @@ -14,6 +14,9 @@ hamiltonian.basis = "sto-3g" # Select the atomic basis functions +# Absolute cutoff applied to spin-orbital and final Pauli coefficients. The default is 1e-8; +# set this to 0.0 to retain every nonzero floating-point coefficient. +hamiltonian.coefficient_threshold = 1e-8 + hamiltonian.num_active_occupied = 4 # Specify the active space hamiltonian.num_active_vacant = 6 - diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 6c82d943..18defe7d 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -8,7 +8,7 @@ import mendeleev import numpy as np -from openfermion import InteractionOperator, MolecularData, bravyi_kitaev, jordan_wigner +from openfermion import MolecularData from openfermion.transforms.opconversions.binary_codes import _encoder_bk from openfermionpyscf import PyscfMolecularData from openfermionpyscf._run_pyscf import compute_integrals, compute_scf, prepare_pyscf_molecule @@ -20,6 +20,11 @@ HamiltonianConfiguration, State, ) +from qhat.hamiltonian_generator.thresholding import ( + DEFAULT_COEFFICIENT_THRESHOLD, + get_molecular_hamiltonian, + map_interaction_operator, +) logger = logging.getLogger(__name__) @@ -245,10 +250,18 @@ def apply_active_space(state, ham1_HartreeFock): # Build an InteractionOperator (from openfermion) with one- and two-body integrals computed state.log_verbose("Build instance of InteractionOperator.") - ham2_ActiveSpace = ham1_HartreeFock.get_molecular_hamiltonian( + ham2_ActiveSpace = get_molecular_hamiltonian( + ham1_HartreeFock, occupied_indices=range(idx_act_occ), active_indices=range(idx_act_occ, idx_frz_vac), + coefficient_threshold=state.config_hamiltonian.coefficient_threshold, ) + stats = ham2_ActiveSpace.threshold_stats + state.log(" ".join([ + f"Coefficient threshold {ham2_ActiveSpace.coefficient_threshold:.6g} retained", + f"{stats['one_body_after']}/{stats['one_body_before']} one-body and", + f"{stats['two_body_after']}/{stats['two_body_before']} two-body tensor entries.", + ])) tstop = time.time() # Pass along metadata that needs to be preserved ham2_ActiveSpace.hf_energy = ham1_HartreeFock.hf_energy @@ -266,13 +279,11 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): tstart = time.time() ham3_Fermion2Qubit = None mapping_name = state.config_hamiltonian.fermion_to_qubit_name() - if mapping_name == "jordan-wigner": - ham3_Fermion2Qubit = jordan_wigner(ham2_ActiveSpace) - elif mapping_name == "bravyi-kitaev": - ham3_Fermion2Qubit = bravyi_kitaev(ham2_ActiveSpace) - else: - mapping = state.config_hamiltonian.f2q_mapping() - raise NotImplementedError(f"invalid fermion-to-qubit mapping \"{mapping}\"") + ham3_Fermion2Qubit = map_interaction_operator( + ham2_ActiveSpace, + mapping_name, + state.config_hamiltonian.coefficient_threshold, + ) tstop = time.time() # Propagate previously-computed metadata that we need to preserve across all run modes ham3_Fermion2Qubit.hf_energy = ham2_ActiveSpace.hf_energy @@ -280,6 +291,8 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): ham3_Fermion2Qubit.basis = ham2_ActiveSpace.basis ham3_Fermion2Qubit.separation = ham2_ActiveSpace.separation ham3_Fermion2Qubit.f2q_mapping = mapping_name + ham3_Fermion2Qubit.coefficient_threshold = ham2_ActiveSpace.coefficient_threshold + ham3_Fermion2Qubit.threshold_stats = ham2_ActiveSpace.threshold_stats ham3_Fermion2Qubit.hf_time = ham2_ActiveSpace.hf_time ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time ham3_Fermion2Qubit.f2q_time = tstop - tstart @@ -366,32 +379,67 @@ def get_ham2(state): # TODO: Python errors get a little weird if you have an exception inside an exception. So # instead the exception clause should _only_ flag that we're going to recompute versus # load, and the actual recomputing should be outside of the except clause. - except FileNotFoundError as err: + except FileNotFoundError: + ham2_ActiveSpace = None + + requested_threshold = state.config_hamiltonian.coefficient_threshold + if ham2_ActiveSpace is not None: + # Pickles created before threshold configurability necessarily used + # OpenFermion's historical 1e-8 default. + cached_threshold = getattr( + ham2_ActiveSpace, + "coefficient_threshold", + DEFAULT_COEFFICIENT_THRESHOLD, + ) + if cached_threshold == requested_threshold: + ham2_ActiveSpace.coefficient_threshold = cached_threshold + if not hasattr(ham2_ActiveSpace, "threshold_stats"): + one_after = int(np.count_nonzero(ham2_ActiveSpace.one_body_tensor)) + two_after = int(np.count_nonzero(ham2_ActiveSpace.two_body_tensor)) + ham2_ActiveSpace.threshold_stats = { + "one_body_before": None, + "one_body_after": one_after, + "one_body_removed": None, + "two_body_before": None, + "two_body_after": two_after, + "two_body_removed": None, + } + state.log(' '.join([ + f"Loaded \"{ham2_filename}\" with coefficient threshold", + f"{cached_threshold:.6g}.", + "Continuing from after the active space is applied."])) + return ham2_ActiveSpace + state.log(' '.join([ - f"Could not load \"{ham2_filename}\".", - f"Trying to load \"{state.filename_ham1()}\"."])) - # Get ham1_HartreeFock (by loading or by recomputing, depending on data availability) - ham1_HartreeFock = get_ham1(state) - # Recompute ham2_ActiveSpace from ham1_HartreeFock - state.log("Apply active space.") - ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) - state.log(f"Pickle to \"{ham2_filename}\" file.") - # Save ham2_ActiveSpace for later re-use - with open(ham2_filename, 'wb') as ham2_file: - pickle.dump(ham2_ActiveSpace, ham2_file) - # Save the one-body and two-body tensors using numpy's `save` function. These files can be - # loaded using numpy's `load` function. - t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" - np.savez_compressed(t_filename, - constant=ham2_ActiveSpace.constant, - one_body=ham2_ActiveSpace.one_body_tensor, - two_body=ham2_ActiveSpace.two_body_tensor) - return ham2_ActiveSpace + f"Cached \"{ham2_filename}\" used coefficient threshold", + f"{cached_threshold:.6g}, but {requested_threshold:.6g} was requested.", + "Recomputing the active-space tensors."])) else: - state.log(' '.join([ - f"Loaded \"{ham2_filename}\".", - "Continuing from after the active space is applied."])) - return ham2_ActiveSpace + state.log(f"Could not load \"{ham2_filename}\".") + + state.log(f"Trying to load \"{state.filename_ham1()}\".") + ham1_HartreeFock = get_ham1(state) + state.log("Apply active space.") + ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) + state.log(f"Pickle to \"{ham2_filename}\" file.") + with open(ham2_filename, 'wb') as ham2_file: + pickle.dump(ham2_ActiveSpace, ham2_file) + + # Store threshold provenance alongside the tensors. + stats = ham2_ActiveSpace.threshold_stats + t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" + np.savez_compressed( + t_filename, + constant=ham2_ActiveSpace.constant, + one_body=ham2_ActiveSpace.one_body_tensor, + two_body=ham2_ActiveSpace.two_body_tensor, + coefficient_threshold=ham2_ActiveSpace.coefficient_threshold, + one_body_entries_before_threshold=stats["one_body_before"], + one_body_entries_after_threshold=stats["one_body_after"], + two_body_entries_before_threshold=stats["two_body_before"], + two_body_entries_after_threshold=stats["two_body_after"], + ) + return ham2_ActiveSpace # ------------------------------------------------------------------------------------------------- @@ -449,13 +497,27 @@ def compute_metadata(state, ham3_Fermion2Qubit): state.metadata["interatomic separation (angstroms)"] = ham3_Fermion2Qubit.separation # fermion-to-qubit transformation state.metadata["fermion-to-qubit operator mapping"] = ham3_Fermion2Qubit.f2q_mapping + state.metadata["coefficient threshold (Hartrees)"] = ( + ham3_Fermion2Qubit.coefficient_threshold + ) + threshold_stats = ham3_Fermion2Qubit.threshold_stats + for tensor_name in ("one_body", "two_body"): + before = threshold_stats[f"{tensor_name}_before"] + after = threshold_stats[f"{tensor_name}_after"] + if before is not None: + state.metadata[ + f"{tensor_name.replace('_', '-')} tensor entries before thresholding" + ] = before + state.metadata[ + f"{tensor_name.replace('_', '-')} tensor entries after thresholding" + ] = after # number of terms in sum of Pauli strings state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms) # one-norm of sum of Pauli strings one_norm = sum(abs(coefficient) for coefficient in ham3_Fermion2Qubit.terms.values()) state.metadata["one-norm of sum of Pauli strings (Hartrees)"] = one_norm # eigenvalue bounds (see Trotter Workflow document for where this comes from) - E0 = ham3_Fermion2Qubit.terms[tuple()] + E0 = ham3_Fermion2Qubit.terms.get(tuple(), 0.0) dE = one_norm - abs(E0) state.metadata["eigenvalue lower bound from shifted one-norm (Hartrees)"] = E0 - dE state.metadata["eigenvalue upper bound from shifted one-norm (Hartrees)"] = E0 + dE diff --git a/hamiltonian_generator/hamgen_types.py b/hamiltonian_generator/hamgen_types.py index 3638463e..c6d6269b 100644 --- a/hamiltonian_generator/hamgen_types.py +++ b/hamiltonian_generator/hamgen_types.py @@ -6,6 +6,11 @@ import subprocess import sys +from qhat.hamiltonian_generator.thresholding import ( + DEFAULT_COEFFICIENT_THRESHOLD, + validate_coefficient_threshold, +) + logger = logging.getLogger(__name__) # ------------------------------------------------------------------------------------------------- @@ -34,9 +39,16 @@ def __init__(self): self._geometry = () self.basis = "sto-3g" self.f2q_mapping = "Jordan-Wigner" + self._coefficient_threshold = DEFAULT_COEFFICIENT_THRESHOLD # The number of active occupied and vacant spin (single-occupancy) orbitals in the molecule self.num_active_occupied = None self.num_active_vacant = None + @property + def coefficient_threshold(self): + return self._coefficient_threshold + @coefficient_threshold.setter + def coefficient_threshold(self, value): + self._coefficient_threshold = validate_coefficient_threshold(value) def add_atom(self, element, x, y, z): # TODO: Clarify the units involved. Right now we just pass numbers along, so the user has # to know the right units for the various choices they make. It would be better to @@ -148,4 +160,3 @@ def filename_ham3(self): ext=self.config_general.ham3_ext()) # ------------------------------------------------------------------------------------------------- - From 8d020d7ab077f116753efef10078dd9a451a7d9d Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Thu, 16 Jul 2026 07:44:03 -0600 Subject: [PATCH 03/16] Refactor Hamiltonian coefficient thresholding documentation and implementation --- hamiltonian_generator/README.md | 20 +-- hamiltonian_generator/THRESHOLDING_HOTFIX.md | 80 +++------- hamiltonian_generator/config.py | 3 +- hamiltonian_generator/hamgen.py | 148 ++++++------------ hamiltonian_generator/thresholding.py | 154 ++----------------- 5 files changed, 97 insertions(+), 308 deletions(-) diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index a53acb4e..ec6c8317 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -74,20 +74,18 @@ Currently allowed values are - for Bravyi-Kitaev: "Bravyi-Kitaev", "Bravyi Kitaev", "BK", "bravyi-kitaev", "bravyi\_kitaev", "bravyi kitaev", or "bk" -The absolute coefficient cutoff is set with **`hamiltonian.coefficient_threshold`**. Its default -is `1e-8`, preserving the historical QHAT/OpenFermion scale. QHAT applies this cutoff when it -builds the active-space spin-orbital tensors and again after all contributions to each final Pauli -coefficient have been summed. Set it to `0.0` to disable magnitude-based pruning: +**`hamiltonian.coefficient_threshold`** controls the absolute cutoff used when QHAT expands +active-space spatial integrals into spin-orbital tensors. Its default is `1e-8`; set it to `0.0` +to disable filtering at this stage: ```python hamiltonian.coefficient_threshold = 0.0 ``` -The value must be finite and non-negative. The comparison is strict, so a coefficient exactly -equal to the threshold is retained. Literal zero also retains floating-point symmetry noise, so -it does not guarantee identical term counts for molecular geometries that differ only in bond -length. See [THRESHOLDING_HOTFIX.md](THRESHOLDING_HOTFIX.md) for the Issue #31 diagnosis, -Li2 results, cache behavior, and verification instructions. +The value must be finite and non-negative. This does not change OpenFermion's standard JW/BK +mapping, which can still discard small coefficients. Changing the value also does not invalidate +an existing active-space pickle; remove that pickle or use a different `general.file_stub` when +comparing thresholds. See [THRESHOLDING_HOTFIX.md](THRESHOLDING_HOTFIX.md) for details. ## Generated Files @@ -106,9 +104,7 @@ locking some high-energy orbitals as being vacant. The results of this step wil active occupied orbitals and the number of active vacant orbitals. The one-body and two-body tensors from this stage of the calculation will also be saved in a file called "[filestub]\_[astag].tensors.npz", which can be loaded into the resource estimation software using -the `load_numpy` function. Both files record the effective coefficient threshold. If a cached -pickle was generated with a different threshold, QHAT rebuilds this stage from the Hartree-Fock -pickle instead of silently reusing incompatible tensors. +the `load_numpy` function. The third stage is to transform the fermionic creation and annihilation operators to qubit operators, yielding some metadata and a set of Pauli strings. This information will be written to diff --git a/hamiltonian_generator/THRESHOLDING_HOTFIX.md b/hamiltonian_generator/THRESHOLDING_HOTFIX.md index 318bd93b..2221eb96 100644 --- a/hamiltonian_generator/THRESHOLDING_HOTFIX.md +++ b/hamiltonian_generator/THRESHOLDING_HOTFIX.md @@ -1,72 +1,42 @@ -# Issue #31 thresholding hotfix +# Issue #31 tensor-stage threshold hotfix -This is a QHAT-only hotfix. It does not edit or monkey-patch OpenFermion. +This is a narrow QHAT-only change. OpenFermion is not edited or monkey-patched. ## What changed -Add one setting to a Hamiltonian-generator configuration: +The Hamiltonian configuration now accepts a finite, non-negative cutoff: ```python -# Existing/default behavior -hamiltonian.coefficient_threshold = 1e-8 - -# Disable magnitude thresholding -hamiltonian.coefficient_threshold = 0.0 +hamiltonian.coefficient_threshold = 1e-8 # default +hamiltonian.coefficient_threshold = 0.0 # no tensor-stage filtering ``` -QHAT now owns two small operations that previously had hidden OpenFermion cutoffs: - -1. Spatial integrals are expanded into spin-orbital tensors in - `hamiltonian_generator/thresholding.py`, using the configured threshold. -2. JW/BK mapping is done one unit-coefficient fermion term at a time. QHAT restores the actual - coefficients, sums matching Pauli strings, and applies the configured threshold only once at - the end. +QHAT uses a local `spinorb_from_spatial` function when it expands active-space spatial integrals +into spin-orbital tensors. The rest of the generator stays on its original path: -The second step matters because OpenFermion also removes small coefficients while incrementally -building a `QubitOperator`. Fixing only the tensor conversion does not disable thresholding in the -final Pauli file. +- OpenFermion's standard Jordan-Wigner and Bravyi-Kitaev functions are unchanged. +- The active-space pickle behavior is unchanged. +- The `.tensors.npz` format remains `constant`, `one_body`, and `two_body`. +- Metadata records only the configured spin-orbital threshold. -## Cache handling +## Limitations -The active-space pickle already contains thresholded tensors. QHAT now records the threshold on -that object and checks it when loading the cache. If the requested value differs, QHAT reloads the -Hartree-Fock pickle and rebuilds the active-space tensors. Old pickles are treated as `1e-8`, which -was the only behavior before this hotfix. +OpenFermion's mapping functions can still remove small coefficients. Setting this option to +`0.0` therefore disables filtering only during tensor construction, not in the final Pauli sum. +Pauli-level thresholding is a separate follow-up. -The output filenames are unchanged. Running the same `general.file_stub` with a new threshold -overwrites the active-space and final files, so use another stub if both versions must be kept. +An existing `[stub]_[astag].pickle` is reused without checking its threshold. When changing the +setting, remove that active-space pickle or use a new `general.file_stub`. Removing only the +`.tensors.npz` file is not enough. -This hotfix is intentionally limited to `hamiltonian_generator/hamgen.py`. Its generated `.dat` -Pauli file follows the configured threshold. The separate analysis path that remaps a -`.tensors.npz` file is unchanged and still uses OpenFermion's normal mapping behavior. +## Li2 check -## Li2 result +For the bundled five-spatial-orbital Li2 case: -Using the bundled five-spatial-orbital Li2 fixture: - -| Threshold | One-body tensor entries | Two-body tensor entries | JW/BK Pauli terms | Ground energy (Ha) | +| Tensor threshold | One-body entries | Two-body entries | Stock JW terms | Stock BK terms | | ---: | ---: | ---: | ---: | ---: | -| `1e-8` | 14 | 436 | 156 | about `-14.6426835699246` | -| `0.0` | 50 | 2500 | 876 | about `-14.6426835699246` | - -This reproduces the Issue #31 observation: `1e-8` removes 36 one-body and 2,064 two-body tensor -entries. The one-norm of the difference between the final Hamiltonians is only about `1.30e-14` -Ha, so their ground energies are numerically equivalent for this fixture. - -## Important zero-threshold caveat - -`0.0` keeps floating-point symmetry and cancellation residue as well as intentionally small terms. -The extra Li2 terms are mostly machine-scale noise, so zero does **not** guarantee identical term -counts at every bond length. A small positive threshold is normally more stable; zero is useful -for auditing and sensitivity checks. - -For `m` dense active spatial orbitals, the threshold-free spin tensors have `2*m^2` one-body and -`4*m^4` two-body entries. For this Li2 case (`m = 5`), that gives 50 and 2,500. - -## Files in the hotfix +| `1e-8` | 14 | 436 | 156 | 156 | +| `0.0` | 50 | 2500 | 156 | 156 | -- `hamiltonian_generator/hamgen.py`: uses the configurable QHAT path and validates cached tensors. -- `hamiltonian_generator/hamgen_types.py`: exposes and validates the configuration setting. -- `hamiltonian_generator/thresholding.py`: contains the local tensor conversion and JW/BK mapper. -- `hamiltonian_generator/config.py`: shows the setting in the example configuration. -- `hamiltonian_generator/README.md`: documents normal usage. +The tensor counts confirm that the setting works at the intended stage. The unchanged mapped +term counts confirm that final Pauli thresholding remains outside this hotfix. diff --git a/hamiltonian_generator/config.py b/hamiltonian_generator/config.py index 679d7be2..aeecb57c 100644 --- a/hamiltonian_generator/config.py +++ b/hamiltonian_generator/config.py @@ -14,8 +14,7 @@ hamiltonian.basis = "sto-3g" # Select the atomic basis functions -# Absolute cutoff applied to spin-orbital and final Pauli coefficients. The default is 1e-8; -# set this to 0.0 to retain every nonzero floating-point coefficient. +# Absolute cutoff applied while building spin-orbital tensors; the default is 1e-8. hamiltonian.coefficient_threshold = 1e-8 hamiltonian.num_active_occupied = 4 # Specify the active space diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 18defe7d..9f16e9d7 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -8,7 +8,7 @@ import mendeleev import numpy as np -from openfermion import MolecularData +from openfermion import InteractionOperator, MolecularData, bravyi_kitaev, jordan_wigner from openfermion.transforms.opconversions.binary_codes import _encoder_bk from openfermionpyscf import PyscfMolecularData from openfermionpyscf._run_pyscf import compute_integrals, compute_scf, prepare_pyscf_molecule @@ -20,11 +20,7 @@ HamiltonianConfiguration, State, ) -from qhat.hamiltonian_generator.thresholding import ( - DEFAULT_COEFFICIENT_THRESHOLD, - get_molecular_hamiltonian, - map_interaction_operator, -) +from qhat.hamiltonian_generator.thresholding import spinorb_from_spatial logger = logging.getLogger(__name__) @@ -250,18 +246,22 @@ def apply_active_space(state, ham1_HartreeFock): # Build an InteractionOperator (from openfermion) with one- and two-body integrals computed state.log_verbose("Build instance of InteractionOperator.") - ham2_ActiveSpace = get_molecular_hamiltonian( - ham1_HartreeFock, - occupied_indices=range(idx_act_occ), - active_indices=range(idx_act_occ, idx_frz_vac), - coefficient_threshold=state.config_hamiltonian.coefficient_threshold, + core_adjustment, one_body_integrals, two_body_integrals = ( + ham1_HartreeFock.get_active_space_integrals( + occupied_indices=range(idx_act_occ), + active_indices=range(idx_act_occ, idx_frz_vac), + ) + ) + one_body_coefficients, two_body_coefficients = spinorb_from_spatial( + one_body_integrals, + two_body_integrals, + state.config_hamiltonian.coefficient_threshold, + ) + ham2_ActiveSpace = InteractionOperator( + ham1_HartreeFock.nuclear_repulsion + core_adjustment, + one_body_coefficients, + 0.5 * two_body_coefficients, ) - stats = ham2_ActiveSpace.threshold_stats - state.log(" ".join([ - f"Coefficient threshold {ham2_ActiveSpace.coefficient_threshold:.6g} retained", - f"{stats['one_body_after']}/{stats['one_body_before']} one-body and", - f"{stats['two_body_after']}/{stats['two_body_before']} two-body tensor entries.", - ])) tstop = time.time() # Pass along metadata that needs to be preserved ham2_ActiveSpace.hf_energy = ham1_HartreeFock.hf_energy @@ -279,11 +279,13 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): tstart = time.time() ham3_Fermion2Qubit = None mapping_name = state.config_hamiltonian.fermion_to_qubit_name() - ham3_Fermion2Qubit = map_interaction_operator( - ham2_ActiveSpace, - mapping_name, - state.config_hamiltonian.coefficient_threshold, - ) + if mapping_name == "jordan-wigner": + ham3_Fermion2Qubit = jordan_wigner(ham2_ActiveSpace) + elif mapping_name == "bravyi-kitaev": + ham3_Fermion2Qubit = bravyi_kitaev(ham2_ActiveSpace) + else: + mapping = state.config_hamiltonian.f2q_mapping() + raise NotImplementedError(f"invalid fermion-to-qubit mapping \"{mapping}\"") tstop = time.time() # Propagate previously-computed metadata that we need to preserve across all run modes ham3_Fermion2Qubit.hf_energy = ham2_ActiveSpace.hf_energy @@ -291,8 +293,6 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): ham3_Fermion2Qubit.basis = ham2_ActiveSpace.basis ham3_Fermion2Qubit.separation = ham2_ActiveSpace.separation ham3_Fermion2Qubit.f2q_mapping = mapping_name - ham3_Fermion2Qubit.coefficient_threshold = ham2_ActiveSpace.coefficient_threshold - ham3_Fermion2Qubit.threshold_stats = ham2_ActiveSpace.threshold_stats ham3_Fermion2Qubit.hf_time = ham2_ActiveSpace.hf_time ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time ham3_Fermion2Qubit.f2q_time = tstop - tstart @@ -379,67 +379,32 @@ def get_ham2(state): # TODO: Python errors get a little weird if you have an exception inside an exception. So # instead the exception clause should _only_ flag that we're going to recompute versus # load, and the actual recomputing should be outside of the except clause. - except FileNotFoundError: - ham2_ActiveSpace = None - - requested_threshold = state.config_hamiltonian.coefficient_threshold - if ham2_ActiveSpace is not None: - # Pickles created before threshold configurability necessarily used - # OpenFermion's historical 1e-8 default. - cached_threshold = getattr( - ham2_ActiveSpace, - "coefficient_threshold", - DEFAULT_COEFFICIENT_THRESHOLD, - ) - if cached_threshold == requested_threshold: - ham2_ActiveSpace.coefficient_threshold = cached_threshold - if not hasattr(ham2_ActiveSpace, "threshold_stats"): - one_after = int(np.count_nonzero(ham2_ActiveSpace.one_body_tensor)) - two_after = int(np.count_nonzero(ham2_ActiveSpace.two_body_tensor)) - ham2_ActiveSpace.threshold_stats = { - "one_body_before": None, - "one_body_after": one_after, - "one_body_removed": None, - "two_body_before": None, - "two_body_after": two_after, - "two_body_removed": None, - } - state.log(' '.join([ - f"Loaded \"{ham2_filename}\" with coefficient threshold", - f"{cached_threshold:.6g}.", - "Continuing from after the active space is applied."])) - return ham2_ActiveSpace - + except FileNotFoundError as err: state.log(' '.join([ - f"Cached \"{ham2_filename}\" used coefficient threshold", - f"{cached_threshold:.6g}, but {requested_threshold:.6g} was requested.", - "Recomputing the active-space tensors."])) + f"Could not load \"{ham2_filename}\".", + f"Trying to load \"{state.filename_ham1()}\"."])) + # Get ham1_HartreeFock (by loading or by recomputing, depending on data availability) + ham1_HartreeFock = get_ham1(state) + # Recompute ham2_ActiveSpace from ham1_HartreeFock + state.log("Apply active space.") + ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) + state.log(f"Pickle to \"{ham2_filename}\" file.") + # Save ham2_ActiveSpace for later re-use + with open(ham2_filename, 'wb') as ham2_file: + pickle.dump(ham2_ActiveSpace, ham2_file) + # Save the one-body and two-body tensors using numpy's `save` function. These files can be + # loaded using numpy's `load` function. + t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" + np.savez_compressed(t_filename, + constant=ham2_ActiveSpace.constant, + one_body=ham2_ActiveSpace.one_body_tensor, + two_body=ham2_ActiveSpace.two_body_tensor) + return ham2_ActiveSpace else: - state.log(f"Could not load \"{ham2_filename}\".") - - state.log(f"Trying to load \"{state.filename_ham1()}\".") - ham1_HartreeFock = get_ham1(state) - state.log("Apply active space.") - ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) - state.log(f"Pickle to \"{ham2_filename}\" file.") - with open(ham2_filename, 'wb') as ham2_file: - pickle.dump(ham2_ActiveSpace, ham2_file) - - # Store threshold provenance alongside the tensors. - stats = ham2_ActiveSpace.threshold_stats - t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" - np.savez_compressed( - t_filename, - constant=ham2_ActiveSpace.constant, - one_body=ham2_ActiveSpace.one_body_tensor, - two_body=ham2_ActiveSpace.two_body_tensor, - coefficient_threshold=ham2_ActiveSpace.coefficient_threshold, - one_body_entries_before_threshold=stats["one_body_before"], - one_body_entries_after_threshold=stats["one_body_after"], - two_body_entries_before_threshold=stats["two_body_before"], - two_body_entries_after_threshold=stats["two_body_after"], - ) - return ham2_ActiveSpace + state.log(' '.join([ + f"Loaded \"{ham2_filename}\".", + "Continuing from after the active space is applied."])) + return ham2_ActiveSpace # ------------------------------------------------------------------------------------------------- @@ -497,27 +462,16 @@ def compute_metadata(state, ham3_Fermion2Qubit): state.metadata["interatomic separation (angstroms)"] = ham3_Fermion2Qubit.separation # fermion-to-qubit transformation state.metadata["fermion-to-qubit operator mapping"] = ham3_Fermion2Qubit.f2q_mapping - state.metadata["coefficient threshold (Hartrees)"] = ( - ham3_Fermion2Qubit.coefficient_threshold + state.metadata["spin-orbital coefficient threshold (Hartrees)"] = ( + state.config_hamiltonian.coefficient_threshold ) - threshold_stats = ham3_Fermion2Qubit.threshold_stats - for tensor_name in ("one_body", "two_body"): - before = threshold_stats[f"{tensor_name}_before"] - after = threshold_stats[f"{tensor_name}_after"] - if before is not None: - state.metadata[ - f"{tensor_name.replace('_', '-')} tensor entries before thresholding" - ] = before - state.metadata[ - f"{tensor_name.replace('_', '-')} tensor entries after thresholding" - ] = after # number of terms in sum of Pauli strings state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms) # one-norm of sum of Pauli strings one_norm = sum(abs(coefficient) for coefficient in ham3_Fermion2Qubit.terms.values()) state.metadata["one-norm of sum of Pauli strings (Hartrees)"] = one_norm # eigenvalue bounds (see Trotter Workflow document for where this comes from) - E0 = ham3_Fermion2Qubit.terms.get(tuple(), 0.0) + E0 = ham3_Fermion2Qubit.terms[tuple()] dE = one_norm - abs(E0) state.metadata["eigenvalue lower bound from shifted one-norm (Hartrees)"] = E0 - dE state.metadata["eigenvalue upper bound from shifted one-norm (Hartrees)"] = E0 + dE diff --git a/hamiltonian_generator/thresholding.py b/hamiltonian_generator/thresholding.py index ef98bd84..d3502df4 100644 --- a/hamiltonian_generator/thresholding.py +++ b/hamiltonian_generator/thresholding.py @@ -1,17 +1,9 @@ -"""Small QHAT-owned hotfix for OpenFermion's fixed coefficient cutoff.""" +"""QHAT-owned spatial-to-spin conversion with a configurable cutoff.""" -from collections import defaultdict import math import numbers import numpy as np -from openfermion import ( - FermionOperator, - InteractionOperator, - QubitOperator, - bravyi_kitaev, - jordan_wigner, -) DEFAULT_COEFFICIENT_THRESHOLD = 1.0e-8 @@ -26,21 +18,6 @@ def validate_coefficient_threshold(value) -> float: return threshold -def _validate_spatial_integrals(one_body_integrals, two_body_integrals): - one_body = np.asarray(one_body_integrals) - two_body = np.asarray(two_body_integrals) - - if one_body.ndim != 2 or one_body.shape[0] != one_body.shape[1]: - raise ValueError("one-body spatial integrals must be a square rank-2 array") - expected_two_body_shape = (one_body.shape[0],) * 4 - if two_body.shape != expected_two_body_shape: - raise ValueError( - "two-body spatial integrals must have shape " - f"{expected_two_body_shape}, got {two_body.shape}" - ) - return one_body, two_body - - def spinorb_from_spatial( one_body_integrals, two_body_integrals, @@ -55,127 +32,20 @@ def spinorb_from_spatial( """ threshold = validate_coefficient_threshold(coefficient_threshold) - one_spatial, two_spatial = _validate_spatial_integrals( - one_body_integrals, two_body_integrals - ) - n_qubits = 2 * one_spatial.shape[0] + n_qubits = 2 * one_body_integrals.shape[0] - one_spin = np.zeros((n_qubits, n_qubits), dtype=one_spatial.dtype) - two_spin = np.zeros((n_qubits,) * 4, dtype=two_spatial.dtype) + one_spin = np.zeros((n_qubits, n_qubits)) + two_spin = np.zeros((n_qubits,) * 4) - one_spin[0::2, 0::2] = one_spatial - one_spin[1::2, 1::2] = one_spatial + one_spin[0::2, 0::2] = one_body_integrals + one_spin[1::2, 1::2] = one_body_integrals - two_spin[0::2, 1::2, 1::2, 0::2] = two_spatial - two_spin[1::2, 0::2, 0::2, 1::2] = two_spatial - two_spin[0::2, 0::2, 0::2, 0::2] = two_spatial - two_spin[1::2, 1::2, 1::2, 1::2] = two_spatial + two_spin[0::2, 1::2, 1::2, 0::2] = two_body_integrals + two_spin[1::2, 0::2, 0::2, 1::2] = two_body_integrals + two_spin[0::2, 0::2, 0::2, 0::2] = two_body_integrals + two_spin[1::2, 1::2, 1::2, 1::2] = two_body_integrals - if threshold > 0.0: - one_spin[np.abs(one_spin) < threshold] = 0.0 - two_spin[np.abs(two_spin) < threshold] = 0.0 + one_spin[np.abs(one_spin) < threshold] = 0.0 + two_spin[np.abs(two_spin) < threshold] = 0.0 return one_spin, two_spin - - -def molecular_hamiltonian_from_spatial( - constant, - one_body_integrals, - two_body_integrals, - coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, -): - """Construct an ``InteractionOperator`` and record cutoff provenance.""" - - threshold = validate_coefficient_threshold(coefficient_threshold) - one_spatial, two_spatial = _validate_spatial_integrals( - one_body_integrals, two_body_integrals - ) - one_before = 2 * int(np.count_nonzero(one_spatial)) - two_before = 4 * int(np.count_nonzero(two_spatial)) - one_spin, two_spin = spinorb_from_spatial(one_spatial, two_spatial, threshold) - - operator = InteractionOperator(constant, one_spin, 0.5 * two_spin) - operator.coefficient_threshold = threshold - operator.threshold_stats = { - "one_body_before": one_before, - "one_body_after": int(np.count_nonzero(one_spin)), - "one_body_removed": one_before - int(np.count_nonzero(one_spin)), - "two_body_before": two_before, - "two_body_after": int(np.count_nonzero(two_spin)), - "two_body_removed": two_before - int(np.count_nonzero(two_spin)), - } - return operator - - -def get_molecular_hamiltonian( - molecule, - occupied_indices=None, - active_indices=None, - coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, -): - """QHAT-local counterpart of ``MolecularData.get_molecular_hamiltonian``.""" - - if occupied_indices is None and active_indices is None: - one_body_integrals, two_body_integrals = molecule.get_integrals() - constant = molecule.nuclear_repulsion - else: - core_adjustment, one_body_integrals, two_body_integrals = ( - molecule.get_active_space_integrals(occupied_indices, active_indices) - ) - constant = molecule.nuclear_repulsion + core_adjustment - - return molecular_hamiltonian_from_spatial( - constant, - one_body_integrals, - two_body_integrals, - coefficient_threshold, - ) - - -def map_interaction_operator( - operator, - mapping_name, - coefficient_threshold=DEFAULT_COEFFICIENT_THRESHOLD, -): - """Map without changing OpenFermion's global ``EQ_TOLERANCE``. - - Each unit-coefficient monomial safely survives OpenFermion's internal - cutoff. QHAT restores the real coefficient, sums matching Pauli terms, - and applies the requested cutoff once at the end. - """ - - threshold = validate_coefficient_threshold(coefficient_threshold) - if mapping_name not in ("jordan-wigner", "bravyi-kitaev"): - raise ValueError(f"invalid fermion-to-qubit mapping {mapping_name!r}") - - contributions = defaultdict(list) - for fermion_term in operator: - coefficient = operator[fermion_term] - unit_term = FermionOperator(fermion_term, 1.0) - if mapping_name == "jordan-wigner": - mapped_unit = jordan_wigner(unit_term) - else: - mapped_unit = bravyi_kitaev(unit_term, n_qubits=operator.n_qubits) - for pauli_term, unit_coefficient in mapped_unit.terms.items(): - contributions[pauli_term].append(coefficient * unit_coefficient) - - terms = {} - for pauli_term, values in contributions.items(): - coefficient = complex( - math.fsum(complex(value).real for value in values), - math.fsum(complex(value).imag for value in values), - ) - # Molecular Hamiltonians are Hermitian. Permit only machine-scale - # imaginary residue from floating-point tensor symmetries. - scale = math.fsum(abs(value) for value in values) - imag_tolerance = 128.0 * np.finfo(float).eps * max(1.0, scale) - if abs(coefficient.imag) > imag_tolerance: - raise ValueError(f"non-Hermitian Pauli coefficient {coefficient!r}") - coefficient = coefficient.real - if coefficient != 0.0 and abs(coefficient) >= threshold: - terms[pauli_term] = coefficient - - result = QubitOperator() - # Do not use += here; it would reintroduce OpenFermion's fixed cutoff. - result.terms = terms - return result From 1e81e15ea16e8e3d4b7c277222f2ca85c9eed6f1 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Thu, 16 Jul 2026 10:20:20 -0600 Subject: [PATCH 04/16] Remove Threshold Hotfix.md and others. --- hamiltonian_generator/README.md | 2 +- hamiltonian_generator/THRESHOLDING_HOTFIX.md | 42 -------------------- 2 files changed, 1 insertion(+), 43 deletions(-) delete mode 100644 hamiltonian_generator/THRESHOLDING_HOTFIX.md diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index ec6c8317..c2787dbc 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -85,7 +85,7 @@ hamiltonian.coefficient_threshold = 0.0 The value must be finite and non-negative. This does not change OpenFermion's standard JW/BK mapping, which can still discard small coefficients. Changing the value also does not invalidate an existing active-space pickle; remove that pickle or use a different `general.file_stub` when -comparing thresholds. See [THRESHOLDING_HOTFIX.md](THRESHOLDING_HOTFIX.md) for details. +comparing thresholds. ## Generated Files diff --git a/hamiltonian_generator/THRESHOLDING_HOTFIX.md b/hamiltonian_generator/THRESHOLDING_HOTFIX.md deleted file mode 100644 index 2221eb96..00000000 --- a/hamiltonian_generator/THRESHOLDING_HOTFIX.md +++ /dev/null @@ -1,42 +0,0 @@ -# Issue #31 tensor-stage threshold hotfix - -This is a narrow QHAT-only change. OpenFermion is not edited or monkey-patched. - -## What changed - -The Hamiltonian configuration now accepts a finite, non-negative cutoff: - -```python -hamiltonian.coefficient_threshold = 1e-8 # default -hamiltonian.coefficient_threshold = 0.0 # no tensor-stage filtering -``` - -QHAT uses a local `spinorb_from_spatial` function when it expands active-space spatial integrals -into spin-orbital tensors. The rest of the generator stays on its original path: - -- OpenFermion's standard Jordan-Wigner and Bravyi-Kitaev functions are unchanged. -- The active-space pickle behavior is unchanged. -- The `.tensors.npz` format remains `constant`, `one_body`, and `two_body`. -- Metadata records only the configured spin-orbital threshold. - -## Limitations - -OpenFermion's mapping functions can still remove small coefficients. Setting this option to -`0.0` therefore disables filtering only during tensor construction, not in the final Pauli sum. -Pauli-level thresholding is a separate follow-up. - -An existing `[stub]_[astag].pickle` is reused without checking its threshold. When changing the -setting, remove that active-space pickle or use a new `general.file_stub`. Removing only the -`.tensors.npz` file is not enough. - -## Li2 check - -For the bundled five-spatial-orbital Li2 case: - -| Tensor threshold | One-body entries | Two-body entries | Stock JW terms | Stock BK terms | -| ---: | ---: | ---: | ---: | ---: | -| `1e-8` | 14 | 436 | 156 | 156 | -| `0.0` | 50 | 2500 | 156 | 156 | - -The tensor counts confirm that the setting works at the intended stage. The unchanged mapped -term counts confirm that final Pauli thresholding remains outside this hotfix. From 661bf09dff05e1e84425cef8db19781b177dda33 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 17 Jul 2026 08:13:14 -0600 Subject: [PATCH 05/16] Fix threshold cache provenance and add tests --- .github/workflows/python-app.yml | 6 + hamiltonian_generator/README.md | 17 +- .../diatomic_lithium_as-004-006_jw.dat | 1 + hamiltonian_generator/hamgen.py | 79 +++-- .../tests/test_thresholding.py | 285 ++++++++++++++++++ hamiltonian_generator/thresholding.py | 7 +- 6 files changed, 363 insertions(+), 32 deletions(-) create mode 100644 hamiltonian_generator/tests/test_thresholding.py diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index ab74dfaa..30129788 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -66,6 +66,12 @@ jobs: script: hamgen.py config_file: '' + - name: hamgen-unit-tests + directory: hamiltonian_generator + script: pytest + config_file: '' + config_setup: '' + - name: analysis-trotter directory: analysis script: driver.py diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index c2787dbc..0d1df022 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -74,18 +74,21 @@ Currently allowed values are - for Bravyi-Kitaev: "Bravyi-Kitaev", "Bravyi Kitaev", "BK", "bravyi-kitaev", "bravyi\_kitaev", "bravyi kitaev", or "bk" -**`hamiltonian.coefficient_threshold`** controls the absolute cutoff used when QHAT expands -active-space spatial integrals into spin-orbital tensors. Its default is `1e-8`; set it to `0.0` -to disable filtering at this stage: +**`hamiltonian.coefficient_threshold`** controls the absolute cutoff, in Hartrees, used when QHAT +expands active-space spatial integrals into spin-orbital tensors. QHAT sets tensor coefficients to +zero when `abs(value) < coefficient_threshold`, so coefficients exactly equal to the threshold are +retained. Its default is `1e-8`; set it to `0.0` to disable filtering at this stage: ```python hamiltonian.coefficient_threshold = 0.0 ``` -The value must be finite and non-negative. This does not change OpenFermion's standard JW/BK -mapping, which can still discard small coefficients. Changing the value also does not invalidate -an existing active-space pickle; remove that pickle or use a different `general.file_stub` when -comparing thresholds. +The value must be finite and non-negative. This setting affects only tensor construction; it does +not change OpenFermion's standard JW/BK mapping, which can still discard small coefficients. +Active-space pickles record the threshold used to construct their tensors. If the configured value +does not match `[filestub]_[astag].pickle`, QHAT recomputes that stage and replaces both the pickle +and `[filestub]_[astag].tensors.npz`. Pickles created before this setting existed are interpreted as +using the historical `1e-8` cutoff. ## Generated Files diff --git a/hamiltonian_generator/diatomic_lithium_as-004-006_jw.dat b/hamiltonian_generator/diatomic_lithium_as-004-006_jw.dat index 62092507..13ca1f11 100644 --- a/hamiltonian_generator/diatomic_lithium_as-004-006_jw.dat +++ b/hamiltonian_generator/diatomic_lithium_as-004-006_jw.dat @@ -7,6 +7,7 @@ # basis set = sto-3g # interatomic separation (angstroms) = 2.0 # fermion-to-qubit operator mapping = jordan-wigner +# spin-orbital coefficient threshold (Hartrees) = 1e-08 # number of terms in sum of Pauli strings = 156 # one-norm of sum of Pauli strings (Hartrees) = 18.125380876421513 # eigenvalue lower bound from shifted one-norm (Hartrees) = -18.125380876421513 diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 9f16e9d7..72a6d305 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -20,7 +20,11 @@ HamiltonianConfiguration, State, ) -from qhat.hamiltonian_generator.thresholding import spinorb_from_spatial +from qhat.hamiltonian_generator.thresholding import ( + DEFAULT_COEFFICIENT_THRESHOLD, + spinorb_from_spatial, + validate_coefficient_threshold, +) logger = logging.getLogger(__name__) @@ -270,6 +274,7 @@ def apply_active_space(state, ham1_HartreeFock): ham2_ActiveSpace.separation = ham1_HartreeFock.separation ham2_ActiveSpace.hf_time = ham1_HartreeFock.hf_time ham2_ActiveSpace.as_time = tstop - tstart + ham2_ActiveSpace.coefficient_threshold = state.config_hamiltonian.coefficient_threshold return ham2_ActiveSpace # ------------------------------------------------------------------------------------------------- @@ -296,6 +301,7 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): ham3_Fermion2Qubit.hf_time = ham2_ActiveSpace.hf_time ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time ham3_Fermion2Qubit.f2q_time = tstop - tstart + ham3_Fermion2Qubit.coefficient_threshold = ham2_ActiveSpace.coefficient_threshold # Return result return ham3_Fermion2Qubit @@ -379,32 +385,57 @@ def get_ham2(state): # TODO: Python errors get a little weird if you have an exception inside an exception. So # instead the exception clause should _only_ flag that we're going to recompute versus # load, and the actual recomputing should be outside of the except clause. - except FileNotFoundError as err: + except FileNotFoundError: state.log(' '.join([ f"Could not load \"{ham2_filename}\".", f"Trying to load \"{state.filename_ham1()}\"."])) - # Get ham1_HartreeFock (by loading or by recomputing, depending on data availability) - ham1_HartreeFock = get_ham1(state) - # Recompute ham2_ActiveSpace from ham1_HartreeFock - state.log("Apply active space.") - ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) - state.log(f"Pickle to \"{ham2_filename}\" file.") - # Save ham2_ActiveSpace for later re-use - with open(ham2_filename, 'wb') as ham2_file: - pickle.dump(ham2_ActiveSpace, ham2_file) - # Save the one-body and two-body tensors using numpy's `save` function. These files can be - # loaded using numpy's `load` function. - t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" - np.savez_compressed(t_filename, - constant=ham2_ActiveSpace.constant, - one_body=ham2_ActiveSpace.one_body_tensor, - two_body=ham2_ActiveSpace.two_body_tensor) - return ham2_ActiveSpace else: - state.log(' '.join([ - f"Loaded \"{ham2_filename}\".", - "Continuing from after the active space is applied."])) - return ham2_ActiveSpace + cached_threshold_value = getattr( + ham2_ActiveSpace, + "coefficient_threshold", + DEFAULT_COEFFICIENT_THRESHOLD, + ) + configured_threshold = state.config_hamiltonian.coefficient_threshold + try: + cached_threshold = validate_coefficient_threshold(cached_threshold_value) + except (TypeError, ValueError): + state.log(' '.join([ + f"Loaded \"{ham2_filename}\" with invalid coefficient threshold", + f"provenance of type {type(cached_threshold_value).__name__}.", + "Recomputing the active space."])) + state.log(f"Trying to load \"{state.filename_ham1()}\".") + else: + if cached_threshold == configured_threshold: + # Active-space pickles created before the threshold was configurable used + # OpenFermion's historical 1e-8 cutoff and have no explicit provenance field. + ham2_ActiveSpace.coefficient_threshold = cached_threshold + state.log(' '.join([ + f"Loaded \"{ham2_filename}\".", + "Continuing from after the active space is applied."])) + return ham2_ActiveSpace + state.log(' '.join([ + f"Loaded \"{ham2_filename}\" with coefficient threshold {cached_threshold},", + f"but the configured threshold is {configured_threshold}.", + "Recomputing the active space."])) + state.log(f"Trying to load \"{state.filename_ham1()}\".") + + # Get ham1_HartreeFock (by loading or by recomputing, depending on data availability) + ham1_HartreeFock = get_ham1(state) + # Recompute ham2_ActiveSpace from ham1_HartreeFock + state.log("Apply active space.") + ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) + state.log(f"Pickle to \"{ham2_filename}\" file.") + # Save ham2_ActiveSpace for later re-use + with open(ham2_filename, 'wb') as ham2_file: + pickle.dump(ham2_ActiveSpace, ham2_file) + # Save the one-body and two-body tensors using numpy's `save` function. These files can be + # loaded using numpy's `load` function. + t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" + np.savez_compressed(t_filename, + constant=ham2_ActiveSpace.constant, + one_body=ham2_ActiveSpace.one_body_tensor, + two_body=ham2_ActiveSpace.two_body_tensor) + return ham2_ActiveSpace # ------------------------------------------------------------------------------------------------- @@ -463,7 +494,7 @@ def compute_metadata(state, ham3_Fermion2Qubit): # fermion-to-qubit transformation state.metadata["fermion-to-qubit operator mapping"] = ham3_Fermion2Qubit.f2q_mapping state.metadata["spin-orbital coefficient threshold (Hartrees)"] = ( - state.config_hamiltonian.coefficient_threshold + ham3_Fermion2Qubit.coefficient_threshold ) # number of terms in sum of Pauli strings state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms) diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py new file mode 100644 index 00000000..22b11025 --- /dev/null +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -0,0 +1,285 @@ +"""Regression tests for QHAT's configurable coefficient threshold.""" + +import pickle +from pathlib import Path + +import numpy as np +import pytest +from openfermion import InteractionOperator +from openfermion.chem.molecular_data import ( + spinorb_from_spatial as openfermion_spinorb_from_spatial, +) + +from qhat.hamiltonian_generator import hamgen +from qhat.hamiltonian_generator.hamgen_types import HamiltonianConfiguration +from qhat.hamiltonian_generator.thresholding import ( + DEFAULT_COEFFICIENT_THRESHOLD, + spinorb_from_spatial, +) + + +HAMGEN_DIRECTORY = Path(__file__).resolve().parents[1] + + +class _ActiveSpaceState: + def __init__(self, threshold=DEFAULT_COEFFICIENT_THRESHOLD): + self.config_hamiltonian = HamiltonianConfiguration() + self.config_hamiltonian.num_active_occupied = 4 + self.config_hamiltonian.num_active_vacant = 6 + self.config_hamiltonian.coefficient_threshold = threshold + + def log(self, *args, **kwargs): + pass + + def log_verbose(self, *args, **kwargs): + pass + + +class _CacheState: + def __init__(self, cache_path, threshold): + self.cache_path = cache_path + self.config_hamiltonian = HamiltonianConfiguration() + self.config_hamiltonian.coefficient_threshold = threshold + + def filename_ham1(self): + return str(self.cache_path.with_name("hartree-fock.pickle")) + + def filename_ham2(self): + return str(self.cache_path) + + def log(self, *args, **kwargs): + pass + + +def _interaction_operator(constant, threshold=None): + operator = InteractionOperator( + constant, + np.diag([constant, 0.0]), + np.zeros((2, 2, 2, 2)), + ) + operator.hf_energy = -1.0 + operator.asmeta = { + "n_frz_occ_so": 0, + "n_act_occ_so": 1, + "n_act_vac_so": 1, + "n_frz_vac_so": 0, + } + operator.basis = "test-basis" + operator.separation = 1.0 + operator.hf_time = 0.0 + operator.as_time = 0.0 + if threshold is not None: + operator.coefficient_threshold = threshold + return operator + + +def test_configuration_uses_openfermion_default_threshold(): + config = HamiltonianConfiguration() + assert config.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD == 1.0e-8 + + +@pytest.mark.parametrize("value", [0, 1.0e-12, np.float64(1.0e-8)]) +def test_configuration_accepts_non_negative_real_thresholds(value): + config = HamiltonianConfiguration() + config.coefficient_threshold = value + assert config.coefficient_threshold == float(value) + + +@pytest.mark.parametrize( + ("value", "exception"), + [ + (-1.0e-8, ValueError), + (float("nan"), ValueError), + (float("inf"), ValueError), + (-float("inf"), ValueError), + pytest.param(10**10000, ValueError, id="overflowing-real"), + (True, TypeError), + ("1e-8", TypeError), + (None, TypeError), + ], +) +def test_configuration_rejects_invalid_thresholds(value, exception): + config = HamiltonianConfiguration() + with pytest.raises(exception): + config.coefficient_threshold = value + + +def test_spinorb_default_parity_strict_boundary_and_zero_threshold(): + threshold = DEFAULT_COEFFICIENT_THRESHOLD + below_threshold = threshold / 2.0 + one_body = np.array( + [ + [below_threshold, threshold], + [-2.0 * threshold, 0.0], + ] + ) + two_body = np.zeros((2, 2, 2, 2)) + two_body[0, 0, 0, 0] = below_threshold + two_body[0, 0, 0, 1] = threshold + two_body[1, 0, 0, 1] = -2.0 * threshold + + expected_one, expected_two = openfermion_spinorb_from_spatial(one_body, two_body) + default_one, default_two = spinorb_from_spatial(one_body, two_body) + zero_one, zero_two = spinorb_from_spatial(one_body, two_body, 0.0) + + np.testing.assert_array_equal(default_one, expected_one) + np.testing.assert_array_equal(default_two, expected_two) + + assert default_one[0, 0] == 0.0 + assert default_one[0, 2] == threshold + assert default_two[0, 1, 1, 0] == 0.0 + assert default_two[0, 1, 1, 2] == threshold + + assert zero_one[0, 0] == below_threshold + assert zero_two[0, 1, 1, 0] == below_threshold + + +def test_active_space_matches_openfermion_and_respects_zero_threshold(): + with open(HAMGEN_DIRECTORY / "diatomic_lithium.pickle", "rb") as file: + molecule = pickle.load(file) + + default_operator = hamgen.apply_active_space(_ActiveSpaceState(), molecule) + expected_operator = molecule.get_molecular_hamiltonian( + occupied_indices=range(1), + active_indices=range(1, 6), + ) + + assert default_operator.constant == expected_operator.constant + np.testing.assert_array_equal( + default_operator.one_body_tensor, + expected_operator.one_body_tensor, + ) + np.testing.assert_array_equal( + default_operator.two_body_tensor, + expected_operator.two_body_tensor, + ) + assert default_operator.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD + + zero_operator = hamgen.apply_active_space(_ActiveSpaceState(0.0), molecule) + assert zero_operator.coefficient_threshold == 0.0 + assert np.count_nonzero(zero_operator.one_body_tensor) > np.count_nonzero( + default_operator.one_body_tensor + ) + assert np.count_nonzero(zero_operator.two_body_tensor) > np.count_nonzero( + default_operator.two_body_tensor + ) + + +def test_mapping_and_metadata_use_effective_operator_threshold(): + state = _ActiveSpaceState(0.0) + state.metadata = {} + active_operator = _interaction_operator( + 1.0, + threshold=DEFAULT_COEFFICIENT_THRESHOLD, + ) + + mapped_operator = hamgen.map_fermions_to_qubits(state, active_operator) + hamgen.compute_metadata(state, mapped_operator) + + assert mapped_operator.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD + assert mapped_operator.coefficient_threshold != state.config_hamiltonian.coefficient_threshold + assert state.metadata["spin-orbital coefficient threshold (Hartrees)"] == ( + DEFAULT_COEFFICIENT_THRESHOLD + ) + + +def test_active_space_cache_tracks_threshold_and_recomputes_on_mismatch( + tmp_path, + monkeypatch, +): + cache_path = tmp_path / "active-space.pickle" + legacy_operator = _interaction_operator(1.0) + with open(cache_path, "wb") as file: + pickle.dump(legacy_operator, file) + + calls = [] + recomputed_operator = _interaction_operator(2.0, threshold=0.0) + + def fake_get_ham1(state): + calls.append("get_ham1") + return object() + + def fake_apply_active_space(state, molecule): + calls.append("apply_active_space") + return recomputed_operator + + monkeypatch.setattr(hamgen, "get_ham1", fake_get_ham1) + monkeypatch.setattr(hamgen, "apply_active_space", fake_apply_active_space) + + default_state = _CacheState(cache_path, DEFAULT_COEFFICIENT_THRESHOLD) + loaded_legacy = hamgen.get_ham2(default_state) + assert loaded_legacy.constant == legacy_operator.constant + assert loaded_legacy.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD + assert calls == [] + + zero_state = _CacheState(cache_path, 0.0) + rebuilt = hamgen.get_ham2(zero_state) + assert rebuilt is recomputed_operator + assert calls == ["get_ham1", "apply_active_space"] + + with open(cache_path, "rb") as file: + persisted = pickle.load(file) + assert persisted.coefficient_threshold == 0.0 + assert persisted.constant == recomputed_operator.constant + with np.load(tmp_path / "active-space.tensors.npz") as tensors: + assert set(tensors.files) == {"constant", "one_body", "two_body"} + assert tensors["constant"] == recomputed_operator.constant + np.testing.assert_array_equal( + tensors["one_body"], + recomputed_operator.one_body_tensor, + ) + np.testing.assert_array_equal( + tensors["two_body"], + recomputed_operator.two_body_tensor, + ) + + loaded_match = hamgen.get_ham2(zero_state) + assert loaded_match.constant == recomputed_operator.constant + assert loaded_match.coefficient_threshold == 0.0 + assert calls == ["get_ham1", "apply_active_space"] + + +@pytest.mark.parametrize( + "invalid_threshold", + [ + False, + float("nan"), + np.array([0.0]), + pytest.param(10**10000, id="overflowing-real"), + ], +) +def test_active_space_cache_recomputes_invalid_threshold_provenance( + tmp_path, + monkeypatch, + invalid_threshold, +): + cache_path = tmp_path / "active-space.pickle" + cached_operator = _interaction_operator(1.0, threshold=invalid_threshold) + with open(cache_path, "wb") as file: + pickle.dump(cached_operator, file) + + recomputed_operator = _interaction_operator( + 2.0, + threshold=DEFAULT_COEFFICIENT_THRESHOLD, + ) + calls = [] + + def fake_get_ham1(state): + calls.append("get_ham1") + return object() + + def fake_apply_active_space(state, molecule): + calls.append("apply_active_space") + return recomputed_operator + + monkeypatch.setattr(hamgen, "get_ham1", fake_get_ham1) + monkeypatch.setattr(hamgen, "apply_active_space", fake_apply_active_space) + + state = _CacheState(cache_path, DEFAULT_COEFFICIENT_THRESHOLD) + rebuilt = hamgen.get_ham2(state) + + assert rebuilt is recomputed_operator + assert calls == ["get_ham1", "apply_active_space"] + with open(cache_path, "rb") as file: + persisted = pickle.load(file) + assert persisted.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD diff --git a/hamiltonian_generator/thresholding.py b/hamiltonian_generator/thresholding.py index d3502df4..2092d54d 100644 --- a/hamiltonian_generator/thresholding.py +++ b/hamiltonian_generator/thresholding.py @@ -12,7 +12,12 @@ def validate_coefficient_threshold(value) -> float: if isinstance(value, bool) or not isinstance(value, numbers.Real): raise TypeError("coefficient threshold must be a real number") - threshold = float(value) + try: + threshold = float(value) + except OverflowError as error: + raise ValueError( + "coefficient threshold must be finite and non-negative" + ) from error if not math.isfinite(threshold) or threshold < 0.0: raise ValueError("coefficient threshold must be finite and non-negative") return threshold From 17a6403857b93c6481d1331af57066d1f003a646 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 17 Jul 2026 08:17:52 -0600 Subject: [PATCH 06/16] Simplify threshold cache validation --- hamiltonian_generator/README.md | 5 +- hamiltonian_generator/hamgen.py | 75 ++++++++----------- .../tests/test_thresholding.py | 68 +++-------------- hamiltonian_generator/thresholding.py | 7 +- 4 files changed, 43 insertions(+), 112 deletions(-) diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index 0d1df022..4994f3d0 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -86,8 +86,9 @@ hamiltonian.coefficient_threshold = 0.0 The value must be finite and non-negative. This setting affects only tensor construction; it does not change OpenFermion's standard JW/BK mapping, which can still discard small coefficients. Active-space pickles record the threshold used to construct their tensors. If the configured value -does not match `[filestub]_[astag].pickle`, QHAT recomputes that stage and replaces both the pickle -and `[filestub]_[astag].tensors.npz`. Pickles created before this setting existed are interpreted as +does not match `[filestub]_[astag].pickle`, QHAT stops rather than silently reusing incompatible +tensors. Remove that active-space pickle and rerun; QHAT will rebuild it and overwrite +`[filestub]_[astag].tensors.npz`. Pickles created before this setting existed are interpreted as using the historical `1e-8` cutoff. ## Generated Files diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 72a6d305..9d17c1d3 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -23,7 +23,6 @@ from qhat.hamiltonian_generator.thresholding import ( DEFAULT_COEFFICIENT_THRESHOLD, spinorb_from_spatial, - validate_coefficient_threshold, ) logger = logging.getLogger(__name__) @@ -301,7 +300,6 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): ham3_Fermion2Qubit.hf_time = ham2_ActiveSpace.hf_time ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time ham3_Fermion2Qubit.f2q_time = tstop - tstart - ham3_Fermion2Qubit.coefficient_threshold = ham2_ActiveSpace.coefficient_threshold # Return result return ham3_Fermion2Qubit @@ -385,57 +383,44 @@ def get_ham2(state): # TODO: Python errors get a little weird if you have an exception inside an exception. So # instead the exception clause should _only_ flag that we're going to recompute versus # load, and the actual recomputing should be outside of the except clause. - except FileNotFoundError: + except FileNotFoundError as err: state.log(' '.join([ f"Could not load \"{ham2_filename}\".", f"Trying to load \"{state.filename_ham1()}\"."])) + # Get ham1_HartreeFock (by loading or by recomputing, depending on data availability) + ham1_HartreeFock = get_ham1(state) + # Recompute ham2_ActiveSpace from ham1_HartreeFock + state.log("Apply active space.") + ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) + state.log(f"Pickle to \"{ham2_filename}\" file.") + # Save ham2_ActiveSpace for later re-use + with open(ham2_filename, 'wb') as ham2_file: + pickle.dump(ham2_ActiveSpace, ham2_file) + # Save the one-body and two-body tensors using numpy's `save` function. These files can be + # loaded using numpy's `load` function. + t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" + np.savez_compressed(t_filename, + constant=ham2_ActiveSpace.constant, + one_body=ham2_ActiveSpace.one_body_tensor, + two_body=ham2_ActiveSpace.two_body_tensor) + return ham2_ActiveSpace else: - cached_threshold_value = getattr( + cached_threshold = getattr( ham2_ActiveSpace, "coefficient_threshold", DEFAULT_COEFFICIENT_THRESHOLD, ) configured_threshold = state.config_hamiltonian.coefficient_threshold - try: - cached_threshold = validate_coefficient_threshold(cached_threshold_value) - except (TypeError, ValueError): - state.log(' '.join([ - f"Loaded \"{ham2_filename}\" with invalid coefficient threshold", - f"provenance of type {type(cached_threshold_value).__name__}.", - "Recomputing the active space."])) - state.log(f"Trying to load \"{state.filename_ham1()}\".") - else: - if cached_threshold == configured_threshold: - # Active-space pickles created before the threshold was configurable used - # OpenFermion's historical 1e-8 cutoff and have no explicit provenance field. - ham2_ActiveSpace.coefficient_threshold = cached_threshold - state.log(' '.join([ - f"Loaded \"{ham2_filename}\".", - "Continuing from after the active space is applied."])) - return ham2_ActiveSpace - state.log(' '.join([ - f"Loaded \"{ham2_filename}\" with coefficient threshold {cached_threshold},", - f"but the configured threshold is {configured_threshold}.", - "Recomputing the active space."])) - state.log(f"Trying to load \"{state.filename_ham1()}\".") - - # Get ham1_HartreeFock (by loading or by recomputing, depending on data availability) - ham1_HartreeFock = get_ham1(state) - # Recompute ham2_ActiveSpace from ham1_HartreeFock - state.log("Apply active space.") - ham2_ActiveSpace = apply_active_space(state, ham1_HartreeFock) - state.log(f"Pickle to \"{ham2_filename}\" file.") - # Save ham2_ActiveSpace for later re-use - with open(ham2_filename, 'wb') as ham2_file: - pickle.dump(ham2_ActiveSpace, ham2_file) - # Save the one-body and two-body tensors using numpy's `save` function. These files can be - # loaded using numpy's `load` function. - t_filename = ham2_filename[:ham2_filename.rfind('.')] + ".tensors.npz" - np.savez_compressed(t_filename, - constant=ham2_ActiveSpace.constant, - one_body=ham2_ActiveSpace.one_body_tensor, - two_body=ham2_ActiveSpace.two_body_tensor) - return ham2_ActiveSpace + if cached_threshold != configured_threshold: + raise ValueError( + f"Active-space cache \"{ham2_filename}\" was built with coefficient threshold " + f"{cached_threshold}, but the configuration requests {configured_threshold}; " + "remove the cache and rerun." + ) + state.log(' '.join([ + f"Loaded \"{ham2_filename}\".", + "Continuing from after the active space is applied."])) + return ham2_ActiveSpace # ------------------------------------------------------------------------------------------------- @@ -494,7 +479,7 @@ def compute_metadata(state, ham3_Fermion2Qubit): # fermion-to-qubit transformation state.metadata["fermion-to-qubit operator mapping"] = ham3_Fermion2Qubit.f2q_mapping state.metadata["spin-orbital coefficient threshold (Hartrees)"] = ( - ham3_Fermion2Qubit.coefficient_threshold + state.config_hamiltonian.coefficient_threshold ) # number of terms in sum of Pauli strings state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms) diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py index 22b11025..62abf920 100644 --- a/hamiltonian_generator/tests/test_thresholding.py +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -92,7 +92,6 @@ def test_configuration_accepts_non_negative_real_thresholds(value): (float("nan"), ValueError), (float("inf"), ValueError), (-float("inf"), ValueError), - pytest.param(10**10000, ValueError, id="overflowing-real"), (True, TypeError), ("1e-8", TypeError), (None, TypeError), @@ -165,25 +164,18 @@ def test_active_space_matches_openfermion_and_respects_zero_threshold(): ) -def test_mapping_and_metadata_use_effective_operator_threshold(): +def test_mapping_and_metadata_record_configured_threshold(): state = _ActiveSpaceState(0.0) state.metadata = {} - active_operator = _interaction_operator( - 1.0, - threshold=DEFAULT_COEFFICIENT_THRESHOLD, - ) + active_operator = _interaction_operator(1.0, threshold=0.0) mapped_operator = hamgen.map_fermions_to_qubits(state, active_operator) hamgen.compute_metadata(state, mapped_operator) - assert mapped_operator.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD - assert mapped_operator.coefficient_threshold != state.config_hamiltonian.coefficient_threshold - assert state.metadata["spin-orbital coefficient threshold (Hartrees)"] == ( - DEFAULT_COEFFICIENT_THRESHOLD - ) + assert state.metadata["spin-orbital coefficient threshold (Hartrees)"] == 0.0 -def test_active_space_cache_tracks_threshold_and_recomputes_on_mismatch( +def test_active_space_cache_rejects_mismatch_and_rebuilds_after_removal( tmp_path, monkeypatch, ): @@ -209,10 +201,14 @@ def fake_apply_active_space(state, molecule): default_state = _CacheState(cache_path, DEFAULT_COEFFICIENT_THRESHOLD) loaded_legacy = hamgen.get_ham2(default_state) assert loaded_legacy.constant == legacy_operator.constant - assert loaded_legacy.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD assert calls == [] zero_state = _CacheState(cache_path, 0.0) + with pytest.raises(ValueError, match="remove the cache and rerun"): + hamgen.get_ham2(zero_state) + assert calls == [] + + cache_path.unlink() rebuilt = hamgen.get_ham2(zero_state) assert rebuilt is recomputed_operator assert calls == ["get_ham1", "apply_active_space"] @@ -237,49 +233,3 @@ def fake_apply_active_space(state, molecule): assert loaded_match.constant == recomputed_operator.constant assert loaded_match.coefficient_threshold == 0.0 assert calls == ["get_ham1", "apply_active_space"] - - -@pytest.mark.parametrize( - "invalid_threshold", - [ - False, - float("nan"), - np.array([0.0]), - pytest.param(10**10000, id="overflowing-real"), - ], -) -def test_active_space_cache_recomputes_invalid_threshold_provenance( - tmp_path, - monkeypatch, - invalid_threshold, -): - cache_path = tmp_path / "active-space.pickle" - cached_operator = _interaction_operator(1.0, threshold=invalid_threshold) - with open(cache_path, "wb") as file: - pickle.dump(cached_operator, file) - - recomputed_operator = _interaction_operator( - 2.0, - threshold=DEFAULT_COEFFICIENT_THRESHOLD, - ) - calls = [] - - def fake_get_ham1(state): - calls.append("get_ham1") - return object() - - def fake_apply_active_space(state, molecule): - calls.append("apply_active_space") - return recomputed_operator - - monkeypatch.setattr(hamgen, "get_ham1", fake_get_ham1) - monkeypatch.setattr(hamgen, "apply_active_space", fake_apply_active_space) - - state = _CacheState(cache_path, DEFAULT_COEFFICIENT_THRESHOLD) - rebuilt = hamgen.get_ham2(state) - - assert rebuilt is recomputed_operator - assert calls == ["get_ham1", "apply_active_space"] - with open(cache_path, "rb") as file: - persisted = pickle.load(file) - assert persisted.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD diff --git a/hamiltonian_generator/thresholding.py b/hamiltonian_generator/thresholding.py index 2092d54d..d3502df4 100644 --- a/hamiltonian_generator/thresholding.py +++ b/hamiltonian_generator/thresholding.py @@ -12,12 +12,7 @@ def validate_coefficient_threshold(value) -> float: if isinstance(value, bool) or not isinstance(value, numbers.Real): raise TypeError("coefficient threshold must be a real number") - try: - threshold = float(value) - except OverflowError as error: - raise ValueError( - "coefficient threshold must be finite and non-negative" - ) from error + threshold = float(value) if not math.isfinite(threshold) or threshold < 0.0: raise ValueError("coefficient threshold must be finite and non-negative") return threshold From 0f71af6c41ee82278b4a4a349c08422e5a0d4265 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 24 Jul 2026 10:14:05 -0600 Subject: [PATCH 07/16] Fix the new test_thresholding.py --- .../tests/test_thresholding.py | 258 +++++++++++++++--- 1 file changed, 224 insertions(+), 34 deletions(-) diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py index 62abf920..21fd5594 100644 --- a/hamiltonian_generator/tests/test_thresholding.py +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -1,11 +1,10 @@ """Regression tests for QHAT's configurable coefficient threshold.""" import pickle -from pathlib import Path import numpy as np import pytest -from openfermion import InteractionOperator +from openfermion import InteractionOperator, MolecularData from openfermion.chem.molecular_data import ( spinorb_from_spatial as openfermion_spinorb_from_spatial, ) @@ -18,9 +17,6 @@ ) -HAMGEN_DIRECTORY = Path(__file__).resolve().parents[1] - - class _ActiveSpaceState: def __init__(self, threshold=DEFAULT_COEFFICIENT_THRESHOLD): self.config_hamiltonian = HamiltonianConfiguration() @@ -57,6 +53,7 @@ def _interaction_operator(constant, threshold=None): np.diag([constant, 0.0]), np.zeros((2, 2, 2, 2)), ) + operator.hf_energy = -1.0 operator.asmeta = { "n_frz_occ_so": 0, @@ -68,20 +65,103 @@ def _interaction_operator(constant, threshold=None): operator.separation = 1.0 operator.hf_time = 0.0 operator.as_time = 0.0 + if threshold is not None: operator.coefficient_threshold = threshold + return operator +def _synthetic_active_space_molecule(): + """ + Build a deterministic MolecularData object without PySCF or pickle files. + + The synthetic molecule has: + + - 6 electrons + - 6 spatial orbitals + - 12 spin orbitals + - 1 frozen occupied spatial orbital + - 5 active spatial orbitals + + This reproduces the active-space ranges used by the original lithium + fixture without requiring a generated binary artifact. + """ + threshold = DEFAULT_COEFFIC_THRESHOLD + n_spatial_orbitals = 6 + + molecule = MolecularData( + geometry=[ + ("H", (0.0, 0.0, float(index))) + for index in range(n_spatial_orbitals) + ], + basis="sto-3g", + multiplicity=1, + charge=0, + ) + + molecule.n_orbitals = n_spatial_orbitals + molecule.n_qubits = 2 * n_spatial_orbitals + molecule.nuclear_repulsion = 1.0 + molecule.hf_energy = -1.0 + + molecule.basis = "sto-3g" + molecule.separation = None + molecule.hf_time = 0.0 + + one_body_integrals = np.zeros( + (n_spatial_orbitals, n_spatial_orbitals), + ) + + # All nonzero entries are inside the active spatial-orbital range [1, 6). + # + # Include coefficients below, exactly at, and above the default threshold. + one_body_integrals[1, 1] = threshold / 2.0 + one_body_integrals[1, 2] = threshold + one_body_integrals[2, 1] = threshold + one_body_integrals[2, 2] = -2.0 * threshold + + two_body_integrals = np.zeros( + ( + n_spatial_orbitals, + n_spatial_orbitals, + n_spatial_orbitals, + n_spatial_orbitals, + ), + ) + + two_body_integrals[1, 1, 1, 1] = threshold / 2.0 + two_body_integrals[1, 1, 1, 2] = threshold + two_body_integrals[2, 1, 1, 2] = -2.0 * threshold + + molecule.one_body_integrals = one_body_integrals + molecule.two_body_integrals = two_body_integrals + + return molecule + + def test_configuration_uses_openfermion_default_threshold(): config = HamiltonianConfiguration() - assert config.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD == 1.0e-8 + assert ( + config.coefficient_threshold + == DEFAULT_COEFFICIENT_THRESHOLD + == 1.0e-8 + ) -@pytest.mark.parametrize("value", [0, 1.0e-12, np.float64(1.0e-8)]) + +@pytest.mark.parametrize( + "value", + [ + 0, + 1.0e-12, + np.float64(1.0e-8), + ], +) def test_configuration_accepts_non_negative_real_thresholds(value): config = HamiltonianConfiguration() config.coefficient_threshold = value + assert config.coefficient_threshold == float(value) @@ -99,6 +179,7 @@ def test_configuration_accepts_non_negative_real_thresholds(value): ) def test_configuration_rejects_invalid_thresholds(value, exception): config = HamiltonianConfiguration() + with pytest.raises(exception): config.coefficient_threshold = value @@ -106,60 +187,104 @@ def test_configuration_rejects_invalid_thresholds(value, exception): def test_spinorb_default_parity_strict_boundary_and_zero_threshold(): threshold = DEFAULT_COEFFICIENT_THRESHOLD below_threshold = threshold / 2.0 + one_body = np.array( [ [below_threshold, threshold], [-2.0 * threshold, 0.0], ] ) + two_body = np.zeros((2, 2, 2, 2)) two_body[0, 0, 0, 0] = below_threshold two_body[0, 0, 0, 1] = threshold two_body[1, 0, 0, 1] = -2.0 * threshold - expected_one, expected_two = openfermion_spinorb_from_spatial(one_body, two_body) - default_one, default_two = spinorb_from_spatial(one_body, two_body) - zero_one, zero_two = spinorb_from_spatial(one_body, two_body, 0.0) + expected_one, expected_two = openfermion_spinorb_from_spatial( + one_body, + two_body, + ) + + default_one, default_two = spinorb_from_spatial( + one_body, + two_body, + ) + + zero_one, zero_two = spinorb_from_spatial( + one_body, + two_body, + 0.0, + ) + + np.testing.assert_array_equal( + default_one, + expected_one, + ) - np.testing.assert_array_equal(default_one, expected_one) - np.testing.assert_array_equal(default_two, expected_two) + np.testing.assert_array_equal( + default_two, + expected_two, + ) + # Coefficients strictly below 1e-8 are removed. assert default_one[0, 0] == 0.0 - assert default_one[0, 2] == threshold assert default_two[0, 1, 1, 0] == 0.0 + + # Coefficients exactly equal to 1e-8 are retained. + assert default_one[0, 2] == threshold assert default_two[0, 1, 1, 2] == threshold + # A zero threshold retains the smaller coefficients. assert zero_one[0, 0] == below_threshold assert zero_two[0, 1, 1, 0] == below_threshold def test_active_space_matches_openfermion_and_respects_zero_threshold(): - with open(HAMGEN_DIRECTORY / "diatomic_lithium.pickle", "rb") as file: - molecule = pickle.load(file) + molecule = _synthetic_active_space_molecule() - default_operator = hamgen.apply_active_space(_ActiveSpaceState(), molecule) expected_operator = molecule.get_molecular_hamiltonian( occupied_indices=range(1), active_indices=range(1, 6), ) + default_operator = hamgen.apply_active_space( + _ActiveSpaceState(), + molecule, + ) + assert default_operator.constant == expected_operator.constant + np.testing.assert_array_equal( default_operator.one_body_tensor, expected_operator.one_body_tensor, ) + np.testing.assert_array_equal( default_operator.two_body_tensor, expected_operator.two_body_tensor, ) - assert default_operator.coefficient_threshold == DEFAULT_COEFFICIENT_THRESHOLD - zero_operator = hamgen.apply_active_space(_ActiveSpaceState(0.0), molecule) + assert ( + default_operator.coefficient_threshold + == DEFAULT_COEFFICIENT_THRESHOLD + ) + + zero_operator = hamgen.apply_active_space( + _ActiveSpaceState(0.0), + molecule, + ) + assert zero_operator.coefficient_threshold == 0.0 - assert np.count_nonzero(zero_operator.one_body_tensor) > np.count_nonzero( + + assert np.count_nonzero( + zero_operator.one_body_tensor + ) > np.count_nonzero( default_operator.one_body_tensor ) - assert np.count_nonzero(zero_operator.two_body_tensor) > np.count_nonzero( + + assert np.count_nonzero( + zero_operator.two_body_tensor + ) > np.count_nonzero( default_operator.two_body_tensor ) @@ -167,12 +292,28 @@ def test_active_space_matches_openfermion_and_respects_zero_threshold(): def test_mapping_and_metadata_record_configured_threshold(): state = _ActiveSpaceState(0.0) state.metadata = {} - active_operator = _interaction_operator(1.0, threshold=0.0) - mapped_operator = hamgen.map_fermions_to_qubits(state, active_operator) - hamgen.compute_metadata(state, mapped_operator) + active_operator = _interaction_operator( + 1.0, + threshold=0.0, + ) + + mapped_operator = hamgen.map_fermions_to_qubits( + state, + active_operator, + ) + + hamgen.compute_metadata( + state, + mapped_operator, + ) - assert state.metadata["spin-orbital coefficient threshold (Hartrees)"] == 0.0 + assert ( + state.metadata[ + "spin-orbital coefficient threshold (Hartrees)" + ] + == 0.0 + ) def test_active_space_cache_rejects_mismatch_and_rebuilds_after_removal( @@ -180,12 +321,18 @@ def test_active_space_cache_rejects_mismatch_and_rebuilds_after_removal( monkeypatch, ): cache_path = tmp_path / "active-space.pickle" + legacy_operator = _interaction_operator(1.0) + with open(cache_path, "wb") as file: pickle.dump(legacy_operator, file) calls = [] - recomputed_operator = _interaction_operator(2.0, threshold=0.0) + + recomputed_operator = _interaction_operator( + 2.0, + threshold=0.0, + ) def fake_get_ham1(state): calls.append("get_ham1") @@ -195,41 +342,84 @@ def fake_apply_active_space(state, molecule): calls.append("apply_active_space") return recomputed_operator - monkeypatch.setattr(hamgen, "get_ham1", fake_get_ham1) - monkeypatch.setattr(hamgen, "apply_active_space", fake_apply_active_space) + monkeypatch.setattr( + hamgen, + "get_ham1", + fake_get_ham1, + ) + + monkeypatch.setattr( + hamgen, + "apply_active_space", + fake_apply_active_space, + ) + + default_state = _CacheState( + cache_path, + DEFAULT_COEFFICIENT_THRESHOLD, + ) - default_state = _CacheState(cache_path, DEFAULT_COEFFICIENT_THRESHOLD) loaded_legacy = hamgen.get_ham2(default_state) + assert loaded_legacy.constant == legacy_operator.constant assert calls == [] - zero_state = _CacheState(cache_path, 0.0) - with pytest.raises(ValueError, match="remove the cache and rerun"): + zero_state = _CacheState( + cache_path, + 0.0, + ) + + with pytest.raises( + ValueError, + match="remove the cache and rerun", + ): hamgen.get_ham2(zero_state) + assert calls == [] cache_path.unlink() + rebuilt = hamgen.get_ham2(zero_state) + assert rebuilt is recomputed_operator - assert calls == ["get_ham1", "apply_active_space"] + assert calls == [ + "get_ham1", + "apply_active_space", + ] with open(cache_path, "rb") as file: persisted = pickle.load(file) + assert persisted.coefficient_threshold == 0.0 assert persisted.constant == recomputed_operator.constant - with np.load(tmp_path / "active-space.tensors.npz") as tensors: - assert set(tensors.files) == {"constant", "one_body", "two_body"} + + tensor_path = tmp_path / "active-space.tensors.npz" + + with np.load(tensor_path) as tensors: + assert set(tensors.files) == { + "constant", + "one_body", + "two_body", + } + assert tensors["constant"] == recomputed_operator.constant + np.testing.assert_array_equal( tensors["one_body"], recomputed_operator.one_body_tensor, ) + np.testing.assert_array_equal( tensors["two_body"], recomputed_operator.two_body_tensor, ) loaded_match = hamgen.get_ham2(zero_state) + assert loaded_match.constant == recomputed_operator.constant assert loaded_match.coefficient_threshold == 0.0 - assert calls == ["get_ham1", "apply_active_space"] + + assert calls == [ + "get_ham1", + "apply_active_space", + ] \ No newline at end of file From b9c5b155aa53acb6d7a5ac0c4b607a88728d459a Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 24 Jul 2026 10:18:29 -0600 Subject: [PATCH 08/16] Fixed typos. --- hamiltonian_generator/tests/test_thresholding.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py index 21fd5594..896294f7 100644 --- a/hamiltonian_generator/tests/test_thresholding.py +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -87,7 +87,7 @@ def _synthetic_active_space_molecule(): This reproduces the active-space ranges used by the original lithium fixture without requiring a generated binary artifact. """ - threshold = DEFAULT_COEFFIC_THRESHOLD + threshold = DEFAULT_COEFFICIENT_THRESHOLD n_spatial_orbitals = 6 molecule = MolecularData( From cf7421d85ad26b9efb7ed59fd51537b2970cb944 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 29 Jul 2026 10:01:39 -0600 Subject: [PATCH 09/16] temp: created a tensor_ground_energy checker. --- .../compare_tensor_ground_energies.py | 298 ++++++++++++++++++ 1 file changed, 298 insertions(+) create mode 100644 hamiltonian_generator/compare_tensor_ground_energies.py diff --git a/hamiltonian_generator/compare_tensor_ground_energies.py b/hamiltonian_generator/compare_tensor_ground_energies.py new file mode 100644 index 00000000..69bc67ff --- /dev/null +++ b/hamiltonian_generator/compare_tensor_ground_energies.py @@ -0,0 +1,298 @@ +"""Compare ground-state energies from two QHAT tensor archives.""" + +from pathlib import Path + +import numpy as np +from openfermion import ( + InteractionOperator, + get_sparse_operator, + jw_get_ground_state_at_particle_number, +) + + +# --------------------------------------------------------------------------- +# Experiment settings +# --------------------------------------------------------------------------- + +MOLECULE_NAME = "H2" + +DEFAULT_THRESHOLD_FILE = Path( + "diatomic_hydrogen_as-002-002.tensors.npz" +) + +ZERO_THRESHOLD_FILE = Path( + "diatomic_hydrogen_0_as-002-002.tensors.npz" +) + +# Number of active electrons. +N_ELECTRONS = 2 + + +# --------------------------------------------------------------------------- +# Load QHAT tensors +# --------------------------------------------------------------------------- + +def load_interaction_operator(npz_path): + """Load a QHAT .tensors.npz file as an InteractionOperator.""" + npz_path = Path(npz_path) + + if not npz_path.is_file(): + raise FileNotFoundError( + f"Tensor file not found: {npz_path}" + ) + + with np.load(npz_path) as tensors: + required_arrays = { + "constant", + "one_body", + "two_body", + } + + available_arrays = set(tensors.files) + missing_arrays = required_arrays - available_arrays + + if missing_arrays: + raise KeyError( + f"{npz_path} is missing arrays: " + f"{sorted(missing_arrays)}" + ) + + constant = np.asarray( + tensors["constant"] + ).item() + + one_body = np.array( + tensors["one_body"], + copy=True, + ) + + two_body = np.array( + tensors["two_body"], + copy=True, + ) + + if one_body.ndim != 2: + raise ValueError( + f"one_body must be a matrix, but its shape is " + f"{one_body.shape}" + ) + + if one_body.shape[0] != one_body.shape[1]: + raise ValueError( + f"one_body must be square, but its shape is " + f"{one_body.shape}" + ) + + n_spin_orbitals = one_body.shape[0] + + expected_two_body_shape = ( + n_spin_orbitals, + n_spin_orbitals, + n_spin_orbitals, + n_spin_orbitals, + ) + + if two_body.shape != expected_two_body_shape: + raise ValueError( + f"two_body has shape {two_body.shape}; expected " + f"{expected_two_body_shape}" + ) + + if N_ELECTRONS > n_spin_orbitals: + raise ValueError( + f"N_ELECTRONS={N_ELECTRONS} exceeds the number of " + f"spin orbitals, {n_spin_orbitals}" + ) + + # The saved two_body array is already the final + # InteractionOperator.two_body_tensor. + # + # Do not multiply it by 0.5 again. + interaction_operator = InteractionOperator( + constant=constant, + one_body_tensor=one_body, + two_body_tensor=two_body, + ) + + return interaction_operator, n_spin_orbitals + + +# --------------------------------------------------------------------------- +# Ground-state calculation +# --------------------------------------------------------------------------- + +def calculate_ground_state( + interaction_operator, + n_spin_orbitals, + n_electrons, +): + """Calculate the ground state in a fixed-electron-number sector.""" + sparse_hamiltonian = get_sparse_operator( + interaction_operator, + n_qubits=n_spin_orbitals, + ) + + energy, state = jw_get_ground_state_at_particle_number( + sparse_hamiltonian, + n_electrons, + ) + + energy = float( + np.real_if_close(energy) + ) + + return energy, state + + +# --------------------------------------------------------------------------- +# Main comparison +# --------------------------------------------------------------------------- + +def main(): + default_operator, default_n_spin_orbitals = ( + load_interaction_operator( + DEFAULT_THRESHOLD_FILE + ) + ) + + zero_operator, zero_n_spin_orbitals = ( + load_interaction_operator( + ZERO_THRESHOLD_FILE + ) + ) + + if default_n_spin_orbitals != zero_n_spin_orbitals: + raise ValueError( + "The tensor files have different active-space sizes: " + f"{default_n_spin_orbitals} and " + f"{zero_n_spin_orbitals} spin orbitals" + ) + + default_energy, default_state = calculate_ground_state( + default_operator, + default_n_spin_orbitals, + N_ELECTRONS, + ) + + zero_energy, zero_state = calculate_ground_state( + zero_operator, + zero_n_spin_orbitals, + N_ELECTRONS, + ) + + signed_energy_difference = ( + zero_energy - default_energy + ) + + absolute_energy_difference = abs( + signed_energy_difference + ) + + state_overlap = abs( + np.vdot( + default_state, + zero_state, + ) + ) + + state_fidelity = state_overlap**2 + + one_body_difference = np.linalg.norm( + zero_operator.one_body_tensor + - default_operator.one_body_tensor + ) + + two_body_difference = np.linalg.norm( + zero_operator.two_body_tensor + - default_operator.two_body_tensor + ) + + print("=" * 72) + print( + f"{MOLECULE_NAME} active-space ground-state " + "energy comparison" + ) + print("=" * 72) + + print( + f"Number of spin orbitals : " + f"{default_n_spin_orbitals}" + ) + + print( + f"Number of electrons : " + f"{N_ELECTRONS}" + ) + + print() + + print("Default threshold tensor file:") + print(f" {DEFAULT_THRESHOLD_FILE}") + print( + f" Ground-state energy = " + f"{default_energy:.16f} Hartree" + ) + + print() + + print("Zero-threshold tensor file:") + print(f" {ZERO_THRESHOLD_FILE}") + print( + f" Ground-state energy = " + f"{zero_energy:.16f} Hartree" + ) + + print() + + print("Energy comparison:") + print( + " E(zero threshold) - E(default threshold)" + f" = {signed_energy_difference:+.16e} Hartree" + ) + + print( + f" Absolute difference = " + f"{absolute_energy_difference:.16e} Hartree" + ) + + print( + f" Absolute difference = " + f"{absolute_energy_difference * 1.0e3:.16e} " + "mHartree" + ) + + print( + f" Absolute difference = " + f"{absolute_energy_difference * 1.0e6:.16e} " + "microHartree" + ) + + print() + + print("Ground-state comparison:") + print( + f" State overlap = " + f"{state_overlap:.16f}" + ) + + print( + f" State fidelity = " + f"{state_fidelity:.16f}" + ) + + print() + + print("Tensor differences:") + print( + f" ||delta one_body||_F = " + f"{one_body_difference:.16e}" + ) + + print( + f" ||delta two_body||_F = " + f"{two_body_difference:.16e}" + ) + + +if __name__ == "__main__": + main() \ No newline at end of file From 445a7312ffd08720c9c9f27be43929fede237603 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 29 Jul 2026 10:55:47 -0600 Subject: [PATCH 10/16] Introduce inline Pauli String Coefficient Threshold in hamgen.py --- hamiltonian_generator/hamgen.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 9d17c1d3..aad90ee5 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -27,6 +27,11 @@ logger = logging.getLogger(__name__) +# Fixed baseline threshold for Pauli-string coefficients. +# This is intentionally not user-configurable in this PR. +# TODO: Consider making the Pauli-string coefficient threshold configurable in a future update. +PAULI_STRING_COEFFICIENT_THRESHOLD = 1.0e-8 + # ------------------------------------------------------------------------------------------------- def load_configuration() -> State: From 1eddfee510bc5c62965173664c7586b169f0eed7 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 29 Jul 2026 11:08:53 -0600 Subject: [PATCH 11/16] Modified the map_fermion_to_qubits to use inline threshold. --- hamiltonian_generator/hamgen.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index aad90ee5..b1c14725 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -285,6 +285,7 @@ def apply_active_space(state, ham1_HartreeFock): def map_fermions_to_qubits(state, ham2_ActiveSpace): # Transform fermionic operators to qubit operators + # A fixed coefficient threshold is applied to the final Pauli strings produced by both the Jordan-Wigner and Bravyi-Kitaev mappings. tstart = time.time() ham3_Fermion2Qubit = None mapping_name = state.config_hamiltonian.fermion_to_qubit_name() @@ -295,6 +296,27 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): else: mapping = state.config_hamiltonian.f2q_mapping() raise NotImplementedError(f"invalid fermion-to-qubit mapping \"{mapping}\"") + + # Apply the same fixed Pauli-string threshold to JW and BK, removing only coefficients with magnitude strictly below the cutoff. + number_of_terms_before = len(ham3_Fermion2Qubit.terms) + + ham3_Fermion2Qubit.terms = { + pauli_string: coefficient + for pauli_string, coefficient + in ham3_Fermion2Qubit.terms.items() + if abs(coefficient) >= PAULI_STRING_COEFFICIENT_THRESHOLD + } + + number_of_terms_removed = ( + number_of_terms_before - len(ham3_Fermion2Qubit.terms) + ) + + state.log( + "Applied fixed Pauli-string coefficient threshold " + f"{PAULI_STRING_COEFFICIENT_THRESHOLD:.1e} Hartrees; " + f"removed {number_of_terms_removed} terms." + ) + tstop = time.time() # Propagate previously-computed metadata that we need to preserve across all run modes ham3_Fermion2Qubit.hf_energy = ham2_ActiveSpace.hf_energy @@ -305,6 +327,12 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): ham3_Fermion2Qubit.hf_time = ham2_ActiveSpace.hf_time ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time ham3_Fermion2Qubit.f2q_time = tstop - tstart + + # Record the fixed Pauli threshold on the generated operator. + ham3_Fermion2Qubit.pauli_coefficient_threshold = ( + PAULI_STRING_COEFFICIENT_THRESHOLD + ) + # Return result return ham3_Fermion2Qubit From c0fd3ee0bab86cc53e5eceab7a08940349f9d204 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 29 Jul 2026 11:38:47 -0600 Subject: [PATCH 12/16] Added the meta data printing for the pauli coefficient threshold --- hamiltonian_generator/hamgen.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index b1c14725..9ba17c2e 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -514,6 +514,10 @@ def compute_metadata(state, ham3_Fermion2Qubit): state.metadata["spin-orbital coefficient threshold (Hartrees)"] = ( state.config_hamiltonian.coefficient_threshold ) + # Fixed threshold applied after the JW or BK mapping + state.metadata["Pauli-string coefficient threshold (Hartrees)"] = ( + ham3_Fermion2Qubit.pauli_coefficient_threshold + ) # number of terms in sum of Pauli strings state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms) # one-norm of sum of Pauli strings From c49feeebe0a0c732a2e6a4cb30d0be73164d85fa Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Wed, 29 Jul 2026 12:11:28 -0600 Subject: [PATCH 13/16] test: add JW and BK Pauli-string threshold regression coverage --- .../tests/test_thresholding.py | 103 +++++++++++++++++- 1 file changed, 101 insertions(+), 2 deletions(-) diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py index 896294f7..52619f57 100644 --- a/hamiltonian_generator/tests/test_thresholding.py +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -1,10 +1,10 @@ -"""Regression tests for QHAT's configurable coefficient threshold.""" +"""Regression tests for QHAT coefficient thresholds.""" import pickle import numpy as np import pytest -from openfermion import InteractionOperator, MolecularData +from openfermion import InteractionOperator, MolecularData, QubitOperator from openfermion.chem.molecular_data import ( spinorb_from_spatial as openfermion_spinorb_from_spatial, ) @@ -315,6 +315,105 @@ def test_mapping_and_metadata_record_configured_threshold(): == 0.0 ) + assert ( + state.metadata[ + "Pauli-string coefficient threshold (Hartrees)" + ] + == hamgen.PAULI_STRING_COEFFICIENT_THRESHOLD + ) + + +@pytest.mark.parametrize( + ("mapping_setting", "mapping_name", "mapper_name"), + [ + ( + "Jordan-Wigner", + "jordan-wigner", + "jordan_wigner", + ), + ( + "Bravyi-Kitaev", + "bravyi-kitaev", + "bravyi_kitaev", + ), + ], +) +def test_pauli_string_threshold_applies_to_jw_and_bk( + monkeypatch, + mapping_setting, + mapping_name, + mapper_name, +): + threshold = hamgen.PAULI_STRING_COEFFICIENT_THRESHOLD + + # Assign terms directly so OpenFermion does not remove the synthetic + # below-threshold coefficient before QHAT receives the mapped operator. + mapped_operator = QubitOperator() + mapped_operator.terms = { + tuple(): 1.0, + ((0, "X"),): 0.5 * threshold, + ((0, "Y"),): threshold, + ((0, "Z"),): 2.0 * threshold, + } + + def fake_mapping(active_operator): + return mapped_operator + + monkeypatch.setattr( + hamgen, + mapper_name, + fake_mapping, + ) + + state = _ActiveSpaceState() + state.metadata = {} + state.config_hamiltonian.f2q_mapping = mapping_setting + + active_operator = _interaction_operator( + 1.0, + threshold=DEFAULT_COEFFICIENT_THRESHOLD, + ) + + result = hamgen.map_fermions_to_qubits( + state, + active_operator, + ) + + # A coefficient strictly below the fixed threshold is removed. + assert ((0, "X"),) not in result.terms + + # A coefficient exactly equal to the threshold is retained. + assert result.terms[((0, "Y"),)] == threshold + + # A coefficient above the threshold is retained unchanged. + assert result.terms[((0, "Z"),)] == 2.0 * threshold + + # The identity term is also retained unchanged. + assert result.terms[tuple()] == 1.0 + + assert result.f2q_mapping == mapping_name + assert result.pauli_coefficient_threshold == threshold + assert len(result.terms) == 3 + + hamgen.compute_metadata( + state, + result, + ) + + assert ( + state.metadata[ + "Pauli-string coefficient threshold (Hartrees)" + ] + == threshold + ) + + assert ( + state.metadata[ + "number of terms in sum of Pauli strings" + ] + == 3 + ) + def test_active_space_cache_rejects_mismatch_and_rebuilds_after_removal( tmp_path, From 8f57bf0d83e2f85e8a57c1b0820ec0103438e2b6 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Thu, 30 Jul 2026 01:09:26 -0600 Subject: [PATCH 14/16] Address coefficient threshold review feedback --- hamiltonian_generator/README.md | 8 +- .../compare_tensor_ground_energies.py | 298 ------------------ 2 files changed, 6 insertions(+), 300 deletions(-) delete mode 100644 hamiltonian_generator/compare_tensor_ground_energies.py diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index 4994f3d0..6ff9d77d 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -83,8 +83,12 @@ retained. Its default is `1e-8`; set it to `0.0` to disable filtering at this s hamiltonian.coefficient_threshold = 0.0 ``` -The value must be finite and non-negative. This setting affects only tensor construction; it does -not change OpenFermion's standard JW/BK mapping, which can still discard small coefficients. +The value must be finite and non-negative. This setting affects only tensor construction. After +OpenFermion's standard JW/BK mapping, QHAT also removes final Pauli strings whose coefficient +magnitude is strictly less than a fixed `1e-8` Hartree cutoff. This second cutoff is currently not +configurable, so setting `hamiltonian.coefficient_threshold = 0.0` preserves small spin-orbital +tensor entries but does not disable filtering of the final Pauli Hamiltonian. Both thresholds are +recorded separately in the generated `.dat` file metadata. Active-space pickles record the threshold used to construct their tensors. If the configured value does not match `[filestub]_[astag].pickle`, QHAT stops rather than silently reusing incompatible tensors. Remove that active-space pickle and rerun; QHAT will rebuild it and overwrite diff --git a/hamiltonian_generator/compare_tensor_ground_energies.py b/hamiltonian_generator/compare_tensor_ground_energies.py deleted file mode 100644 index 69bc67ff..00000000 --- a/hamiltonian_generator/compare_tensor_ground_energies.py +++ /dev/null @@ -1,298 +0,0 @@ -"""Compare ground-state energies from two QHAT tensor archives.""" - -from pathlib import Path - -import numpy as np -from openfermion import ( - InteractionOperator, - get_sparse_operator, - jw_get_ground_state_at_particle_number, -) - - -# --------------------------------------------------------------------------- -# Experiment settings -# --------------------------------------------------------------------------- - -MOLECULE_NAME = "H2" - -DEFAULT_THRESHOLD_FILE = Path( - "diatomic_hydrogen_as-002-002.tensors.npz" -) - -ZERO_THRESHOLD_FILE = Path( - "diatomic_hydrogen_0_as-002-002.tensors.npz" -) - -# Number of active electrons. -N_ELECTRONS = 2 - - -# --------------------------------------------------------------------------- -# Load QHAT tensors -# --------------------------------------------------------------------------- - -def load_interaction_operator(npz_path): - """Load a QHAT .tensors.npz file as an InteractionOperator.""" - npz_path = Path(npz_path) - - if not npz_path.is_file(): - raise FileNotFoundError( - f"Tensor file not found: {npz_path}" - ) - - with np.load(npz_path) as tensors: - required_arrays = { - "constant", - "one_body", - "two_body", - } - - available_arrays = set(tensors.files) - missing_arrays = required_arrays - available_arrays - - if missing_arrays: - raise KeyError( - f"{npz_path} is missing arrays: " - f"{sorted(missing_arrays)}" - ) - - constant = np.asarray( - tensors["constant"] - ).item() - - one_body = np.array( - tensors["one_body"], - copy=True, - ) - - two_body = np.array( - tensors["two_body"], - copy=True, - ) - - if one_body.ndim != 2: - raise ValueError( - f"one_body must be a matrix, but its shape is " - f"{one_body.shape}" - ) - - if one_body.shape[0] != one_body.shape[1]: - raise ValueError( - f"one_body must be square, but its shape is " - f"{one_body.shape}" - ) - - n_spin_orbitals = one_body.shape[0] - - expected_two_body_shape = ( - n_spin_orbitals, - n_spin_orbitals, - n_spin_orbitals, - n_spin_orbitals, - ) - - if two_body.shape != expected_two_body_shape: - raise ValueError( - f"two_body has shape {two_body.shape}; expected " - f"{expected_two_body_shape}" - ) - - if N_ELECTRONS > n_spin_orbitals: - raise ValueError( - f"N_ELECTRONS={N_ELECTRONS} exceeds the number of " - f"spin orbitals, {n_spin_orbitals}" - ) - - # The saved two_body array is already the final - # InteractionOperator.two_body_tensor. - # - # Do not multiply it by 0.5 again. - interaction_operator = InteractionOperator( - constant=constant, - one_body_tensor=one_body, - two_body_tensor=two_body, - ) - - return interaction_operator, n_spin_orbitals - - -# --------------------------------------------------------------------------- -# Ground-state calculation -# --------------------------------------------------------------------------- - -def calculate_ground_state( - interaction_operator, - n_spin_orbitals, - n_electrons, -): - """Calculate the ground state in a fixed-electron-number sector.""" - sparse_hamiltonian = get_sparse_operator( - interaction_operator, - n_qubits=n_spin_orbitals, - ) - - energy, state = jw_get_ground_state_at_particle_number( - sparse_hamiltonian, - n_electrons, - ) - - energy = float( - np.real_if_close(energy) - ) - - return energy, state - - -# --------------------------------------------------------------------------- -# Main comparison -# --------------------------------------------------------------------------- - -def main(): - default_operator, default_n_spin_orbitals = ( - load_interaction_operator( - DEFAULT_THRESHOLD_FILE - ) - ) - - zero_operator, zero_n_spin_orbitals = ( - load_interaction_operator( - ZERO_THRESHOLD_FILE - ) - ) - - if default_n_spin_orbitals != zero_n_spin_orbitals: - raise ValueError( - "The tensor files have different active-space sizes: " - f"{default_n_spin_orbitals} and " - f"{zero_n_spin_orbitals} spin orbitals" - ) - - default_energy, default_state = calculate_ground_state( - default_operator, - default_n_spin_orbitals, - N_ELECTRONS, - ) - - zero_energy, zero_state = calculate_ground_state( - zero_operator, - zero_n_spin_orbitals, - N_ELECTRONS, - ) - - signed_energy_difference = ( - zero_energy - default_energy - ) - - absolute_energy_difference = abs( - signed_energy_difference - ) - - state_overlap = abs( - np.vdot( - default_state, - zero_state, - ) - ) - - state_fidelity = state_overlap**2 - - one_body_difference = np.linalg.norm( - zero_operator.one_body_tensor - - default_operator.one_body_tensor - ) - - two_body_difference = np.linalg.norm( - zero_operator.two_body_tensor - - default_operator.two_body_tensor - ) - - print("=" * 72) - print( - f"{MOLECULE_NAME} active-space ground-state " - "energy comparison" - ) - print("=" * 72) - - print( - f"Number of spin orbitals : " - f"{default_n_spin_orbitals}" - ) - - print( - f"Number of electrons : " - f"{N_ELECTRONS}" - ) - - print() - - print("Default threshold tensor file:") - print(f" {DEFAULT_THRESHOLD_FILE}") - print( - f" Ground-state energy = " - f"{default_energy:.16f} Hartree" - ) - - print() - - print("Zero-threshold tensor file:") - print(f" {ZERO_THRESHOLD_FILE}") - print( - f" Ground-state energy = " - f"{zero_energy:.16f} Hartree" - ) - - print() - - print("Energy comparison:") - print( - " E(zero threshold) - E(default threshold)" - f" = {signed_energy_difference:+.16e} Hartree" - ) - - print( - f" Absolute difference = " - f"{absolute_energy_difference:.16e} Hartree" - ) - - print( - f" Absolute difference = " - f"{absolute_energy_difference * 1.0e3:.16e} " - "mHartree" - ) - - print( - f" Absolute difference = " - f"{absolute_energy_difference * 1.0e6:.16e} " - "microHartree" - ) - - print() - - print("Ground-state comparison:") - print( - f" State overlap = " - f"{state_overlap:.16f}" - ) - - print( - f" State fidelity = " - f"{state_fidelity:.16f}" - ) - - print() - - print("Tensor differences:") - print( - f" ||delta one_body||_F = " - f"{one_body_difference:.16e}" - ) - - print( - f" ||delta two_body||_F = " - f"{two_body_difference:.16e}" - ) - - -if __name__ == "__main__": - main() \ No newline at end of file From c43de925d92e88040fa4d58ffc93a8b14dfeeddb Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Thu, 13 Aug 2026 07:19:13 -0600 Subject: [PATCH 15/16] Removal of Pauli Threshold. --- hamiltonian_generator/README.md | 10 +++--- hamiltonian_generator/hamgen.py | 35 ------------------- .../tests/test_thresholding.py | 31 ++++++---------- 3 files changed, 15 insertions(+), 61 deletions(-) diff --git a/hamiltonian_generator/README.md b/hamiltonian_generator/README.md index 6ff9d77d..588028e3 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -83,12 +83,10 @@ retained. Its default is `1e-8`; set it to `0.0` to disable filtering at this s hamiltonian.coefficient_threshold = 0.0 ``` -The value must be finite and non-negative. This setting affects only tensor construction. After -OpenFermion's standard JW/BK mapping, QHAT also removes final Pauli strings whose coefficient -magnitude is strictly less than a fixed `1e-8` Hartree cutoff. This second cutoff is currently not -configurable, so setting `hamiltonian.coefficient_threshold = 0.0` preserves small spin-orbital -tensor entries but does not disable filtering of the final Pauli Hamiltonian. Both thresholds are -recorded separately in the generated `.dat` file metadata. +The value must be finite and non-negative. This setting affects only tensor construction; QHAT +does not apply an additional coefficient cutoff after OpenFermion's standard JW/BK mapping. The +tensor threshold is recorded in the generated `.dat` file metadata. + Active-space pickles record the threshold used to construct their tensors. If the configured value does not match `[filestub]_[astag].pickle`, QHAT stops rather than silently reusing incompatible tensors. Remove that active-space pickle and rerun; QHAT will rebuild it and overwrite diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 9ba17c2e..4df29f47 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -27,11 +27,6 @@ logger = logging.getLogger(__name__) -# Fixed baseline threshold for Pauli-string coefficients. -# This is intentionally not user-configurable in this PR. -# TODO: Consider making the Pauli-string coefficient threshold configurable in a future update. -PAULI_STRING_COEFFICIENT_THRESHOLD = 1.0e-8 - # ------------------------------------------------------------------------------------------------- def load_configuration() -> State: @@ -285,7 +280,6 @@ def apply_active_space(state, ham1_HartreeFock): def map_fermions_to_qubits(state, ham2_ActiveSpace): # Transform fermionic operators to qubit operators - # A fixed coefficient threshold is applied to the final Pauli strings produced by both the Jordan-Wigner and Bravyi-Kitaev mappings. tstart = time.time() ham3_Fermion2Qubit = None mapping_name = state.config_hamiltonian.fermion_to_qubit_name() @@ -297,26 +291,6 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): mapping = state.config_hamiltonian.f2q_mapping() raise NotImplementedError(f"invalid fermion-to-qubit mapping \"{mapping}\"") - # Apply the same fixed Pauli-string threshold to JW and BK, removing only coefficients with magnitude strictly below the cutoff. - number_of_terms_before = len(ham3_Fermion2Qubit.terms) - - ham3_Fermion2Qubit.terms = { - pauli_string: coefficient - for pauli_string, coefficient - in ham3_Fermion2Qubit.terms.items() - if abs(coefficient) >= PAULI_STRING_COEFFICIENT_THRESHOLD - } - - number_of_terms_removed = ( - number_of_terms_before - len(ham3_Fermion2Qubit.terms) - ) - - state.log( - "Applied fixed Pauli-string coefficient threshold " - f"{PAULI_STRING_COEFFICIENT_THRESHOLD:.1e} Hartrees; " - f"removed {number_of_terms_removed} terms." - ) - tstop = time.time() # Propagate previously-computed metadata that we need to preserve across all run modes ham3_Fermion2Qubit.hf_energy = ham2_ActiveSpace.hf_energy @@ -328,11 +302,6 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace): ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time ham3_Fermion2Qubit.f2q_time = tstop - tstart - # Record the fixed Pauli threshold on the generated operator. - ham3_Fermion2Qubit.pauli_coefficient_threshold = ( - PAULI_STRING_COEFFICIENT_THRESHOLD - ) - # Return result return ham3_Fermion2Qubit @@ -514,10 +483,6 @@ def compute_metadata(state, ham3_Fermion2Qubit): state.metadata["spin-orbital coefficient threshold (Hartrees)"] = ( state.config_hamiltonian.coefficient_threshold ) - # Fixed threshold applied after the JW or BK mapping - state.metadata["Pauli-string coefficient threshold (Hartrees)"] = ( - ham3_Fermion2Qubit.pauli_coefficient_threshold - ) # number of terms in sum of Pauli strings state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms) # one-norm of sum of Pauli strings diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py index 52619f57..ecb92ab0 100644 --- a/hamiltonian_generator/tests/test_thresholding.py +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -315,12 +315,7 @@ def test_mapping_and_metadata_record_configured_threshold(): == 0.0 ) - assert ( - state.metadata[ - "Pauli-string coefficient threshold (Hartrees)" - ] - == hamgen.PAULI_STRING_COEFFICIENT_THRESHOLD - ) + assert "Pauli-string coefficient threshold (Hartrees)" not in state.metadata @pytest.mark.parametrize( @@ -338,13 +333,13 @@ def test_mapping_and_metadata_record_configured_threshold(): ), ], ) -def test_pauli_string_threshold_applies_to_jw_and_bk( +def test_mapping_does_not_apply_additional_pauli_string_threshold( monkeypatch, mapping_setting, mapping_name, mapper_name, ): - threshold = hamgen.PAULI_STRING_COEFFICIENT_THRESHOLD + threshold = DEFAULT_COEFFICIENT_THRESHOLD # Assign terms directly so OpenFermion does not remove the synthetic # below-threshold coefficient before QHAT receives the mapped operator. @@ -379,8 +374,9 @@ def fake_mapping(active_operator): active_operator, ) - # A coefficient strictly below the fixed threshold is removed. - assert ((0, "X"),) not in result.terms + # QHAT preserves all terms returned by OpenFermion, including coefficients + # below the tensor-construction threshold. + assert result.terms[((0, "X"),)] == 0.5 * threshold # A coefficient exactly equal to the threshold is retained. assert result.terms[((0, "Y"),)] == threshold @@ -392,26 +388,21 @@ def fake_mapping(active_operator): assert result.terms[tuple()] == 1.0 assert result.f2q_mapping == mapping_name - assert result.pauli_coefficient_threshold == threshold - assert len(result.terms) == 3 + assert not hasattr(result, "pauli_coefficient_threshold") + assert len(result.terms) == 4 hamgen.compute_metadata( state, result, ) - assert ( - state.metadata[ - "Pauli-string coefficient threshold (Hartrees)" - ] - == threshold - ) + assert "Pauli-string coefficient threshold (Hartrees)" not in state.metadata assert ( state.metadata[ "number of terms in sum of Pauli strings" ] - == 3 + == 4 ) @@ -521,4 +512,4 @@ def fake_apply_active_space(state, molecule): assert calls == [ "get_ham1", "apply_active_space", - ] \ No newline at end of file + ] From 851d9ed13e5e848d1d3bbe40d107d1aa572bfe9b Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Tue, 25 Aug 2026 05:14:48 -0600 Subject: [PATCH 16/16] Fix Threshold merge regressions --- hamiltonian_generator/config.py | 4 ---- hamiltonian_generator/hamgen_types.py | 5 +---- hamiltonian_generator/tests/test_thresholding.py | 9 +++++++++ 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/hamiltonian_generator/config.py b/hamiltonian_generator/config.py index 3fe04e64..b3ad9a5e 100644 --- a/hamiltonian_generator/config.py +++ b/hamiltonian_generator/config.py @@ -24,12 +24,8 @@ hamiltonian.basis = params.get("basis", "sto-3g") # Select the atomic basis functions -<<<<<<< HEAD # Absolute cutoff applied while building spin-orbital tensors; the default is 1e-8. hamiltonian.coefficient_threshold = 1e-8 -hamiltonian.num_active_occupied = 4 # Specify the active space -======= hamiltonian.num_active_occupied = 4 # Specify the active space ->>>>>>> upstream/main hamiltonian.num_active_vacant = 6 diff --git a/hamiltonian_generator/hamgen_types.py b/hamiltonian_generator/hamgen_types.py index 4944c17e..559bea1d 100644 --- a/hamiltonian_generator/hamgen_types.py +++ b/hamiltonian_generator/hamgen_types.py @@ -6,14 +6,11 @@ import subprocess import sys -<<<<<<< HEAD +from qhat.common.git_utils import get_git_hash from qhat.hamiltonian_generator.thresholding import ( DEFAULT_COEFFICIENT_THRESHOLD, validate_coefficient_threshold, ) -======= -from qhat.common.git_utils import get_git_hash ->>>>>>> upstream/main logger = logging.getLogger(__name__) diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py index ecb92ab0..4049ec75 100644 --- a/hamiltonian_generator/tests/test_thresholding.py +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -31,9 +31,18 @@ def log_verbose(self, *args, **kwargs): pass +class _PassthroughPaths: + def get_cache_path(self, filename): + return filename + + def get_output_path(self, filename): + return filename + + class _CacheState: def __init__(self, cache_path, threshold): self.cache_path = cache_path + self.config_general = _PassthroughPaths() self.config_hamiltonian = HamiltonianConfiguration() self.config_hamiltonian.coefficient_threshold = threshold