diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 98d14ec4..aadb93ad 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 b56645e9..588028e3 100644 --- a/hamiltonian_generator/README.md +++ b/hamiltonian_generator/README.md @@ -74,6 +74,25 @@ 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, 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 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 +`[filestub]_[astag].tensors.npz`. Pickles created before this setting existed are interpreted as +using the historical `1e-8` cutoff. + ## Generated Files This script will generate various files for intermediate and final stages of the process. These diff --git a/hamiltonian_generator/config.py b/hamiltonian_generator/config.py index 0f6ff3d1..b3ad9a5e 100644 --- a/hamiltonian_generator/config.py +++ b/hamiltonian_generator/config.py @@ -24,6 +24,8 @@ hamiltonian.basis = params.get("basis", "sto-3g") # Select the atomic basis functions +# 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_vacant = 6 - diff --git a/hamiltonian_generator/hamgen.py b/hamiltonian_generator/hamgen.py index 3b1e4214..b488494f 100644 --- a/hamiltonian_generator/hamgen.py +++ b/hamiltonian_generator/hamgen.py @@ -33,6 +33,10 @@ HamiltonianConfiguration, State, ) +from qhat.hamiltonian_generator.thresholding import ( + DEFAULT_COEFFICIENT_THRESHOLD, + spinorb_from_spatial, +) logger = logging.getLogger(__name__) @@ -287,9 +291,21 @@ 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( - occupied_indices=range(idx_act_occ), - active_indices=range(idx_act_occ, idx_frz_vac), + 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, ) tstop = time.time() # Pass along metadata that needs to be preserved @@ -299,6 +315,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 # ------------------------------------------------------------------------------------------------- @@ -315,6 +332,7 @@ def map_fermions_to_qubits(state, 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 @@ -325,6 +343,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 + # Return result return ham3_Fermion2Qubit @@ -437,6 +456,18 @@ def get_ham2(state): two_body=ham2_ActiveSpace.two_body_tensor) return ham2_ActiveSpace else: + cached_threshold = getattr( + ham2_ActiveSpace, + "coefficient_threshold", + DEFAULT_COEFFICIENT_THRESHOLD, + ) + configured_threshold = state.config_hamiltonian.coefficient_threshold + 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_cache_path}\".", "Continuing from after the active space is applied."])) @@ -498,6 +529,9 @@ 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["spin-orbital coefficient threshold (Hartrees)"] = ( + 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) # one-norm of sum of Pauli strings diff --git a/hamiltonian_generator/hamgen_types.py b/hamiltonian_generator/hamgen_types.py index a8b2cdd5..559bea1d 100644 --- a/hamiltonian_generator/hamgen_types.py +++ b/hamiltonian_generator/hamgen_types.py @@ -7,6 +7,10 @@ import sys from qhat.common.git_utils import get_git_hash +from qhat.hamiltonian_generator.thresholding import ( + DEFAULT_COEFFICIENT_THRESHOLD, + validate_coefficient_threshold, +) logger = logging.getLogger(__name__) @@ -40,9 +44,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 @@ -208,4 +219,3 @@ def filename_ham3(self): ext=self.config_general.ham3_ext()) # ------------------------------------------------------------------------------------------------- - diff --git a/hamiltonian_generator/tests/test_thresholding.py b/hamiltonian_generator/tests/test_thresholding.py new file mode 100644 index 00000000..4049ec75 --- /dev/null +++ b/hamiltonian_generator/tests/test_thresholding.py @@ -0,0 +1,524 @@ +"""Regression tests for QHAT coefficient thresholds.""" + +import pickle + +import numpy as np +import pytest +from openfermion import InteractionOperator, MolecularData, QubitOperator +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, +) + + +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 _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 + + 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 _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_COEFFICIENT_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 + ) + + +@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), + (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, + ) + + # Coefficients strictly below 1e-8 are removed. + assert default_one[0, 0] == 0.0 + 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(): + molecule = _synthetic_active_space_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 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_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, + ) + + assert ( + state.metadata[ + "spin-orbital coefficient threshold (Hartrees)" + ] + == 0.0 + ) + + assert "Pauli-string coefficient threshold (Hartrees)" not in state.metadata + + +@pytest.mark.parametrize( + ("mapping_setting", "mapping_name", "mapper_name"), + [ + ( + "Jordan-Wigner", + "jordan-wigner", + "jordan_wigner", + ), + ( + "Bravyi-Kitaev", + "bravyi-kitaev", + "bravyi_kitaev", + ), + ], +) +def test_mapping_does_not_apply_additional_pauli_string_threshold( + monkeypatch, + mapping_setting, + mapping_name, + mapper_name, +): + threshold = DEFAULT_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, + ) + + # 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 + + # 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 not hasattr(result, "pauli_coefficient_threshold") + assert len(result.terms) == 4 + + hamgen.compute_metadata( + state, + result, + ) + + assert "Pauli-string coefficient threshold (Hartrees)" not in state.metadata + + assert ( + state.metadata[ + "number of terms in sum of Pauli strings" + ] + == 4 + ) + + +def test_active_space_cache_rejects_mismatch_and_rebuilds_after_removal( + 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 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", + ] + + with open(cache_path, "rb") as file: + persisted = pickle.load(file) + + assert persisted.coefficient_threshold == 0.0 + assert persisted.constant == recomputed_operator.constant + + 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", + ] diff --git a/hamiltonian_generator/thresholding.py b/hamiltonian_generator/thresholding.py new file mode 100644 index 00000000..d3502df4 --- /dev/null +++ b/hamiltonian_generator/thresholding.py @@ -0,0 +1,51 @@ +"""QHAT-owned spatial-to-spin conversion with a configurable cutoff.""" + +import math +import numbers + +import numpy as np + + +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 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) + n_qubits = 2 * one_body_integrals.shape[0] + + one_spin = np.zeros((n_qubits, n_qubits)) + two_spin = np.zeros((n_qubits,) * 4) + + 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_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 + + one_spin[np.abs(one_spin) < threshold] = 0.0 + two_spin[np.abs(two_spin) < threshold] = 0.0 + + return one_spin, two_spin