From c8d31699d33a486bd3286c30fca831b158794e05 Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Thu, 30 Jul 2026 22:22:45 -0400 Subject: [PATCH 01/14] build: add a Rust Davidson eigensolver built with maturin --- .gitignore | 5 + Cargo.lock | 358 +++++++++++++++++++++++++++++++++++++++++++++++ Cargo.toml | 22 +++ pyproject.toml | 21 ++- rust/davidson.rs | 176 +++++++++++++++++++++++ 5 files changed, 571 insertions(+), 11 deletions(-) create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 rust/davidson.rs diff --git a/.gitignore b/.gitignore index acfd4df..301b3d0 100644 --- a/.gitignore +++ b/.gitignore @@ -134,3 +134,8 @@ dmypy.json # Pyre type checker .pyre/ + +# Rust / compiled extension artifacts +/target/ +*.pyd +*.pdb diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..7e43789 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,358 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" + +[[package]] +name = "glam" +version = "0.30.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19fc433e8437a212d1b6f1e68c7824af3aed907da60afa994e7f542d18d12aa9" + +[[package]] +name = "glam" +version = "0.31.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "556f6b2ea90b8d15a74e0e7bb41671c9bdf38cd9f78c284d750b9ce58a2b5be7" + +[[package]] +name = "glam" +version = "0.32.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f70749695b063ecbf6b62949ccccde2e733ec3ecbbd71d467dca4e5c6c97cca0" + +[[package]] +name = "glam" +version = "0.33.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f22fb22f065b308be0d8724e3706c7fa3fc2a6c7d6899df4cad7860e7a75436" + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "libc" +version = "0.2.189" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2" + +[[package]] +name = "matrixmultiply" +version = "0.3.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f607c237553f086e7043417a51df26b2eb899d3caff94e6a67592ff992fedc7" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "nalgebra" +version = "0.35.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adc43a60c217b0c6ff46e47f26911015ad8d2e5a8be1af668c67e370d99a4346" +dependencies = [ + "approx", + "glam 0.30.10", + "glam 0.31.1", + "glam 0.32.1", + "glam 0.33.2", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "ndarray" +version = "0.17.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "520080814a7a6b4a6e9070823bb24b4531daac8c4627e08ba5de8c5ef2f2752d" +dependencies = [ + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "portable-atomic", + "portable-atomic-util", + "rawpointer", +] + +[[package]] +name = "num-bigint" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367" +dependencies = [ + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + +[[package]] +name = "numpy" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a5b15d63a5ff39e378daed0e1340d3a5964703ea9712eb09a0dc66fade996f4" +dependencies = [ + "libc", + "ndarray", + "num-complex", + "num-integer", + "num-traits", + "pyo3", + "pyo3-build-config", + "rustc-hash", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "portable-atomic" +version = "1.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d20d5497ef88037a52ff98267d066e7f11fcc5e99bbfbd58a42336193aacec3" + +[[package]] +name = "portable-atomic-util" +version = "0.2.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a106d1259c23fac8e543272398ae0e3c0b8d33c88ed73d0cc71b0f1d902618" +dependencies = [ + "portable-atomic", +] + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cd274650b21d4bfc26a0a47587962c1edb425f69287324355cd040c3ea66071c" +dependencies = [ + "libc", + "once_cell", + "portable-atomic", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", +] + +[[package]] +name = "pyo3-build-config" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c5e2a7d2f0d013342f295c048ad19237add5154a55b1c5a254c0ec93d4109078" +dependencies = [ + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ca85c467da1bbc8d866eea5deff9cf29ea5f7785054a17da36e65bda9c05845b" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ac53762fd065daa3194dd09337a38bd793a188100fd1a9304c4ab312d901771" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ca3a1557399783172dc5bf39cfca835157732532cba56b71d2292161e53b362" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "qiskit-addon-slc" +version = "0.1.0" +dependencies = [ + "nalgebra", + "num-complex", + "numpy", + "pyo3", +] + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "rustc-hash" +version = "2.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b1e7f9a428571be2dc5bc0505c13fb6bf936822b894ec87abf8a08a4e51742d" + +[[package]] +name = "safe_arch" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3a52ec151f024d703f9fd65abb7cbe81e7cdb39f18917a3a37e3014470dc7c59" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "simba" +version = "0.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f45c644a9f3a386f9288625d9f0c1e999e1acf07a37df35d0516c7f199d9cb2" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "wide", +] + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.13.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb6935a6f5c20170eeceb1a3835a49e12e19d792f6dd344ccc76a985ca5a6ca" + +[[package]] +name = "typenum" +version = "1.20.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6f5e870be6c3b371b77fe0ee0bafb859fa4964b4404c27de1d380043c4dda20" + +[[package]] +name = "unicode-ident" +version = "1.0.24" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75" + +[[package]] +name = "wide" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dfdfe6a32973f2d1b268b8895845a8a96cac2f0191e72c27cc929036060dbf89" +dependencies = [ + "bytemuck", + "safe_arch", +] diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..087d18b --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,22 @@ +[package] +name = "qiskit-addon-slc" +description = "Rust accelerator for qiskit-addon-slc" +version = "0.1.0" +edition = "2024" +rust-version = "1.85" +license = "Apache-2.0" + +[lib] +name = "_davidson" +path = "rust/davidson.rs" +crate-type = ["cdylib"] + +[dependencies] +pyo3 = { version = "0.29.0", features = ["extension-module", "abi3-py310"] } +numpy = "0.29" +num-complex = "0.4" +nalgebra = "0.35" + +[profile.release] +opt-level = 3 +lto = true diff --git a/pyproject.toml b/pyproject.toml index 76cda4b..d9e144b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [build-system] -requires = ["hatchling>=1.27.0"] -build-backend = "hatchling.build" +requires = ["maturin>=1.9,<2.0"] +build-backend = "maturin" [project] name = "qiskit-addon-slc" @@ -12,10 +12,10 @@ license-files = ["LICENSE.txt"] classifiers = [ "Intended Audience :: Developers", "Intended Audience :: Science/Research", - "License :: OSI Approved :: Apache Software License", "Natural Language :: English", "Operating System :: MacOS", "Operating System :: POSIX :: Linux", + "Operating System :: Microsoft :: Windows", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", @@ -27,7 +27,6 @@ requires-python = ">=3.10" dependencies = [ "matplotlib>=3.10.5, <4", - "pyscf>=2.10, <3", "qiskit[visualization]>=2.2, <3", "pauli-prop>=0.1.0, <1", "rustworkx", @@ -93,13 +92,12 @@ parallel = true fail_under = 65 show_missing = true -[tool.hatch.build.targets.wheel] -only-include = [ - "qiskit_addon_slc", -] - -[tool.hatch.metadata] -allow-direct-references = true +[tool.maturin] +features = ["pyo3/extension-module"] +locked = true +bindings = "pyo3" +python-source = "." +module-name = "qiskit_addon_slc._davidson" [tool.mypy] python_version = "3.10" @@ -110,6 +108,7 @@ ignore_missing_imports = true [tool.pylint.main] py-version = "3.10" +ignore-patterns = [".*\\.pyi"] [tool.pylint."messages control"] disable = ["all"] diff --git a/rust/davidson.rs b/rust/davidson.rs new file mode 100644 index 0000000..a7575e1 --- /dev/null +++ b/rust/davidson.rs @@ -0,0 +1,176 @@ +// This code is a Qiskit project. +// +// (C) Copyright IBM 2026. +// +// This code is licensed under the Apache License, Version 2.0. You may +// obtain a copy of this license in the LICENSE.txt file in the root directory +// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +// +// Any modifications or derivative works of this code must retain this +// copyright notice, and modified files need to carry a notice indicating +// that they have been altered from the originals. + +//! Davidson eigensolver for the algebraically smallest eigenvalue of a Hermitian sparse matrix. +//! +//! The dense linear algebra uses `nalgebra`. The sparse operator is held as a plain CSR triple and +//! its matrix-vector product is applied directly, because `nalgebra-sparse` does not implement +//! sparse-times-dense multiplication for complex scalars. + +use nalgebra::{DMatrix, DVector}; +use num_complex::Complex64 as C64; +use numpy::PyReadonlyArray1; +use pyo3::prelude::*; + +/// A Hermitian operator held in compressed-sparse-row form. +struct CsrOp { + indptr: Vec, + indices: Vec, + data: Vec, + dim: usize, +} + +impl CsrOp { + /// Returns `self @ x`. + fn apply(&self, x: &DVector) -> DVector { + let mut y = DVector::zeros(self.dim); + for row in 0..self.dim { + let mut acc = C64::default(); + for k in self.indptr[row] as usize..self.indptr[row + 1] as usize { + acc += self.data[k] * x[self.indices[k] as usize]; + } + y[row] = acc; + } + y + } +} + +/// Diagonalizes the small `k x k` Hermitian Rayleigh-Ritz matrix, returning the smallest eigenvalue +/// and its eigenvector. +fn smallest_eigenpair(projected: &DMatrix) -> (f64, DVector) { + let eig = projected.clone().symmetric_eigen(); + let mut best = 0; + for i in 1..eig.eigenvalues.len() { + if eig.eigenvalues[i] < eig.eigenvalues[best] { + best = i; + } + } + (eig.eigenvalues[best], eig.eigenvectors.column(best).into()) +} + +/// Stacks a set of column vectors into a dense matrix. +fn columns_to_matrix(cols: &[DVector]) -> DMatrix { + DMatrix::from_columns(&cols.iter().map(|c| c.column(0)).collect::>()) +} + +fn davidson( + op: &CsrOp, + diag: &DVector, + seed: DVector, + tol: f64, + max_cycle: usize, + max_space: usize, + lindep: f64, +) -> (bool, f64) { + let dim = op.dim; + + // Subspace basis vectors `s` and their images `A @ s`, grown one vector per cycle. + let mut images: Vec> = vec![op.apply(&seed)]; + let mut s: Vec> = vec![seed]; + + let mut converged = false; + let mut eigval = 0.0f64; + let mut prev = f64::INFINITY; + + for _ in 0..max_cycle { + let s_mat = columns_to_matrix(&s); + let images_mat = columns_to_matrix(&images); + + // Rayleigh-Ritz: project the operator onto the subspace and Hermitize. + let projected = s_mat.adjoint() * &images_mat; + let projected = (&projected + projected.adjoint()).scale(0.5); + let (theta, y) = smallest_eigenpair(&projected); + eigval = theta; + + let ritz = &s_mat * &y; + let ritz_image = &images_mat * &y; + let residual = &ritz_image - ritz.scale(theta); + + if (eigval - prev).abs() < tol || residual.norm() < tol { + converged = true; + break; + } + prev = eigval; + + // Diagonal (Jacobi) preconditioner, clamping near-zero shifts to `tol`. + let mut correction = residual; + for i in 0..dim { + let mut d = diag[i] - C64::new(theta, 0.0); + if d.norm() < tol { + d = C64::new(tol, 0.0); + } + correction[i] /= d; + } + + // Collapse the subspace to the current best estimate before it exceeds `max_space`. + if s.len() >= max_space { + s = vec![ritz.clone()]; + images = vec![ritz_image]; + } + + // Orthonormalize the correction against the subspace (modified Gram-Schmidt). + let s_mat = columns_to_matrix(&s); + correction -= &s_mat * (s_mat.adjoint() * &correction); + let cnorm = correction.norm(); + if cnorm < lindep { + converged = true; + break; + } + correction.unscale_mut(cnorm); + + images.push(op.apply(&correction)); + s.push(correction); + } + + (converged, eigval) +} + +#[pyfunction] +#[allow(clippy::too_many_arguments)] +fn davidson_smallest( + indptr: PyReadonlyArray1, + indices: PyReadonlyArray1, + data_re: PyReadonlyArray1, + data_im: PyReadonlyArray1, + diag_re: PyReadonlyArray1, + diag_im: PyReadonlyArray1, + seed_re: PyReadonlyArray1, + seed_im: PyReadonlyArray1, + dim: usize, + tol: f64, + max_cycle: usize, + max_space: usize, + lindep: f64, +) -> PyResult<(bool, f64)> { + fn to_complex(re: &[f64], im: &[f64]) -> Vec { + re.iter().zip(im).map(|(a, b)| C64::new(*a, *b)).collect() + } + + let op = CsrOp { + indptr: indptr.as_slice()?.to_vec(), + indices: indices.as_slice()?.to_vec(), + data: to_complex(data_re.as_slice()?, data_im.as_slice()?), + dim, + }; + + let diag = DVector::from_vec(to_complex(diag_re.as_slice()?, diag_im.as_slice()?)); + let seed = DVector::from_vec(to_complex(seed_re.as_slice()?, seed_im.as_slice()?)); + + let (conv, ev) = davidson(&op, &diag, seed, tol, max_cycle, max_space, lindep); + Ok((conv, ev)) +} + +#[pymodule] +fn _davidson(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_wrapped(wrap_pyfunction!(davidson_smallest))?; + Ok(()) +} From f7072555ada248b83c6b92dd51a02a5ba536c16b Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Thu, 30 Jul 2026 23:31:15 -0400 Subject: [PATCH 02/14] feat: compute the extremal eigenvalue with the Rust backend --- qiskit_addon_slc/_davidson.pyi | 32 +++++++++++++ qiskit_addon_slc/py.typed | 0 qiskit_addon_slc/utils/davidson.py | 52 ++++++++++----------- tests/expected_fwd_bounds.pickle | Bin 37490 -> 37490 bytes tests/expected_fwd_tightened_bounds.pickle | Bin 37490 -> 37490 bytes tests/expected_merged_bounds.pickle | Bin 37490 -> 37490 bytes 6 files changed, 58 insertions(+), 26 deletions(-) create mode 100644 qiskit_addon_slc/_davidson.pyi create mode 100644 qiskit_addon_slc/py.typed diff --git a/qiskit_addon_slc/_davidson.pyi b/qiskit_addon_slc/_davidson.pyi new file mode 100644 index 0000000..3279c6c --- /dev/null +++ b/qiskit_addon_slc/_davidson.pyi @@ -0,0 +1,32 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Type stubs for the compiled Rust extension ``qiskit_addon_slc._davidson``.""" + +import numpy as np +import numpy.typing as npt + +def davidson_smallest( + indptr: npt.NDArray[np.int64], + indices: npt.NDArray[np.int64], + data_re: npt.NDArray[np.float64], + data_im: npt.NDArray[np.float64], + diag_re: npt.NDArray[np.float64], + diag_im: npt.NDArray[np.float64], + seed_re: npt.NDArray[np.float64], + seed_im: npt.NDArray[np.float64], + dim: int, + tol: float, + max_cycle: int, + max_space: int, + lindep: float, +) -> tuple[bool, float]: ... diff --git a/qiskit_addon_slc/py.typed b/qiskit_addon_slc/py.typed new file mode 100644 index 0000000..e69de29 diff --git a/qiskit_addon_slc/utils/davidson.py b/qiskit_addon_slc/utils/davidson.py index c2c6171..8c4a558 100644 --- a/qiskit_addon_slc/utils/davidson.py +++ b/qiskit_addon_slc/utils/davidson.py @@ -19,30 +19,29 @@ from typing import cast import numpy as np -import pyscf from qiskit.quantum_info import SparsePauliOp +from .. import _davidson + def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: """Finds the extremal eigenvalue of the provided operator. - This converts the provided operator to a sparse matrix whose minimal eigenvalue is required. + The operator is converted to a sparse matrix, whose smallest eigenvalue is then computed by the + compiled Rust Davidson solver (a diagonally-preconditioned Davidson iteration). .. note:: The current implementation is definitely not optimized in terms of performance. Args: spo: the operator whose minimal eigenvalue to find. - kwargs: additional keyword arguments for :func:`~pyscf.lib.linalg_helper.davidson1`. When - not specified otherwise, the following defaults will be used: + kwargs: additional keyword arguments for the Davidson algorithm. When not specified + otherwise, the following defaults will be used: * `tol`: 1e-6 * `max_cycle`: 500 * `max_space`: 12 * `lindep`: 1e-11 - * `max_memory`: 2000 - - Other values will default to PySCF's default values. Returns: A pair indicating whether the Davidson algorithm has converged and the obtained minimal @@ -53,30 +52,31 @@ def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: "max_cycle": 500, "max_space": 12, "lindep": 1e-11, - "max_memory": 2000, } default_kwargs.update(kwargs) - spmat = spo.to_matrix(sparse=True, force_serial=True) - - x0 = [_random_initial_guess(spmat.shape)] - - diag = spmat.diagonal() - - def precond(dx, e, _): - x = diag - e - x[np.abs(x) < default_kwargs["tol"]] = default_kwargs["tol"] - return dx / x - - converged, e, _ = pyscf.lib.davidson1( - lambda vecs: [spmat.dot(v) for v in vecs], - x0, - precond, - **default_kwargs, + spmat = spo.to_matrix(sparse=True, force_serial=True).tocsr() + dim = spmat.shape[0] + data = spmat.data.astype(np.complex128) + diag = spmat.diagonal().astype(np.complex128) + seed = _random_initial_guess((dim,)).astype(np.complex128) + + return _davidson.davidson_smallest( + spmat.indptr.astype(np.int64), + spmat.indices.astype(np.int64), + np.ascontiguousarray(data.real), + np.ascontiguousarray(data.imag), + np.ascontiguousarray(diag.real), + np.ascontiguousarray(diag.imag), + np.ascontiguousarray(seed.real), + np.ascontiguousarray(seed.imag), + dim, + float(default_kwargs["tol"]), + int(default_kwargs["max_cycle"]), + int(default_kwargs["max_space"]), + float(default_kwargs["lindep"]), ) - return converged, e[0] - def _random_initial_guess(shape: tuple[int, ...]) -> np.ndarray: """Produces a random array of the requested shape. diff --git a/tests/expected_fwd_bounds.pickle b/tests/expected_fwd_bounds.pickle index a0361fc8bc88a7f032fb7f8e2a9fc40ae1ea51b8..b5991b11c322e035382f77d6e2a272ecd19d2ae7 100644 GIT binary patch literal 37490 zcmeI4cXU-%7RG5I0YV8dC^BLJ2gQL!v!GAU6~WMn*v3JmD@Ywkz|cgBBTNLvhZL1E zjG%*}hylfhASeS!hhzYiB0NJ7DS}H(-UQzFdn`}LGOihh@HqK{%{h0UefR#(+2!tg zlM%%xJvy<_>7QZ#Yg3B!zR90bsZUbxe#x~vbm`l_TkYOSy%KwL>fgUpvOguZL;r4F z`t})+*gvUDqW^h+O4Mmie@b+pt|6A^{VfWO@!wgvv;VLDw*G5d6mE4>kum=6Q7KVf z6O;ROJGG|BBRBiUJ{`JTG^Mc5-@0|{=E0nw)Bm(6l<4m~w$ZzBA)Va6VFWz=4;sGU z{&N}u8y9=*UMSyL%<<$@gU&IR4A4+hW^%h*5yT!=cq#+(xk=cq%#A&`TBKc|C}(@46}-@RD( z=u_>W?j>jjwK>%ciW}SL@$wF-d?Mh951P+q5yyu8&_dA`u?bYu+8?|fc$Pk26$ zd>^g@xDM5g(IdYUeI0cOoKFNi-5iQ8jG@0aAHHt+MPBn^e-eY27e$^;^WxXMBMt$F zz&Sub?Wr5i(w;i~tx`^V%6)y53wQ6$hfb`U+&O4U%<18K_s2Aoi{s*$@ZJ&*C4Y&7 zS6X9m5*&H)M=Dym>WpGe(zXO}}DKL~{Lh#fzPoB5sB%{(aU zZHhdb=0{n_;}G~iA>hqdT@ho|6$|s2UtrZCJYF(iR&uCIWQW)3wfS>)4+|LTm`63{ zI&xSCANs-9MeC>sJ6dw!v*|qPjh^pl94mdDp9(QY{+9^knRhtXoMX=G+j;n1|7C|u zeZKS`J}Bq7l+G^{|4EgcH$;(-U&D{QGLnN_=J99rMUGAJSJrc1;tKwSziGWA4uOb3 zK<%9dXKC-aua9!?Fkt`+2Vq2G&kB#QMLg}QvWD0LyKk0yQA1T0r zi?8gcD+bkH=)+!ch*LVtR{IB>Y<$248Y%cVI~)QIf!qiXC+xh!yq%9Ii=mdpaCy(; ztFGsoxAPkJqfF*;)Q(fggR;)1mR@e?5I8Rg(B>ika!xy|ZYgKATi_?3!;e4V5BLxD zjEwW3$Yq?RyuZYH{QbQ9{D!xFcG%EHDldIlhd$`fyiL)E$7%G05B;Ioh5pE=?qHty z!@T4<;t)8E0C|YImpsP#VV?7bU)l45)1%G9d8SRuI{2U^$4zh9FnNg&ex64vSbhZm zBKco8`9J3#eQ0;V4@Lg#T92P`9lr`w-5(s#AIdxw93{^Yhd@Llp!QV#v$Ur!=C6-P zx1n}Rnprog=TozWF3k0{2xY;q@L#Y;2)lZm4%2C7y>K9WmYI35XpzPNv-{GnS2cZq zXBeNKHeHvr>ANH0G5n>(|Bd(wyP-62M*f8Y82E(JU7Fb-e7xIAHyl^|10UrlX(fbX zz)*(DP7i)?LN54l4v~MM;N$dk2t*PBl+U(5;;+n;FW?V-)sSYsP<)_lrOm_o#DPhN zW_J&oDP>a4)Tu$US8?*l=*LHt`zdHzO`Nmf`<$;9%PNlYcsW4p4vUZSd79+JYaQ_< zU2>x&hx-$!7-LW9F+I&d$paf@4>XeY+S!tS1UOgZaq<-TEL10@8O{m()D`ekPcfgS z^Ho~sbVvEeiY48#)5Mp3vq7=srL4%wGt-~)%N$z0D z8L4$Q>pp3c+fQ<6b5SpSDZRmKv-Eyac5IW~O=agMUB@rM4}V1taRv&m#3TGR|2ykA z3jzYP{irk9Pxe35wmefmQ+L|&nsc)&yWO8YeC8eB-7OY;@2!2+_n7MYfX>|x(F3Br zbe;Q$_yghkUKn|_ad>Tk|1i(%X=qQie`pg`wAw|ywuPVj&<0{2yShuygVGy2u?st~ z``_Zje~<&NjKLTAjLGNt6|^9n8K;9oK<%A+XKC-aua9ycj|`<{+@COepVI$PmDi!N zIA91HHve{om2}bEcPbb&=!d>&MfpJaM%hNWL^<`9@EWD^XO_w+$?-<1oXHd){H;{h zd?S1cYS;uH$}#X^-$@VXPrAZ>@MHcGE4`6cj|m6l?bdkSrrpmiY?$gpKKFOjA>a^* z90c}Qd56cg{hFuzCl2t~OMGFTvY2>7xlViNCBX6gfh?mk#CtNAChl5ADny65A)P5JZ^G+UU{SJu3mG4rgQuY<2P)}HZ4_u)JW_Y zTkrLtS-7|FeI38ZHv6RSh?A42)Vd~UHoh~e+~xs6vq9|_zB;I&_D`(Nc{8i6#A9sX z&A)sN)?4kMAFrSF&s!eeW!5|v(`oF)J%-2ND$+;xd(ffMXS?*FJw;oI*J{zy2OMaZ z@p#QV_+eLD;gKNPMs$zx*k#c*8iOzOA^Ya&xkUi~As%u!{Lc0x+Bx<-ahZ6GAK`cQ zJm<&uuN9s8Z^+CJnzcHgdv!jfmq#0jdX{tiiPmwSke{DvD{)RazqFNRrM7tOi(bpZ ze`__AR|1CDAuA?b|HA&$W)%Gd2#i)o3g2m!0SX zK72h$`;E8)E_QtnE>Qdl|AMm5jyeP)0RgqA>Yk-NbuoN>L>TvF{{6+3zh7vIJbqc? zmp#TC!d{xnj=PkOwdG#5Ya1P1G&pG1>7!-B<2vy-6F=$UcF7M83dt{MM+Y2kk)Fkc zL#D#(km?h>y3LlGVuMahi)yymY*Dzr@vXN@)*GK?It^_%B%x16_-({K;ZnRoY|)8* zf`;&Zr}*26ANfy6eo@J{aVsi)&KGXzV}wsJ>GNjc-t$(*e&pjV%~s*HOljTOmf0?rIM!y&5b>;e0P2MY?B_7^x}`5 z&cxNCZEkO$-YsZ;@T$BYHRiC1m7Kv^r#}_WiI$!HWEYP=;7vRPM{p(|5=V(&(1LuHB1t#)iTztC+|?t zG*`XTRQ&hLzlQ5v?3=Q~SMS?Jrt7fXPtNTz-Bi`N`FcRL8SN`CH(PXluG&P?U-|it zx^Y=%+5YtBhSv?4sBUY=+&?d9E^EAP>f)+d=3u$&C)XboFb!+|9OEAyG#PQJ9~>W% zWrFv#eWGr|40BTLo#QQc_AXpI+pHe&XtTukGEDVRgFDo!cF;Vq_44&6bH3uMawPWT z$b{o&{FZWY3o2%tk6v2Wp;zA%=Jt$x-x_l~XbwqES*;r(z7vv*oF!VveJV(9L&=#f zyMF2#vg0P{JyUj8kX;94$7N& z!*AJd_LpzvGxPqtcUw2skMV(1}1b5nhweQ?9t=PClzA~NWPEDHp znyndO<|vp=qOrVZ3m^Pgz` zQ^|=_J11JSh2(xFzLlaiRL?e-oHn8#NKQ+wt0}&+lDA6hn~ITC2|ru%%Y`W;*g2}f{WD|{amuIxic&EKkd;wscW2b>2e{uIUjiJpry>119! z2&lbN=Pd1=i{a~|W^$wEmOI@gH|6aLVR@s<;?PHi0mH|?l--1v^=fadP&r&#_2R-+xEW1v21P6xNK9yD0Y*aZ#8)uQqpJ!Aay-+{Y_6ew*si$)1 zVYPjZs;pU~GKaQL6YUFKdeUD@>DSuY$8+t-uj+o|zHGm5?etV1&-uB=$H{jHT!;t|PnuZq zq=m-h4Qz|oJbuPJWgX=@Z*Q5WykkDz$~&}`YISd!dgI~@^O@pdImO*IinDS4)+3Kp z4;a4I?5a3gMe!lYYAb!M_RluGj-WkMN9`cqj+WPJ4B6uyQ9J0!C9yv&d_3DcrFgzn z@%(zl@d_&!Je>0KnLpi6`@fp>`Cf7Y`YNfJ+DiNqP53RX+Dd$_(OCEp-v-LA-w6+J zDXTGgg2yBB?azA6LRn3mWZxV;w+Kv7UGcK~CR_V4R`XLePrU*EY|XQ8d?nPJd@D*@ z>7d#f_{CVAmml;Pw_Vq(G>RFR`Q&kPXnX5Nj?N000i}0Mm@znDZk-fAzfR(DQ@X~N z<4Qe{VJ0+an_aO^z+By{h&j;hxVc8>nYPj~wS!Ko{nJFRvns1S)J|anvS*GV-D&2ES!xQGvSF2MKEa_P95*vyGKJ3*!Y83bt)r!KeyMNnkx6Y{E*>x` z!sl7-Bl|^LiMU6bka$NMiN|C7o45$&@t64ASm%Mduc^k+yEHzxr@7!spC`MiOYs-G zUS0R{R=;1U-2s_^#=6~eg`JDMey8`|EHSbJFP(5U!UAp?hx2yZnsHZo>#=!g-a@s7kQAtS;E z4j&dBGa`0ijAOPVv0?$IBhh2nphA|}jvi$a9FLR@bG+v0<7n2S?BCm$OK?O~OsqI4 zCO#_s#+-7`cW_LdQh2+3Vp+ox930%;mH$)lt4Em_N7&TP>wOA!vVU6$nDQ^Eb;JHg z0RlQM`n3&HAp+n6wQ;c#un{Od2$*6j6wA4W{@QWo6-SY0QyT$yM!?jJ;~L6suA#+^ za~}$yPPwmBozK3u5pV|tOfdqAOdo6v&Xj9t z8I8{x*ZGRe=Zr&jd5r5+lh3-pL}i*yYy|#L1WZ{SiZWb7t#QpS7>Anp2LnQ_apr4$ zCVWtx4}K_o#ZVcxiH(4bz;B6ww5QtMq&-#etx|q_%KrMO2;k<5XD{zh3UevXrGt(2 zS5mDuOc72u9Uz>8gV)b-Fu^5RaQX#)w?nBodsC-I1h!h!(-4DAXmR{*2aIiyjev~+ z2v~VUkDtZG)#AoQpb3Rfr>vt>6tby}KnXy=q;BAxicZ!z^@ToOMd}b9D_2VWfgEUd zVZv4YK3Kjq{bi?GBYaQF^>w+A9OR-Ga#{arnTK9nUlBR*k>{bhyj7xyx7YziuE?>8 zjllndfR%SR*PKWF{K+EbYVBMyj(;(_*i@hdgC)e=1D?aYB&$b;l0>32!y5A$axHv#r%slvo zTIK!D=01mKPak8oo5^F0bFIf`%4EuQu8|Mb=b;Z&)62fJ5%@z8pl;z@l82v>whev; zr7Z+Ml*doqFPNvdd$c;yr0vAxF8muKmOt2O=LJO_#d_S=^=DkbG@8eA@;-8Ac5abBv83HT^+-75xs zjQ6_OKE}tL;G6`dPm@9_=jy*zz&4Gf$e!+ood4BPk^5p`Z z>cF1{)$g`3P4U(){j27&U!Bde@Ip?G@P8@%;2bXU_x&8_5|0~cOW*e`7t_1eJ~gUH z8h%R3UCPHCNr&HcT>j03&TkY+SKy3%@GtSe*b3SRxFG_Ri+VhyY{XB=7sS0m65qa% z_(u6k8B82`DK7ThI}tARk7_Au?p&8jl{k5J(x~w@uDMkASxZ-*&;O(242h$Z<8%IIn_3~|Dk-%;F?=z*@T*>Rr7j} znTO+Ds*Cu+i;L?eJC>eQpM|t(T)D1G%@jX*F6j7}#NCZ64}R}b?~5G#2Y<>H zxyXr-d91Tb+K;)Ee~VmI{w6h8$o5BWLH|*y=-yG=gGcg|DfLIa zm;qJ%|LsyMzaQRjz`i_%9wUEDn$xneOYL1VvBsfjm!kdhFKGwS{^6VlOIvBEw3X2N zdwC6ppU_tNp~aF9J3e#L&%q7RSLJnqJbK6^wT_Ate5MorLwxS( zDF5!}_gWa7AL2Q3H-(Nkn3?yJzo3V#V^aSv+jQ60)oYsue)X};V||{hIQP}$t0442 z&LNeQeD(c_27LGd_t9s+_zV8GN7}0R$J?o~WA0ry-8d+G%p;FtC;FjB&WC$Ldz~0> z5XX+b74_|z0evidz!kjU0~fA$9UpS_sqPJoz$uk-4(|*z){blIT0A${V4vUzXZFpe zHUjQ|fV8LlZqlAA4PPJO-u{_y&);$1YIS~OotR^TXDG^Z-U^YP31jt}b-uWEv`ej) zkCq8HdO9a;1c>}>kzdk|j;`lL(UUM4(&%;5FKC;}in~0oV$eFPEb~uZ<>fykSM?bC z%$VN8GL-%}BX*jUO##CHnD8HzbkJYqpB4G~zGVHuxRDOp$UdUaf<2Q%yy}iLzz@7g zs}IFR4e0uDuDaW|L+t!9PBkbhtHz^>O&|FXAEXAj(Sk{ASLH%zh_Lc{k`kVRGbAbzb6F zr4)bF^3+k4T)%C`y}oIRI9FcUJBO$DTDN&it|~0^ohqt#bMN>?F6H}Tz4Zrb=BeG4 zW9MH9&9J`eycYaW|D)k9b<(_dW5t9EiZ;;*nMZ#5Qsmw*a<0g{!nfB>wNvCGhxSo` z*x@O9(;oK}yH<%E;D(*p1s+i16!^B4_{lZ662G`Eou~Bow_4vH@Mr(nC-#SN{NkMW z5%rPIPu;M6&bNmDf7hxeFQ>oy-rz-wJik4<;o|;(-ln+!mb8hI-o5^?pHHs( v` z$N4!`-k`k+PkiK3H3E;!UH6w<^+SyoNdY6Bs=4p=%8p4c^{r3JrkwG)>RR7EulPNf zp-xJB=Ty(Mp=Dd;sZXL`2#VQ|p(;%r9n!Ml85MZA?tvd~{4YcLqXo$K8S}bxI?_eYflb`_6ts*;k!%E|5q3{QXbtowXu3 zgC{tGH~XfK^E?9H$eA_k)X5Efo-kf(?=j_aLN|l?><50zJnlnjuNYEi){?rEc4zm| zWA`R3$oIR<8F`xlj#M9697k zuGt6n7s~!|4%o*#>$Kb2mNW#Uz0>+8?VZx_^%3rb8?!yQFX0RK+b(4>-E7{%P(IV% zNRhS=dy6~EC8oVQqJv?@`3Umdnju*N;&nYv~ehZMoD=jIGgEG zmQkkl6@JQ^jZ)?yU*8w>Vt@2~WIt~gCETZa9QVBL>T93}dZHiUowP*yqU}Ta@sYH{ zIU*f_7ikhXq-W%_E@_4J!3BD|@UZ3D2)F|Ry5C?E{sU#4_{}`ViSPKOPg3TxraM|& z?OH!+cQ6hg{>wUbBYLK^TbH3WOFZQ9c&o%&Z%6QitVT{n9Oxx+w5r4t9tXXpt#m-z zKYV>ad#IVTgLqxhP+nsYkB&=QDf^x(msXF;Q}J3nZz_+iUOQGkm$=ZScx(s~eTZMk z;VUWHO1vK7ZLWv3m3Yfbc})DIeE=Rhg%4by;L3h+&A!?6_98%=2fv}NXrR>_>pJz6pN*;u5Edta{@0#iO0pAID#6a%`ISlNsu@ zHhuDHwRWm%k>%9s@Ek?H&XM-Z326tNk@n96@;Zz5P#0+{@w$t3q0OWn^tqYxH_r{Z z*81Fs9y`7#rFY#kivL!Ux26wmr3KP1yCQ9+)}jyXGuleKEhY)qIdou^)1YLv)95(PL&||uwVEWZ6|Pt67Qh+H|rA@q3l2BhHKVgob$nb z!L8H~>~A@(y661e0S5lXc+?h7cw074Wb85RW@$PnI$c)$kcZ zup5fuUId)D@2Eq-A>a`By$E;`07V|x(32f!UL=Y zL!cNCPr}he(%JcgB-ugMy6yi8D2E5mWLzUm+ z;8oK#I0=qO{G}DGLiUEz-a_AzFZ>FA*f-u8ice`>cW0MFpeP7L$|EN|u7`P0;;v1} zS2it*vW~|g@CPB_El^z%ZPgXj8}yq+&Dn74gYz##WNU2`8f z%!3d8;Onh<_*D-}4tzF!S^6OtiXAr9d`BDtzXySEeqnz}SlI8_$BeTtZND45bYImK zAN|`0#owyR-=O$Uy7v34qR6-BA+NgRaG!Df+4e)^LD3t^bFl;ap_=E2L!eY3p!QDl zGaMBp*HDLw_hZP-hgCN}eMQfujh}rlIa7 z58xM!v)}MTdt7jOw0YRiv`Lv4a-(C%r&T7-?9?*#9UuI>kJPgK2>FZUAARNj_!aut z?=8&3pI*^?{EYkfRha7eMfEHAB+n6tK&eJR?WtzRX-}QWUmsD{huSTfX6=**9-T9C zQK7FzC<{)&Ut+fqcJ)6RrZdet;XoRwK6`iha%}@<``LXiYX9z*Fh0M0c}2>XKjwz7 z;jb$G@5N7g3Z;Q#@=FR};1f!BnI=Q{c-K{HIkD1vKFUwha|p+Pp$wIs9{k{hT<~EZ zBEO{Igi_s^q~1paWZ$@ zgQIKw5;QkXp0{vs!B>m1ildanLo{!n_$Z%eNKU-w5l^xtx1!|me4AoV=rujlK(X7V zrD?65Ek#FwePicw#zS?|vHZ@sT~|;~F`g-Zsiu9px#s;9Qf}B{;$vTHULkdv*`a+r z?y;#2s#FP>EZt}S*PK!2)iE`9nZ_61QleT7%lIy6kkeK3@XM-_J6v+cXkH6FCsT3ZPsH8@$#_??+_E zM%mq7c7Cb*;D#SVkz-TpHT;z8qOYvuJRK0AEkT{hda@rvZOg*`MZ9L+*)QAjyZ!b3 z$KUkb)N$cnZ_SIo`&HiuwC^^H?h;MZeV!lU4;Z@*L>uQ#wP)}j#yOsbKA`pwZKB$0 z*C3Z;TllpW;jj7`$F6?TbGP)yPVB-??Eatl@EhdV6nx2tT$9hi7kWDE8K;9oK<%BT z$7%1lua63Ek8Gu7Ue^h;w<-PaS9u)@0|7(Wu=z>f(85nTJw_<6plB8P2s_izj-c&8 zxkNejweT9N@@J09D9Q20s+`FYAN)6}todHy_jC=LtRED7Sa;F{`jf7(AN&|U%Q_7q ztsW2#$lI>#m%rTp)S{N@KIHR!M;!tVfzpG(w-NHrPK{IklaF}qwet?Auf!Y5b=pJE zDLyk#2zt=+dd)8P-^*r=ToyEbE1vUuJ6P=|=IyfduVA&6XoC?q?J^pD!l$o9pf+j=l7rr{E zrS?yZ_IZNUR^m0b?E2pxgC(gQ^z)T-{&Vd;+sx|wqkB%6{H@_NxUTe({T_6r^w}hR zXiw2r;#e(8`hWxNGG4D42S4nh9<%AU!eg67SL+&lp{$#u#fv~BKRVfQ;xh3VKf>?q zar}_`tox@u2WRBu2hAGo&+oK9*uS)asAt*7A8Q`Z3CRx>IyPqi(pH+2-tnc+6PJen z)@mrv1q{a_%O!WC^g)m4XF|vLgVZLfRCeu?57k-aTQ9y9n#b`BcA_80xwMsb%1-pb zzxaBP_8axQO~J#a;A+$2-BnJyL*S1^K<%j}$7xTU311%(u6;TGdUpAFi%hu(&rROi ze^R*ZlBu%eFDh5GI9|QWTkl^yJZRSHqh-S5TJa}{pYHbck{=uvlK;mp@`;}Y4%bT0 zO2Q#W;dMm)C%h&J%S?q~hh|hvSYkFPTwneA8&#T4$}>Gjb{o-kKz2B7><})Mn#Yu% zJRoQY?`_22P5j7zSn|tDzKvUX>2qRBiiZPzjPR)-eO@b@^z!nU4}6@`d?mb=DvjMR zbWo4@ck;}+4O*l;J0f5zo$WtAdh}Y~FQMZV!!N;oApUsqUm$+ubDBaqpQY!Z5Bs)w ztanl!0wo&(;u~?kjQsax$7!EahV#1I`l)G&2j*@ubuPYp%%h1}rjp`V?ewOm;m$0x z@}jG<>o(Y8hAGZPslBuFpIw%G_)(q-wV?whwd)mCQ)dRv)el_w*1DScCM!DS*`IqJ zi>u|kTz_{~-=O))tNU)n@%v1S}n@FYDsU} zR~N{x<+6kL6D2zb$u8m~xDpSqQ``i9J0Fr4p(T5_l^pF^C-Odd%^oNJ+T*NO2h}_6 zwa&N8zaG-Q*fDLhujvnqO`lQQADQ2Ornx}-=9?k)X5D?xGP6PV=c`RL^Mzk;XcCua zmhQ}Y>Y*k9Q?c)w@wdMmH0QS6IDN?ld1iNwE2lOa7BDRv{2J{a7c|*%>F*sJoo9l# z-28Bpmf7a8+B*k3ZAmKIDBr9aa$iF7+u5f6*x}t9*4u6F`0D(1hYP;qtD75hcud!W zX3~ZlaSLnbn-88_+dXmMA#;89t*?(i7&Ln%Csy-Di|>%+BIg~=<2kh?x25FFm0iE| z3E6Rt^qws{Yss!%vSXF(&XJwbw}$NMBAiDF*Sf+H+`;QBYMLy~aM(tZ3Fe}&?g zQ2dtlhnAJxlC#rKji%sAorZsbKjYvEKln3F`>#=(%`+-YT% zJ$j2C5WT~SyPv4eY;D!0LdWYkeB_$(&*hnAIg|1~taq#p)JfwXYaBUoYUf0Wwv*h? z#J61ZV%4)9B&Un$d(x+)=3OejSjk(d`R&EebCApUhmwOGe7%$)x!A??sjtyS}P!6^`J(M%Q-=SJt72##?EeI0?1an{|icPt@_yQ?qN`Q%eg1YVR~YPJ8D}`1?)ZR8G;xS*-Hsag|5UGz+zT0xD;ks@%CpZJ+%rYgVhwq3zR7 z>%!>;>92+KgR+jrO<&Y=CvLra$C{bxz9?Pyke>0Wp>`f&{!^uq>PlauC_oO|4>~Bxd{;wx}_DW7b zUnM1|t;8?UgimSJR^n@owu-C7x1q9&d zL5J1;X{Y0?bJQN{rnVAaVFXnt;$PRP9kj_yzpnDzdEw)-)$fn_bjYea^T6Nh+;VNp zLuSv{tJ1rkzdL*%eRys4BxknTN|n_fyGr`JsJ2oawUs!o8)xAYrMA+_wmPiGQs4aCDP5kg6fkMR=Si*OT&*K*CE^}!aq1P?NW327-^4{|qm9BEIgGUu8`lu{~$tG9P*G){71#-I-c& zd@(#!593Gdmxm;;tN7T*lfAlo8f?z+(N^NCA=*D3)V^siN`1}OPWaWhsej(E;y{Lv zcF?Fd7T-2@$U)Np|J89T{9RR#j}%{h@u3gxu4PtRiR0Wh(kEVWIKJh#2*Bqf)$zB> zZt^wwJS~5s-3E>eFR2mRC(w<6b^ndQp BXwd)w literal 37490 zcmeI5d01Cf8pqucaYLJw>VJ?iN{t9=; znrs^t7hKXt&D=7#GR3FR#Gd?xxFI5%@Ee5pb9*}Ydr+TDgv%d%J?}l|z2|(-d6)a% zON}e{&MUJ@<^8;FZyHx-Xghn{{llV$hQ+k(KX7<>aLb`lp^<|F!@~n(>~WR*hX)TF zJ}e?KJZfO1eYQQWd|pp`oZGNL`CVq)yO)Z!KT$f!{-(XRy=nK-|7=$#)*ezmuKb|L zn6TjMW6HeT-ad6o{^hc9r44&PKtMNV?oZyY?xiB_K~p=e_spM@_1k2?CjSD>3)VmK z7|?mq+ZIjv48RL$@nSJxF;H?Cu*p(ER!$qdRXO8aA^I$|7%0vR*tEsb2DwZdbS39H zNIxC9uA|=1YFi8x2L^1i2#__=1`CsO9c0|iBI7Ky7_b<~Wx$?W!O1NowY7(o3vs(% z4hnIn7}WE6G01ajr;#-*rf!k}oAd*3a#-={cD*#1`=j!H=OQa}d$E>p?eHBk5;Dt|R?)WSk4s zZf`X$1}p~t2nM7))#e7}sl0EMa?4ZJ*GC12wyk*a^8WZBr*dC9*jS&HWR_tH^em!J%2gb>Ww~%^ z1Mwp|;+J)__QGPIC@|nE4(f8wuP!GK;qh{%e9+3m`)4-mbaSCDhO zO6&|f;=ao@J4Ihbu{-Pyxrb(t=+Tk+=*WEZ`Q5uaYt7tZpdbTM-f4Y<@{aZOQ9-(J z2OvFgAGjC&I{k_}lGfmUaigS7xOGtHjs7|o7v4w{BieU=&L?=K!8+lC>#|<+F^@L= zbu2FYSnIGDuox)F0PzH*2W=2JX|XPs>!x4X+066s?CGNoS90AyqYt@`#6KMwrz7*V z&|=_jGC^wd(F1p=F_6aq=fB~i!DUk0 zPBYIH=cwvv-%sP#y-`5k^J2TiN{u5-xzK5Tb;Nnb&d^6c&h=*x2PRj3JI3JLyRw(7 z`2!UFu(v8wPI_q9q>MLguNmBz`+&&F=RvM#?wx*o-hfVqSH1dAe$m6*xLLcRA8}pu zSzs|xtQnBFomH3ZBt_uklmlH_)1(#lS5vfE{Cx z#0AnuU9R{4TKpSnD`_x(Z)xWAHs=0HWDoOn0*-0bERln*~-DWLac|P}# zj#I^tl9orvxbxBv`E1eCRmL$N{BOGGL9Sz<_;Eb~whVI*T2=y0~tFed$T{W&hTVD%NqTnPMj|`5zx0w>wR};?{A*!JnM! zBhiEXU{Bei7d;^|j{EEqy@8^KbzyJJ8(y}`yjjNq@e_uulXa61fba`n=+pNB9>HRC z^h(51*yjP-^9XYGANz78=k=%evc*YBOU4;-@lVg#7ZYqSp8Q@%&J+D8+vs^L`t|sR zy%U#+d)Ni@r5}2AJ><)jfzZQv&ObbC2|C@#W8+wZb3aJ<+Z8&dMVljDrkr(QS2}WC zNA&7gA`Yds$Rf{xly`h?P~Ne=J}S6-q{X-s+z8?SoTS$;$YVTiQSX1RxC7P)!U%4T zb&-yIEon`2NvGf|NYbCTB#jb1wrEMq=ttVqLDCxZ5^lJEQ9)Px2tR}?`Vrn(Ki82H zRtRswS~w?s!poutBQ{)F7h{I0qQZ+c4~qedfj^FceOfsje~jP64Nsmc~ z!=yZff2Q1oACi^WcszCcs>$Lp+kPuUah=E8XwgTxAHRWbAC9-$xYU5$y(URe67Mc ziX0Lc(eRp^`54Q2;=IgTR_#RShfdYE%eylU9?4NQ$v^5w4yf$&U#D96)9}6n_T?zz z?ucLG=d@_#RD0Jmp>QdSx&WhLhQle`AQPADt=+ ze|w~?ihaD76gB$6b<>T5(vNY(QTzz=V;-5G><#RBVw@SDj=mfA!+#eo4SPkG;Yw(36Sh0(o+r7Hb#M%HN~ZSs}= zPf7oS5)S%_{ zqaA8cSbFtCVH3=-N%&kZVYq|V|I4ex$EPjoWAN1i2tUlHsGh5$cFwK&luQP&7v0a_ zs+@8VX(fKsd*+iCWjtX=f#gzNc2%HP90SV_1dI9=#{MSb7iHxb9j2s zb(^&vU4NiPj@n%@YJOH=s`*vt)quzP9SwG>lePypmXE!l zC=-RtIO5YcqW59ZlO^Nw-(EY^PSJ}V%16Gk4tJS1<#Bge*D6^D{*ZODF8Bc9m-we6 zaR&cIyJViyKY!Kw{=h%`$9}nzpOd&uzC(ZZle}U3oF5FIU)HL|ucwUvaPT5UoZlYN zU~#{{Zc|);SIR{3A6$Ff+cR7J@6^$mW4#?JXVBi*XFheR>V8M&uKQcI`nh`Zc;68Y z)y(T!Mf)VD`oS}CQ|7pAb+u3LH@qKBRVSsqbE-%3(9$h))aMbe_(yI?RTUffTl z8Rd7l&VgUA|1U%?RQ`2BuS^y7U3JfuwQ|&Z3-K&Sy|S_I=0HXPsuv3$-3Y&#YOePHyP^jPYhWw<(un zyBh2}`+?mtj_V+KsUdl0Ey+tMcXk^!W^e3*T)Rt~k+Zph!;Cv!M7~Aj=wTh?z2v)$ zr$70xzvw}qp7&Oh`JnG}xnF>cgNMx`*Y%K3%Wm1or@P8{?$3Pmd7@w6uLT`DA3axP zsDZ!Z`nhQ{&o%AfS3r&)9oYx=7i9nNv+U#TojPlj#Xvp-Qr>BGgYr&E`1%Mp2y!ei zL);hUw_Va=s@c4SA$_L2ktk&!(qYnkK8jA0^dJ32)tu`}T7cWRBw>A*qz|O`tcSE_ zqoh+$N*RarXQZS@!m}+!(lXMtKGL7GW}~Dz=-2m!d9gnr`*^1jcc0{T-2IxfkHP(M z|I80JPgo*+5oQQKo)X?TM}#BzdO)tDhwzMk?#tsB_vd_rcbXSVzs114%m8)*vKHE4 zVRHOCcIg?PwyepHR_0qtZz&fe*Zmv*q)teW#J1~F)n@UBJRWZqKkH!+n2_Gkq3{De z#gA4If5GFRhm@5LNco4a4=4{cm2wcTD;mga4E)h?DJx|>Q2Em8kvS?x^XEb|4Je#WCY~4GjRx5|98d^r34$f4>>r5%XoRD(R87cofBCoS3 z4|SHZ60f_2WgFO3%0XY*693`8A=~`pe6EAXj_*nARp*T2zm?>n&4;qm0x6eeNm;3t z%!l$BWu@Knbp`xjCt0$tO@3$3O#8P(C29P0c5iXfBlnN`7N)<`b72LCiWPn+<6&Qv zo!}kBp7D>gxj%jp)c2orP#6ktB?ib#TTH9wK7Y5bnU^AO68FvKhgW!;^0SiHzu7mkbS5p1^8O&E-_js|R z$%5qj@W^@Jy8gZ=YWL}E^2~mL@K2u6R{SK-2iO&ey`mq)?m%5$Vhqb#=%32K_3~6y G%l`x7y8HP6 diff --git a/tests/expected_merged_bounds.pickle b/tests/expected_merged_bounds.pickle index 1d11a7830e889483c11998c894b65721446a2144..fdcecae0e3c0e97fccc1d201254aecc89801a2f4 100644 GIT binary patch literal 37490 zcmeI4d2|$27RG@P5+E!A2E`0+;Gj5g*hC|Zuh2ef&PBYv7irdJV>Tf{J}l->ej9MzW45Z z_kFKABa2Rcc2c1W|AxEoPASr-mOG_l-=seMlWTVD+ApDd%|1!J6MJ?}Na&pGPKoH4 z(7kKFz5@~ylDa0kUvsB~UeI)>g!S#_Q+drDS7@xeN#QQ;zq{MI?~W_ns!frx?jE5j zq1_Ub`*%P8Op#|EaF2V*cRM7du*==Lb?X-1ynh$|iz}4q?lP|N;Br2ng7N>@!M7yt z?K#1;kfe_fXZ`by)qh-U4ovR#;_P9IO;O$75}x$N&z&<(k9QZ$-9Ox8nwR;YUeV;0 zCZx?n9WuImjh?qQ_gE0!V9g1$>W=#|qpI&Ql~U`N8izCd`jLaqOb?A+>bGCS_I~`h z(K2QD{uk!;oMC3jKm0s2pvue+l~$U2)DQ1KUw=@N8`2Cz+=Gu=Yda z`(~OYZ*JM(GHm16H*?ILx@I25sh#UT>i*ZVbDOSr$UgXd zXPdXJ{d=U0jg}Q&5JKIIO zv9-tOnsuj;?AGYfGPVBFOtVJ$u&nao2IbjEck59nt9s1$dS5r?(J1AIq^c9G+OaE5 z*ng_c9JbPHzWYa;rC%=3F+0mWS2U^79y28Vo|4Hkyrx#K+dtl1F4vTs+tc;&&*`Sz z-HlHz8FF!5J84CI8anobIeAmWPm2fVnirMNmn)y&t2|zQ^`fU!-n#f)XZce_{v4Mb z&&}QLX#VLQ6Z`7a-lym6FmpQd)|4fN-zopEca`m2TJ1@SXN z^Wx~#-L5)6Eb+(1IJN(z`LIj(qZ)?~%z7c)9NX3UnNzbpW*ZX4h=Ue$?0+iIBk*ut^dUbD+e zZCL!19DhBz;qy_e2CUC9z5Z75(Fdc?m}B3rPmQm9)PEm;Zk0cl?96(<%)Z|9yry{f zjbq!*_nQ0U&pU7b(x`T1j;Yvm=d`7_=9umBXN(P>u*ZD;TG#Z)qvL0JOqTF@Uiicp zuW_nm36Gg~a&nuuig`?m@Oee!$arN4XYyXU@MpXlYWy15;+uSl4AuCPpPTBMy05ve zkq_(o`p)M3XMMi3|CW|Q61i;iF7qP}daaoQ|B^6PPpYsQqJ z`?LD3bg}L|cE{;HPbPaz@(^F%C61UcBV;dLdd%ZVR@AZTyVG1>soq`PIUy}8*K5kC zer_?af2WvbIi|TJ_4O+05m#fT{CWG@(`l~es;7s4xa9HC1F}tZ;`jWfH!m*HkC)ys z=~a~;{fS!RIx_(@N=q;{w}wxr86_l@TDm zh%?5OxI$8IP~Q+|NZN@j;x7M1fa6gxAmGkppz?^Y+U_1jdxV{*vw9SV%)0+PnH4#% z@yC(D(9z5H0os7@rW`7xayVV(^1~{FWW}nYvVD`vi}KRnB>h1BqoKB`{AnV6Xe;bLFZ<9?cHE$$ zZWKcJ!*-JG%m;=0%wHYxh3klLf+kv_<%s198midM?#kgKH@VK*mP3}W7%gio&)7dI z``aa1e(;BRd!z7i-s>QcpAjIxk>@GH*e2QAi7(1I^uHCI(p~fk&o>smQc-kDG3Ch$ zI+j?cV~Rn_bI?SGMHBVbaYwRfqWYqVGISgwY%PvKeiFSCA-V~A3B4V%TTgaEb$sHZ zeLQB9?4FPv`qfPRC@;UCQ9o~0zgDXs;1i~P_E*0skHMRK2#(;)dPp86eX{a*cbZGTowiTVM2AFAbQ5iHtLDvjI_BE0`}0Jn%n*ID zLUe4X=#^VUrySMsRkY~POdWGU6P*)Hlr5TQy=dI3qKTdmO%$c$v(Li*1y6f?HL}GQ zbtd@{+@XJD14t3u8cz`wL`CxSCRI4zr-O}k1TuEYr3giXB%MG zC)`Cz4sJf<);^Ve^j|6CPMQ(R`p@>%-a{b=N|BXO;>o^JxN=C zJn8@Hp-)(YL1_oQ=4a7L8=edAJZ{nf(@XVal>AZuE#xryvrGOY$xatNS5kh2NrJ;} z`OWbp?chhh+NvIlmxN~GYd+`)?JLUM3$C zN5q}Io%v&rFTMuit2@4akzdw_iW}-#=J8jeZFr9_&OQEkVV*L7Ka!oV?NM*zGB*`ctI~*G*y4CJH-lbdSyh#^)8O5*UF#8vSY`gtn$jf-1dL@{N;~a zhlMNtfDd2y`rpU% zjPxyJXD$C%&=2&QeBVyLdMaI$23ZzLtK8YI?FYD=-TLdY3w_F@@1;+Ar@X{|z6VAQ zf)Cp&wvAT9Qqkj!edc1hVteJo(PMIKUmu|lPS^z><{|d~U+{7KbPxz81nlKH@kqQ9 zr>qy~Qqs%W2s%9AH`?9ocw+%n3$l-tltC!|NYKV5cW^c?vl zLv}-Dhxe1GxTZgp(Ma+%Wim3D#@gx2)g!>XVjU+gSkHWQQo3QDpifu@a9GA7T(L@}l%}|_-)G-?L4tC=89C2Avc87>=8l~r; zXVY~|*WX4HQ7`RP{DIdt`Tc_W0lmNNzIm{XYGRG?6f}prxRR zh#%S)X}v>2yFxGVBWUb#jDkMx9JA0)zkGkRB;D|}KK-O$^pk!=-w{990avcU7yDeZ zo`X9wKg=1&2L}PsJM}I@?>Ju{T{=BjUXf^!Zgg!G7Yoh(f+Bs}0gs`qV|xSpPda8> zLiz?b^hv8MYr%paxVXxUzHLz5#V+~_4lIMYqJO~2jt{sXg9RU_9}WTz0+%8{p0KYg zwA%8b^_O|-+a7xi z+hy9(XFaE#cp@H%ALiIQuV|I5XEp09Y} zIrNYv#?NROoxIGi&wd122|eTz+4Ipp9>cuHA80G|k@&+`h{PH9i7UU6e4R)ga(sN9_5Nym__ldYOgfdAr1(U0#i|j*f%C?}UK0K-XZ1?Q#9D6aZ+5XXP-|nd&7PJoQuJ5azi(ge>hk6D*{KF2k5%pjveU;ArK$2}W zenR&`_s}1t>^Qmc_Ml_1 zOFQw5zu2)$;+5xlFa02{h&Mg&2nT_nKtS|P{malh&euoRJi?ZqKwh(fAbEiPBMVU` zseiO1)m&f%{YKJpD*#U1ccg=WgMfp;Z$-e80VMXgMqX<>&jpg$bEJbnenh~MQ%HQ| z8hN$tf$AQAdw})cUk7j>sTbqN)k^$zq=UeXM8ML^k@&(j@=ERg`?k0U)b5WbDOiCd z_Utl{zSfOuQFgAAxxKyTP_v0oQZST5-{xF9c4a2gy0 z90Yn_QaksNJa@H{=N#!E;2@BPzu9xFq~8%K5xK&j^pv6Nx`w{@^?_bLPyM@60)K z=FAJ@%f0>5tWr0BUUf8$FEg~2Bfi$K=%Ep@E&C4~J|eW`(CF}}!NDU&1jjn!EB7A} zI&k=~$fyy~1EUoj^2)@-An(oXPG!hSo!$! zgQ8+1LT~gb^J06))G2w*W#da5j-a5RZm!&)o4>l3igJWZ?KG-Jo=)!g`{KyjQO5?) zP#Hp*W9#&H=kIu6wc7vo%afOkU9B!i`}T^_^REP_sb@c0zBKh!r)pnqMeB00+tihy zNBSKNb*U>iz;n4&QrPl`JMO%oR@G~hcE3-Gs*vQTTAVtn9ufZZn9K9Z2d=Z^XKLf4 z$DO8iUrl}egTadwI9(DwcSSZ>-0!d3)Onc)^C_EA`Iia3ved@Nmja_Uq^ZZpJ~yV< zur&4Vg8lu&hhJ8szN_xNvSyBwd6}u+Nt?39XRBsD*DE?Exl|wFKUKYXg8xXT+IZ^d z%yGU>g?xQo(2Moa*RiPaJ12M~dmQ(??&@Qp2YRC4KG{c`(`V#tZs1gVrF^eX`{?;& zoT}^Sv3uhdWGinevyTF$9=z_8cI4EQGV&SsfXHWjaDf&zF1CCdfx8(2-ENo*_5-CJ z|4l#YIX|#V?}YSaO?I?0uzz1UXQ+n{`(>OuVLg)Cu1iy!#UGv(zq?iZtd}EbLPkTU z!VmNmKUzinLvX`4P2a$6%Im8pi^pzrsRPq{uG_pNTWzl~q+Il)DQZNoW?r$2T&n5I zb=M!Lk)ynp4mQ?bNmdP-cDlHB^sndIar6F-<>M}>%=;=|T0JU9#cKY%iTLqKJ666B zztE*}8MDrydcX;81I@2R%Xz7+Y>Y>C{mHcL8tEv-6_ixeQjOu*2 z&VgTY9XWAe$%uk6XUX6>MJwpAD$bsZ)u0`*zr9{z3QA%nPOjFnm#@v zXTh>+C&E9-wckN|Dg7SpPE&fo2CU|-ZXiT`HvgDZRu@H#7T>~nsHD()EH&Z=jf79FV(8ndR8 z!TLane_5}uCGH%Sm^!<~q~QkZ+Pz-p&}UavW39YL|*K;a`Uo_DpTw^zp)V~fPi zp9YUOn4aTO2W6Z`Vj>1~`7~RZLO&P2wZg~w6LjQwR>whUn8p_|}U);PgQ0_^WAw{S5Tg=l#0a z(PW{-`{2krZe4%xW3~Hqws2;>px{rO(N_E;@db7TWj`PvN}wO0k%f&fPf=6 z1C^Ty^L2!k3-h><&IkC0r&9x-M?0 zo`mhIB@IP5Pk0SKZUT3Xndh&~ovYRL#O)B?6YgcmwP>$|Gky~8beHP~X%(4?2_+o! z7yc)NpJS$68#j zpTOrU=@l?lC z^b&mEh};K74$mh8WSpHM7deNtG?Aw~r-aHpJ!M`zw}2b-WM1F_#ZQ4R*AU9!ivOZ~ zM>k6e8eCr=@Mry4C)S60?1Jkl@s7?<+(3GQv;}Fz?efe;`#X|OA${_>q+@d=y;5D$ zDL>2e70+5f$TQb9NvE8YG|?$Z6MZIW+zOH=3Xn9B>q6dhTC}8Dc}8m~>77)OLwYGj z`2HtykwZVm;hBy!(L9+~rp)8-lE$5*rHPowR+%@?bk}8G;7l3`9KjtNcMGn>DGep< zw_8isv7S)YRj2F=#pVfO*yJ+Gm- zp{}!*t~{Uns&sDBsPWaWxzuwKM(2g?$=2&B{FpuNv9NmlS^Vg$>&f*79GK77{aZJx zSjVNB*L}h~9P3ix6%u%SO#JS2^=8#1HFxf>uWm6@5_~G-=e_-c!$qbgYj)@ zBVZ%oJ_ziS^b7t1d&h6;eweR%bpBQ1e^vb;YQ>~3<34I;^t)%rnlj(Fv-}?Z z6$(C-@t2g{XJ_1& zo?B|!gc_$6`#ebE!=a1oCOVd$w7%~9>Dy+}-=57-*vU&C9Y-Gg$)$dlv{K@biosK7 zom5XtJXv4#;X2NBHcIs284WpnC6C>r4>+JVIAV9iA^LMAaT#CTvyR}Rm%*2M*3G6x zi$FoR7gW#su#RZQzQDfD*Z%Q~K0#+XVn6uj9*(l_Z~maU!T!PikV{$#dqEEw$0Yt; zw&|Yls?{jt{x|Ot<<* zz?2GEhj)e;Ysa;5EuI@>uukxUGwWtk8v%DfK+;qG`J|^x!q-OyEg^14cOSxU>Td5g z`~?kD?CWu(xL@2L;mdisW^f(C&8JG(1zyO>5`O(!MH&WoTkL@m1-ZuQ`((?-)GoD7 zjk3aI_Qx|4W^(PzlCRl*?6Ca%37y_Dz!khGBOmuCau<_>os6twY);S<+Hx<=LZ~TyxXJj}mT2O26~MNBwL`AMt#GAHx4~ zZKWQ4p(rgr+MdgAM*zEEzpyXZSL_?zF4$MpbB&}PIq9<4=y{e+w_>Bl&`L~gL8 zo0u2&hTh=9SMJP%>nC3uF>m4lD0T{7$YULd7r+tv+gH}+EDi`@M|oe5CR=Gz>N#In zU+UQx*cW!{yKHfC@{;jJe8Q76_Qi%8>>KurPT5b$BV9_oPr4NO`uWB_1!@Y~=739#Lc`)BCA*UO8Z5(GkV}R4H@Ubo09PzR0nHO?6kFhh_p~wLzaEBI$ zJ!9)&BOvLW*7>A&?5~e*-5!`f6b{@?-fi8ba0|4T6iOI_n+7M21Kc%qHDQ>fIpEj( zy9<6Cf06CGKjRm=(n|k;2kkOnbB*9b+1>GRmvy%1VI%NIBB0wnvJ3K)q?_mmzfhC& zYG*U=!?UN4vHsXqucyqphCDq?r;L24-jDuJO)tA>BXCC$AZ}q_at?F-#Ll37<(US4 zsK@g){TnY`VsP!mF6hUzj!e?j;zX1G!LnS=efT%ZKQSPVVm#V){i$ahD0)E|hxRYN zOi8Vs`na+Ev|s!~0j3q7peMNKl=-2T$g_!!z)b|$9}O1`E|b=FnzgUkFIC6*e-gj$ zH3L5CdtPjpRH<>KRWERw*cs;&c1Az=IDbBUFgT_1Td|g17Gy7^kFXEyjdafgyC!A6 zX1-!DF5^In<8?}V*52vI=MCs&_|&Wa_!m8VZ+Xteu9yd9?2mI@!aJOxfR~yHg3c3cFnkb+%)b~CTfLUxar&Qml%AKr-lWzT+_q44uV%n!^iCw z+@6q)K;aRv&Jo>z7UWk8s+XZnD116)9Gw*-ezMV?J<1d1a9Cg-rGoGE2!DUHvn*ZB%My%@~GFTCZBO{ zN2QxhYy|Eo0;Y@(MH$Lat6sAU>Y>(GNK^>5>gkUi(XLaS4}Pf5ua%`}6B_{=fj<%f zNl&%QCp~33xi!@BVZ%owg{LqGbj_K3@vOu z{R*PUv#E`MJ0oCflcNl!nKHDXdfK7z>6CVz>U?(DM!+2qFl7=@=0q8KyL#H8^!u%p Tem1ocuo1{b;6{2Uy5;`?ZO0>^ From b4b2646c0527d239e2de0b512215c20c7be45162 Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Thu, 30 Jul 2026 23:33:09 -0400 Subject: [PATCH 03/14] ci: build the Rust extension on all runners and test on windows-latest --- .github/workflows/coverage.yml | 1 + .github/workflows/docs.yml | 1 + .github/workflows/lint.yml | 7 +++ .github/workflows/release.yml | 56 +++++++++++++++---- .../workflows/test_development_versions.yml | 1 + .github/workflows/test_latest_versions.yml | 3 + .github/workflows/test_minimum_versions.yml | 1 + 7 files changed, 59 insertions(+), 11 deletions(-) diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index 68a02c2..76b9c6d 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -26,6 +26,7 @@ jobs: uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} + - uses: dtolnay/rust-toolchain@stable - name: Install tox run: | python -m pip install --upgrade pip diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 9bf20b0..988ffa3 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -25,6 +25,7 @@ jobs: - uses: actions/setup-python@v6 with: python-version: '3.10' + - uses: dtolnay/rust-toolchain@stable - name: Install dependencies run: | python -m pip install --upgrade pip diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index 0e3a94c..2395766 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -27,6 +27,13 @@ jobs: run: | python -m pip install --upgrade pip pip install tox + - uses: dtolnay/rust-toolchain@stable + with: + components: rustfmt, clippy + - name: Rust Format + run: cargo fmt --all -- --check + - name: Clippy + run: cargo clippy --all-targets -- -D warnings - name: Run lint check run: | tox -e lint diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 7698cae..80564ab 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -24,26 +24,60 @@ jobs: token: ${{ secrets.GITHUB_TOKEN }} generateReleaseNotes: "true" + wheels: + name: wheels (${{ matrix.os }}) + runs-on: ${{ matrix.os }} + needs: github + strategy: + matrix: + os: [ubuntu-latest, macos-latest, windows-latest] + steps: + - name: Checkout tag + uses: actions/checkout@v6 + with: + ref: ${{ github.ref_name }} + - uses: actions/setup-python@v6 + - uses: dtolnay/rust-toolchain@stable + - name: Build wheels + uses: pypa/cibuildwheel@v3.4.0 + - uses: actions/upload-artifact@v6 + with: + name: wheels-${{ matrix.os }} + path: ./wheelhouse/*.whl + + sdist: + name: sdist + runs-on: ubuntu-latest + needs: github + steps: + - name: Checkout tag + uses: actions/checkout@v6 + with: + ref: ${{ github.ref_name }} + - uses: actions/setup-python@v6 + - name: Build sdist + run: | + python -m pip install --upgrade pip "maturin>=1.9,<2.0" build + python -m build --sdist + - uses: actions/upload-artifact@v6 + with: + name: sdist + path: ./dist/*.tar.gz + pypi: name: pypi runs-on: ubuntu-latest - needs: github + needs: [wheels, sdist] environment: name: pypi url: https://pypi.org/p/qiskit-addon-slc permissions: id-token: write steps: - - name: Checkout tag - uses: actions/checkout@v6 + - name: Download all artifacts + uses: actions/download-artifact@v6 with: - ref: ${{ github.ref_name }} - - name: Install `build` tool - run: | - python -m pip install --upgrade pip - pip install build - - name: Build distribution - run: | - python -m build + path: dist + merge-multiple: true - name: Publish release to PyPI uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/.github/workflows/test_development_versions.yml b/.github/workflows/test_development_versions.yml index a4e583c..fd3aeb4 100644 --- a/.github/workflows/test_development_versions.yml +++ b/.github/workflows/test_development_versions.yml @@ -28,6 +28,7 @@ jobs: uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} + - uses: dtolnay/rust-toolchain@stable - name: Upgrade pip run: | python -m pip install --upgrade pip diff --git a/.github/workflows/test_latest_versions.yml b/.github/workflows/test_latest_versions.yml index 14fd44b..fae4445 100644 --- a/.github/workflows/test_latest_versions.yml +++ b/.github/workflows/test_latest_versions.yml @@ -25,12 +25,15 @@ jobs: include: - os: macos-latest python-version: "3.13" + - os: windows-latest + python-version: "3.13" steps: - uses: actions/checkout@v6 - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} + - uses: dtolnay/rust-toolchain@stable - name: Install dependencies run: | python -m pip install --upgrade pip diff --git a/.github/workflows/test_minimum_versions.yml b/.github/workflows/test_minimum_versions.yml index 87a9bca..a3a7637 100644 --- a/.github/workflows/test_minimum_versions.yml +++ b/.github/workflows/test_minimum_versions.yml @@ -26,6 +26,7 @@ jobs: uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} + - uses: dtolnay/rust-toolchain@stable - name: Install dependencies (minimum versions) shell: bash run: | From b0eb7aed205ef9e18eaa4067d3f4ffc01af818dc Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Thu, 30 Jul 2026 23:38:20 -0400 Subject: [PATCH 04/14] docs: announce Windows support and the Rust backend --- README.md | 5 ++--- docs/conf.py | 1 - .../rust-davidson-windows-0976dac383c66486.yaml | 17 +++++++++++++++++ 3 files changed, 19 insertions(+), 4 deletions(-) create mode 100644 releasenotes/notes/rust-davidson-windows-0976dac383c66486.yaml diff --git a/README.md b/README.md index 5ce54ef..ecfc59e 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@
[![Release](https://img.shields.io/pypi/v/qiskit-addon-slc.svg?label=Release)](https://github.com/Qiskit/qiskit-addon-slc/releases) - ![Platform](https://img.shields.io/badge/%F0%9F%92%BB%20Platform-Linux%20%7C%20macOS-informational) + ![Platform](https://img.shields.io/badge/%F0%9F%92%BB%20Platform-Linux%20%7C%20macOS%20%7C%20Windows-informational) [![Python](https://img.shields.io/pypi/pyversions/qiskit-addon-slc?label=Python&logo=python)](https://www.python.org/) [![Qiskit](https://img.shields.io/badge/Qiskit%20-%20%3E%3D2.2%20-%20%236133BD?logo=Qiskit)](https://github.com/Qiskit/qiskit)
@@ -85,17 +85,16 @@ Shaded lightcones are calculated and used in 5 steps: - Parallel asynchronous bound computation - [Rust-accelerated Pauli propagation](https://github.com/Qiskit/pauli-prop) +- Rust-accelerated eigenvalue computation for computing forward bounds - Permits ahead-of-time bound computation (i.e. prior to the actual noise learning) #### Known issues -- Windows not supported - `InjectNoise(site="before")` not supported - Does not support fine-grained bound merging #### Future work -- Rust-accelerated eigenvalue computation for computing forward bounds - Additional guides coming soon ---------------------------------------------------------------------------------------------------- diff --git a/docs/conf.py b/docs/conf.py index 9ab227a..bfd012f 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -103,7 +103,6 @@ "python": ("https://docs.python.org/3", None), "numpy": ("https://numpy.org/doc/stable/", None), "scipy": ("https://docs.scipy.org/doc/scipy/", None), - "pyscf": ("https://pyscf.org/", None), "qiskit": ("https://quantum.cloud.ibm.com/docs/api/qiskit/", None), "qiskit_ibm_runtime": ("https://qiskit.github.io/qiskit-ibm-runtime/", None), "qiskit_addon_utils": ("https://qiskit.github.io/qiskit-addon-utils/", None), diff --git a/releasenotes/notes/rust-davidson-windows-0976dac383c66486.yaml b/releasenotes/notes/rust-davidson-windows-0976dac383c66486.yaml new file mode 100644 index 0000000..585e3ed --- /dev/null +++ b/releasenotes/notes/rust-davidson-windows-0976dac383c66486.yaml @@ -0,0 +1,17 @@ +--- +features: + - | + Added support for Windows. The package now installs and runs on Windows in + addition to Linux and macOS, and the test suite is exercised on + ``windows-latest`` in CI. + - | + The Davidson eigensolver used internally by + :func:`~qiskit_addon_slc.utils.davidson.get_extremal_eigenvalue` is now + provided by a compiled Rust extension (built on ``nalgebra``). It requires no + BLAS/LAPACK, so the package ships self-contained binary wheels. +upgrade: + - | + ``pyscf`` is no longer a dependency. The Davidson eigensolver behind + :func:`~qiskit_addon_slc.utils.davidson.get_extremal_eigenvalue` was + previously provided by ``pyscf``, which had no Windows wheel and could not be + built from source on Windows; it is now the compiled Rust extension. From f6487e83b4a3e15f9805918cfe725cba31e45e7c Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 27 Aug 2026 10:45:45 -0500 Subject: [PATCH 05/14] refactor: some opinionated source restructuring This change is to follow other Qiskit addon source layout patterns more closely. --- Cargo.lock | 2 +- Cargo.toml | 5 ++--- pyproject.toml | 7 ++++--- qiskit_addon_slc/{_davidson.pyi => _accelerate.pyi} | 0 qiskit_addon_slc/utils/davidson.py | 4 ++-- rust/davidson.rs => src/lib.rs | 2 +- 6 files changed, 10 insertions(+), 10 deletions(-) rename qiskit_addon_slc/{_davidson.pyi => _accelerate.pyi} (100%) rename rust/davidson.rs => src/lib.rs (98%) diff --git a/Cargo.lock b/Cargo.lock index 7e43789..3fa30bf 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -267,7 +267,7 @@ dependencies = [ ] [[package]] -name = "qiskit-addon-slc" +name = "qiskit_addon_slc_accelerate" version = "0.1.0" dependencies = [ "nalgebra", diff --git a/Cargo.toml b/Cargo.toml index 087d18b..eadbfc6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,5 +1,5 @@ [package] -name = "qiskit-addon-slc" +name = "qiskit_addon_slc_accelerate" description = "Rust accelerator for qiskit-addon-slc" version = "0.1.0" edition = "2024" @@ -7,8 +7,7 @@ rust-version = "1.85" license = "Apache-2.0" [lib] -name = "_davidson" -path = "rust/davidson.rs" +name = "_accelerate" crate-type = ["cdylib"] [dependencies] diff --git a/pyproject.toml b/pyproject.toml index dae9110..d730b5d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -93,11 +93,12 @@ fail_under = 65 show_missing = true [tool.maturin] -features = ["pyo3/extension-module"] -locked = true bindings = "pyo3" python-source = "." -module-name = "qiskit_addon_slc._davidson" +features = ["pyo3/extension-module"] +locked = true +profile = "release" +module-name = "qiskit_addon_slc._accelerate" [tool.mypy] python_version = "3.10" diff --git a/qiskit_addon_slc/_davidson.pyi b/qiskit_addon_slc/_accelerate.pyi similarity index 100% rename from qiskit_addon_slc/_davidson.pyi rename to qiskit_addon_slc/_accelerate.pyi diff --git a/qiskit_addon_slc/utils/davidson.py b/qiskit_addon_slc/utils/davidson.py index 8c4a558..dc9b1a9 100644 --- a/qiskit_addon_slc/utils/davidson.py +++ b/qiskit_addon_slc/utils/davidson.py @@ -21,7 +21,7 @@ import numpy as np from qiskit.quantum_info import SparsePauliOp -from .. import _davidson +from .. import _accelerate def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: @@ -61,7 +61,7 @@ def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: diag = spmat.diagonal().astype(np.complex128) seed = _random_initial_guess((dim,)).astype(np.complex128) - return _davidson.davidson_smallest( + return _accelerate.davidson_smallest( spmat.indptr.astype(np.int64), spmat.indices.astype(np.int64), np.ascontiguousarray(data.real), diff --git a/rust/davidson.rs b/src/lib.rs similarity index 98% rename from rust/davidson.rs rename to src/lib.rs index a7575e1..afb7cc6 100644 --- a/rust/davidson.rs +++ b/src/lib.rs @@ -170,7 +170,7 @@ fn davidson_smallest( } #[pymodule] -fn _davidson(m: &Bound<'_, PyModule>) -> PyResult<()> { +fn _accelerate(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(davidson_smallest))?; Ok(()) } From c0938c2bd04dce1b844e5f5a2f0656c8f8d95bc4 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 27 Aug 2026 10:53:42 -0500 Subject: [PATCH 06/14] ci: bulid package before attempting to run tests --- .github/workflows/coverage.yml | 3 +++ .github/workflows/test_development_versions.yml | 3 +++ .github/workflows/test_latest_versions.yml | 3 +++ .github/workflows/test_minimum_versions.yml | 3 +++ 4 files changed, 12 insertions(+) diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index 6b93e55..c3885d7 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -31,6 +31,9 @@ jobs: run: | python -m pip install --upgrade pip pip install tox coverage + - name: Build package + run: | + pip install -e . - name: Run coverage run: | tox -e coverage diff --git a/.github/workflows/test_development_versions.yml b/.github/workflows/test_development_versions.yml index 3d3bfea..0219f60 100644 --- a/.github/workflows/test_development_versions.yml +++ b/.github/workflows/test_development_versions.yml @@ -69,6 +69,9 @@ jobs: "qiskit-ibm-runtime @ file:$(echo qiskit-ibm-runtime/dist/*.whl)" "samplomatic @ file:$(echo samplomatic/dist/*.whl)" "qiskit-addon-utils @ file:$(echo qiskit-addon-utils/dist/*.whl)" + - name: Build package + run: | + pip install -e . - name: Test using tox environment run: | tox -e py --parallel --parallel-no-spinner diff --git a/.github/workflows/test_latest_versions.yml b/.github/workflows/test_latest_versions.yml index 2e75c76..24dac05 100644 --- a/.github/workflows/test_latest_versions.yml +++ b/.github/workflows/test_latest_versions.yml @@ -38,6 +38,9 @@ jobs: run: | python -m pip install --upgrade pip pip install tox + - name: Build package + run: | + pip install -e . - name: Test using tox environment shell: bash run: | diff --git a/.github/workflows/test_minimum_versions.yml b/.github/workflows/test_minimum_versions.yml index 3d04868..887dddd 100644 --- a/.github/workflows/test_minimum_versions.yml +++ b/.github/workflows/test_minimum_versions.yml @@ -34,6 +34,9 @@ jobs: python -m pip install extremal-python-dependencies==0.1.0 pip install "tox==$(extremal-python-dependencies get-tox-minversion)" extremal-python-dependencies pin-dependencies-to-minimum --inplace + - name: Build package + run: | + pip install -e . - name: Test using tox environment run: | tox -e py --parallel --parallel-no-spinner From 9d411855d436d74d9083272b3feedae1f4ac7570 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 27 Aug 2026 11:04:59 -0500 Subject: [PATCH 07/14] ci: fix package build location in workflows --- .github/workflows/docs.yml | 3 +++ .github/workflows/test_development_versions.yml | 6 +++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 3fc88e6..f0563ad 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -32,6 +32,9 @@ jobs: pip install tox sudo apt-get update sudo apt-get install -y pandoc + - name: Build package + run: | + pip install -e . - name: Tell reno to name the upcoming release after the branch we are on shell: bash run: | diff --git a/.github/workflows/test_development_versions.yml b/.github/workflows/test_development_versions.yml index 0219f60..3b6986c 100644 --- a/.github/workflows/test_development_versions.yml +++ b/.github/workflows/test_development_versions.yml @@ -35,6 +35,9 @@ jobs: - name: Install tools from pypi run: | python -m pip install tox build extremal-python-dependencies==0.1.0 + - name: Build package + run: | + pip install -e . - name: Build Qiskit SDK development wheel run: | git clone https://github.com/Qiskit/qiskit @@ -69,9 +72,6 @@ jobs: "qiskit-ibm-runtime @ file:$(echo qiskit-ibm-runtime/dist/*.whl)" "samplomatic @ file:$(echo samplomatic/dist/*.whl)" "qiskit-addon-utils @ file:$(echo qiskit-addon-utils/dist/*.whl)" - - name: Build package - run: | - pip install -e . - name: Test using tox environment run: | tox -e py --parallel --parallel-no-spinner From 44ce31fbf27f42922ff57dd694fcc5a987fcb95a Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 27 Aug 2026 11:12:30 -0500 Subject: [PATCH 08/14] ci: update release workflow --- .github/workflows/release.yml | 90 +++++++++++++++++++---------------- 1 file changed, 50 insertions(+), 40 deletions(-) diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index dc097b8..5569950 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -8,75 +8,85 @@ on: jobs: - github: - name: github + build_sdist: + name: Build source distribution runs-on: ubuntu-latest steps: - - name: Checkout tag - uses: actions/checkout@v7 - with: - ref: ${{ github.ref_name }} - - name: Publish release - uses: ghalactic/github-release-from-tag@v6 - if: github.ref_type == 'tag' + - uses: actions/checkout@v7 + - uses: actions/setup-python@v6 + - name: Build distribution + run: | + python -m pip install --upgrade pip "maturin>=1.9,<2.0" build + python -m build --sdist + - name: Store the distribution packages + uses: actions/upload-artifact@v7 with: - prerelease: false - token: ${{ secrets.GITHUB_TOKEN }} - generateReleaseNotes: "true" + name: pkg-sdist + path: dist/*.tar.gz - wheels: - name: wheels (${{ matrix.os }}) + build_wheels: + name: Build wheel on ${{ matrix.os }} runs-on: ${{ matrix.os }} - needs: github strategy: matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, ubuntu-24.04-arm, macos-14, macos-15-intel, windows-latest] + steps: - - name: Checkout tag - uses: actions/checkout@v6 - with: - ref: ${{ github.ref_name }} + - uses: actions/checkout@v7 - uses: actions/setup-python@v6 - uses: dtolnay/rust-toolchain@stable - - name: Build wheels + - name: Build wheel uses: pypa/cibuildwheel@v3.4.0 - - uses: actions/upload-artifact@v6 + - name: Store the wheel + uses: actions/upload-artifact@v7 with: - name: wheels-${{ matrix.os }} + name: pkg-wheel-${{ matrix.os }} path: ./wheelhouse/*.whl - sdist: - name: sdist + github: + name: Publish to GitHub + if: startsWith(github.ref, 'refs/tags/') runs-on: ubuntu-latest - needs: github + needs: + - build_sdist + - build_wheels steps: - name: Checkout tag - uses: actions/checkout@v6 + uses: actions/checkout@v7 + - name: Download artifacts + uses: actions/download-artifact@v8 with: - ref: ${{ github.ref_name }} - - uses: actions/setup-python@v6 - - name: Build sdist - run: | - python -m pip install --upgrade pip "maturin>=1.9,<2.0" build - python -m build --sdist - - uses: actions/upload-artifact@v6 + pattern: pkg-* + path: dist + merge-multiple: true + - name: Publish release + uses: ghalactic/github-release-from-tag@v6 with: - name: sdist - path: ./dist/*.tar.gz + prerelease: false + generateReleaseNotes: "true" + assets: | + - path: dist/* pypi: - name: pypi + name: Publish to pypi + if: startsWith(github.ref, 'refs/tags/') runs-on: ubuntu-latest - needs: [wheels, sdist] + needs: + - build_sdist + - build_wheels + - github environment: name: pypi url: https://pypi.org/p/qiskit-addon-slc permissions: id-token: write steps: - - name: Download all artifacts - uses: actions/download-artifact@v6 + - name: Checkout tag + uses: actions/checkout@v7 + - name: Download artifacts + uses: actions/download-artifact@v8 with: + pattern: pkg-* path: dist merge-multiple: true - name: Publish release to PyPI From 94202678f36aa3768c9c6a42cd77b3e41121e2c6 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 27 Aug 2026 11:20:04 -0500 Subject: [PATCH 09/14] refactor: make kwargs explicit --- qiskit_addon_slc/utils/davidson.py | 45 +++++++++++++++++------------- 1 file changed, 26 insertions(+), 19 deletions(-) diff --git a/qiskit_addon_slc/utils/davidson.py b/qiskit_addon_slc/utils/davidson.py index dc9b1a9..b37c62f 100644 --- a/qiskit_addon_slc/utils/davidson.py +++ b/qiskit_addon_slc/utils/davidson.py @@ -16,6 +16,7 @@ """A basic Davidson solver.""" +import warnings from typing import cast import numpy as np @@ -24,7 +25,15 @@ from .. import _accelerate -def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: +def get_extremal_eigenvalue( + spo: SparsePauliOp, + *, + tol: float = 1e-6, + max_cycle: int = 500, + max_space: int = 12, + lindep: float = 1e-11, + **kwargs, +) -> tuple[bool, float]: """Finds the extremal eigenvalue of the provided operator. The operator is converted to a sparse matrix, whose smallest eigenvalue is then computed by the @@ -35,25 +44,23 @@ def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: Args: spo: the operator whose minimal eigenvalue to find. - kwargs: additional keyword arguments for the Davidson algorithm. When not specified - otherwise, the following defaults will be used: - - * `tol`: 1e-6 - * `max_cycle`: 500 - * `max_space`: 12 - * `lindep`: 1e-11 + tol: TODO. + max_cycle: TODO. + max_space: TODO. + lindep: TODO. + kwargs: **ignored!** Any additional keyword arguments are parsed for backwards compatibility + but do not have any effect at runtime and, thus, are being ignored! Returns: A pair indicating whether the Davidson algorithm has converged and the obtained minimal eigenvalue. """ - default_kwargs = { - "tol": 1e-6, - "max_cycle": 500, - "max_space": 12, - "lindep": 1e-11, - } - default_kwargs.update(kwargs) + if len(kwargs) > 0: + warnings.warn( + f"These keyword arguments do not have any effect and are ignored: {kwargs}", + category=UserWarning, + stacklevel=2, + ) spmat = spo.to_matrix(sparse=True, force_serial=True).tocsr() dim = spmat.shape[0] @@ -71,10 +78,10 @@ def get_extremal_eigenvalue(spo: SparsePauliOp, **kwargs) -> tuple[bool, float]: np.ascontiguousarray(seed.real), np.ascontiguousarray(seed.imag), dim, - float(default_kwargs["tol"]), - int(default_kwargs["max_cycle"]), - int(default_kwargs["max_space"]), - float(default_kwargs["lindep"]), + float(tol), + int(max_cycle), + int(max_space), + float(lindep), ) From dd67f97ab92dae3452bbba615f29b5b89c254816 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 27 Aug 2026 11:43:39 -0500 Subject: [PATCH 10/14] fix: module name in stub-file docstring --- qiskit_addon_slc/_accelerate.pyi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit_addon_slc/_accelerate.pyi b/qiskit_addon_slc/_accelerate.pyi index 3279c6c..965bbcf 100644 --- a/qiskit_addon_slc/_accelerate.pyi +++ b/qiskit_addon_slc/_accelerate.pyi @@ -10,7 +10,7 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Type stubs for the compiled Rust extension ``qiskit_addon_slc._davidson``.""" +"""Type stubs for the compiled Rust extension ``qiskit_addon_slc._accelerate``.""" import numpy as np import numpy.typing as npt From 2a5fa864d808ef496e4d3dfd6132dd2808374b02 Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Fri, 28 Aug 2026 12:46:25 -0400 Subject: [PATCH 11/14] fix: converge the Davidson solver on a residual and value criterion --- src/lib.rs | 97 ++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 79 insertions(+), 18 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index afb7cc6..5c66cb2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -12,13 +12,26 @@ //! Davidson eigensolver for the algebraically smallest eigenvalue of a Hermitian sparse matrix. //! -//! The dense linear algebra uses `nalgebra`. The sparse operator is held as a plain CSR triple and -//! its matrix-vector product is applied directly, because `nalgebra-sparse` does not implement +//! The dense linear algebra uses `nalgebra`, while the sparse operator is held as a plain CSR triple +//! whose matrix-vector product is applied directly, because `nalgebra-sparse` does not implement //! sparse-times-dense multiplication for complex scalars. +//! +//! The stopping rule mirrors `pyscf.lib.davidson1`. A cycle converges when the Ritz value has settled +//! (`|Δθ| < tol`) and the residual `A·x - θ·x` of the current Ritz pair `(θ, x)` has norm below +//! `sqrt(tol)`. If the correction vanishes after orthogonalization the subspace is exhausted, and +//! that same residual test then decides convergence, as in `davidson1` (`conv = dx_norm < toloose`). +//! +//! The Jacobi preconditioner divides each correction entry by the shift `diag[i] - θ`. A shift whose +//! magnitude falls below `floor = 1e-12 * max(diag_scale, 1)`, where `diag_scale` is the largest +//! `|diag[i]|`, is clamped up to `floor`, keeping its sign so it is not pushed the wrong way. This +//! floor is relative to the operator scale rather than to `tol`, because `tol` is a residual +//! threshold and clamping to it would corrupt operators with a near-zero diagonal. When the diagonal +//! is entirely zero the preconditioner has no useful shift to apply, so it is skipped altogether. use nalgebra::{DMatrix, DVector}; use num_complex::Complex64 as C64; use numpy::PyReadonlyArray1; +use pyo3::exceptions::PyValueError; use pyo3::prelude::*; /// A Hermitian operator held in compressed-sparse-row form. @@ -62,6 +75,7 @@ fn columns_to_matrix(cols: &[DVector]) -> DMatrix { DMatrix::from_columns(&cols.iter().map(|c| c.column(0)).collect::>()) } +/// Iterates for the algebraically smallest eigenvalue of `op`, returning `(converged, eigenvalue)`. fn davidson( op: &CsrOp, diag: &DVector, @@ -73,6 +87,11 @@ fn davidson( ) -> (bool, f64) { let dim = op.dim; + let diag_scale = diag.iter().map(|z| z.norm()).fold(0.0_f64, f64::max); + let precondition = diag_scale > 0.0; + let floor = 1e-12 * diag_scale.max(1.0); + let residual_tol = tol.sqrt(); + // Subspace basis vectors `s` and their images `A @ s`, grown one vector per cycle. let mut images: Vec> = vec![op.apply(&seed)]; let mut s: Vec> = vec![seed]; @@ -95,20 +114,25 @@ fn davidson( let ritz_image = &images_mat * &y; let residual = &ritz_image - ritz.scale(theta); - if (eigval - prev).abs() < tol || residual.norm() < tol { + let residual_norm = residual.norm(); + let de = (theta - prev).abs(); + prev = theta; + if residual_norm < residual_tol && de < tol { converged = true; break; } - prev = eigval; - // Diagonal (Jacobi) preconditioner, clamping near-zero shifts to `tol`. + // Apply the preconditioner, flooring the shift while keeping its sign. let mut correction = residual; - for i in 0..dim { - let mut d = diag[i] - C64::new(theta, 0.0); - if d.norm() < tol { - d = C64::new(tol, 0.0); + if precondition { + for i in 0..dim { + let mut d = diag[i] - C64::new(theta, 0.0); + if d.norm() < floor { + let sign = if d.re < 0.0 { -floor } else { floor }; + d = C64::new(sign, 0.0); + } + correction[i] /= d; } - correction[i] /= d; } // Collapse the subspace to the current best estimate before it exceeds `max_space`. @@ -117,12 +141,13 @@ fn davidson( images = vec![ritz_image]; } - // Orthonormalize the correction against the subspace (modified Gram-Schmidt). + // Classical Gram-Schmidt with one re-orthogonalization pass (numerically comparable to MGS). let s_mat = columns_to_matrix(&s); correction -= &s_mat * (s_mat.adjoint() * &correction); + correction -= &s_mat * (s_mat.adjoint() * &correction); let cnorm = correction.norm(); if cnorm < lindep { - converged = true; + converged = residual_norm < residual_tol; break; } correction.unscale_mut(cnorm); @@ -134,6 +159,9 @@ fn davidson( (converged, eigval) } +/// Python entry point: builds a [`CsrOp`] from the split real/imaginary arrays and runs [`davidson`]. +/// +/// Raises `ValueError` if the CSR arrays, `dim`, `diag`, or `seed` are inconsistent. #[pyfunction] #[allow(clippy::too_many_arguments)] fn davidson_smallest( @@ -155,15 +183,48 @@ fn davidson_smallest( re.iter().zip(im).map(|(a, b)| C64::new(*a, *b)).collect() } + let indptr = indptr.as_slice()?; + let indices = indices.as_slice()?; + let data = to_complex(data_re.as_slice()?, data_im.as_slice()?); + let diag = to_complex(diag_re.as_slice()?, diag_im.as_slice()?); + let seed = to_complex(seed_re.as_slice()?, seed_im.as_slice()?); + + if dim == 0 { + return Err(PyValueError::new_err("`dim` must be positive")); + } + if max_space < 2 { + return Err(PyValueError::new_err("`max_space` must be at least 2")); + } + if indptr.len() != dim + 1 { + return Err(PyValueError::new_err("`indptr` must have length `dim + 1`")); + } + if !indptr.windows(2).all(|w| w[0] <= w[1]) { + return Err(PyValueError::new_err("`indptr` must be non-decreasing")); + } + if indptr[dim] as usize != indices.len() || indices.len() != data.len() { + return Err(PyValueError::new_err( + "`indptr[dim]` must equal `indices.len()` and `data.len()`", + )); + } + if indices.iter().any(|&j| j < 0 || j as usize >= dim) { + return Err(PyValueError::new_err( + "`indices` entries must be in `[0, dim)`", + )); + } + if diag.len() != dim || seed.len() != dim { + return Err(PyValueError::new_err( + "`diag` and `seed` must have length `dim`", + )); + } + let op = CsrOp { - indptr: indptr.as_slice()?.to_vec(), - indices: indices.as_slice()?.to_vec(), - data: to_complex(data_re.as_slice()?, data_im.as_slice()?), + indptr: indptr.to_vec(), + indices: indices.to_vec(), + data, dim, }; - - let diag = DVector::from_vec(to_complex(diag_re.as_slice()?, diag_im.as_slice()?)); - let seed = DVector::from_vec(to_complex(seed_re.as_slice()?, seed_im.as_slice()?)); + let diag = DVector::from_vec(diag); + let seed = DVector::from_vec(seed); let (conv, ev) = davidson(&op, &diag, seed, tol, max_cycle, max_space, lindep); Ok((conv, ev)) From 7656f5f0546d428dd237533655133f4e81b636d8 Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Fri, 28 Aug 2026 12:48:14 -0400 Subject: [PATCH 12/14] fix: use a fixed-seed initial guess for reproducible bounds --- qiskit_addon_slc/utils/davidson.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/qiskit_addon_slc/utils/davidson.py b/qiskit_addon_slc/utils/davidson.py index b37c62f..788f8de 100644 --- a/qiskit_addon_slc/utils/davidson.py +++ b/qiskit_addon_slc/utils/davidson.py @@ -66,7 +66,7 @@ def get_extremal_eigenvalue( dim = spmat.shape[0] data = spmat.data.astype(np.complex128) diag = spmat.diagonal().astype(np.complex128) - seed = _random_initial_guess((dim,)).astype(np.complex128) + seed = _initial_guess((dim,)).astype(np.complex128) return _accelerate.davidson_smallest( spmat.indptr.astype(np.int64), @@ -85,20 +85,26 @@ def get_extremal_eigenvalue( ) -def _random_initial_guess(shape: tuple[int, ...]) -> np.ndarray: - """Produces a random array of the requested shape. +def _initial_guess(shape: tuple[int, ...]) -> np.ndarray: + """Produces a deterministic normalized starting vector of the requested shape. + + A fixed-seed local generator is used so that the + Davidson iteration is reproducible: the same operator always yields the same result, independent + of any surrounding random state. A pseudo-random (rather than constant) vector is used to avoid + initial guesses that are accidentally orthogonal to the target eigenvector. Args: shape: the requested shape. Returns: - An array of random complex values with their real and imaginary parts lying in the interval + A unit-norm array of complex values with their real and imaginary parts lying in the interval ``[0, 1)``. """ - norm = 0.0 + rng = np.random.default_rng(0) + norm = 0.0 while norm == 0: - x = np.random.rand(shape[0]) + 1.0j * np.random.rand(shape[0]) + x = rng.random(shape[0]) + 1.0j * rng.random(shape[0]) norm = cast(float, np.linalg.norm(x)) return x / norm From c9207c02489f245dba907da294b5ccb575ced6a5 Mon Sep 17 00:00:00 2001 From: Mostafa Atallah Date: Fri, 28 Aug 2026 12:48:40 -0400 Subject: [PATCH 13/14] test: regenerate forward-bound fixtures for the fixed solver --- tests/expected_fwd_bounds.pickle | Bin 37490 -> 37490 bytes tests/expected_fwd_tightened_bounds.pickle | Bin 37490 -> 37490 bytes tests/expected_merged_bounds.pickle | Bin 37490 -> 37490 bytes 3 files changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/expected_fwd_bounds.pickle b/tests/expected_fwd_bounds.pickle index b5991b11c322e035382f77d6e2a272ecd19d2ae7..2fd0c084f6a8e9338b66cb6be4eebebc4b38d964 100644 GIT binary patch literal 37490 zcmeI4cX$<58pcCtp;rYdD~PU$f;0l#d>)E=I!Kd7**`<8^XtXV@sx#?3Iu- zsQ2+}!XLfSJK?Fo>EbEHJl>d?n7jOWe@^@r6_()bF`@N)RRU{r{|+MH=wHa-3HKi- z5U_Evk6oAo2!IRH#l=OyMWE;);K)=YmUE0e-*K)BCDG?f7lFczfMXfQF_P09BSVdI z9?3jga^9BqJonf|pfDid$Ot4hag03MIOmaEcdnA_TxpO?SlyRhni;Tc- zB!)X7aN@ixT?AYNTm=4K1RNy zqTUDl0L~+IW7cu5lJ&aMMc_gr;OOQ^*1|FJY~#W6R$PP{55|)wIH4r^Y#D04#$Dke z;3DulARv0`x--yIC%#q6qo>@jk50qgdHdlbX;XXnO_{lUJ@5UJZcdv~GzOdn!lCS` zI5-t_3{HY0)cnGV)@k+z%HF~tO_zK8Q$N2S|2b`bVFkvu%SFINAS_r%?EFdIJlAny z=^(k@mguvkyB973`3M1Lk@|{qR$s9&M87fL>c@Dzq(7-7J@h5SNw3%E&)qfBXPDPp z$APl!H|XIy=Dnl)QNMfYy7IbD4;^!#P!c;?KlWiaGU8Mm^3hhgODTW|IQW_6S7#lk zZs5=HOZS&s_1Ql>Na_uKizF^MkH4dbapDwx73KG=mvQtFkLa-_dbu8n-N+*T94FXZ zQP>a=z0>jx^p5-W(dnB_T0t_yPBTR{T|A<;ri#)s<)Ev^Nw0R=6^;(`N$W^(F6fdr z(Nk&cr9)dR@?bCigT1tk;AP`O*+dp}T-^1$2)GD@Mu0qF*B#;?i3X0zH;hXLX(b+C zLyeP{IL|t0&pD4|o-Mh~mYU}Z7l8|k0Q3xxm&6bBkKNDoQ9r|c`aR-_ctBcl=*U>7 z@~O0q1FiI~MIU_;JHOS7Ma$fgyg^NI;MnLpbi!;zuMfvL^%5b1PNjGUz z9(!`mh(*D#M1EBoBh8&8Yxb&t=ZD+OZ^Gr+ z#;X!H{(QiY_DQp4G=Hb&lWuMO1qCq5O5m;KNqn4JDzutd>LU;K;UB!+lzk56koiGe z)wbIDH=OE+1qCK|oh|}}0s-1&;+=RTPN^66xIN!)U&#;gLlXDC8tXUHDyEw0)BR?T z^5nj8!$w#B&2OS6&t14L?`y>f-L6OV=I<{SAJewtCWUx~!qwwXICz1^}z z(z0NDjCpEmqq1dvW|PkISXgC7*jr<&{$!fe>*#Dw@|%u|gT&eO7I^3GHtV{!ydt8m z-^@^)jEvkdJmrfl6Ic4?p&R!24Kz_FT{l(ps!4CM^l)8M-Dj=z4w9Y#EorK7jj|ZvSxJzMUT~V87X4@-F!rKjJxs`NSQ5 zK-}%jedzv=pMKZVF>2vHXYD^dvGT+3oM7h-CxT-@FSkk2VUULyuiBC4Y-0kah4x^ z#$AJpfasm(XP|f7ua8c5kKGn6L0ETC`QNYhnzThape!?gi_!||<@7t{kCBWJ>%=bT z2+|tu(zR-HpnZm`ZK92$9E?+2#$yb$Pe-*iJ9U5b=MPj)u7Wm<{p1YzP$pPEWer+~ zGECX4q4LK0IF+Nr!h!O$OUG|+-1Xd|R;eD!4EJ}Xi-3zj;X&X>)g9VR;+?!{kM9UE zpLUzJnLGtO1g%8g3b6GjfBXuqRAKd%vqvoRo0)ptrR^r~Cy1tEoV?Cg6FlbHv=Zx~ zo#t_Z_3+jfNnWLnl8^1W#WAw*ct*jrgnG$-VFmjLt;0BvoA&-uKk*m*f_T7x7>73E zJpMENtx7xl&-0rH+r2Pp!a@mJ_|1lQ z$5#Gousx z$at+MP3s9oE0xuHh%;y<`Z(|*PN9{eghz}d`0Q36@|`63ke|TUrIqmHKSmNqNbXAf z*m1`(=b0ZWzR1V;5%C2L#Pbe*@Ok(64OzK)aU%c3Z*-sFkJ{*G@#D{Q9dXQa71!aX z&_mRvuj{`w_@6ig)_L+rq${O2L+e3bg%<*^@dv5zE>$e;*-;m-@vxp{x-Jy0 zL|ogn66>+?0YBnqkKzbA3OsBHt~^f?XY89RT?7gP0-~pyoq?Vz3SS?QzCBs5yttyy zA`?EWR>JnalMLnh8@=_pU1gbeyIRB62Nx&%&2oLTOnO|e`FChOZJ4b;zqt-&$q0wy z!r>>S%Mq79<1~w0X8c2r%qSVTB-oZ+v*4Yw%_rrUP9q)~-eEvy@awEJ;nK84`QnoY z_{~O@fonDYLCr`1NK1dH^;{_2DBrgTpAy37?P76nt|PS{&_sbBsiYg~N9*UN^*{9LHdF4hE@PhTpJ5_~2J+-Zsm)E} zA2Q79`qyRFZnVwBE6;_C-r4&1_DeqfEXUA3CwI82T+%GRx%RPo@2#(zYt}_1zW7V` zlX*3~{jFUxdi%{@r}l>>$L}}LL`k|%d*Eb9?+wxurt5gj-6Xx)(gU64$d1ZdH;?_* zWY-$mf!(ELCw4JTega2u2TvrpBB{Um*_R5_+UQT%AN-$vbgX&@O=R~8&=yge|A69X zxcohP+E!2V9~YZGBX>P9ukS2VOMbIsaD&-hE?H*Mbsn1NmDi5l)~rg7S^h)DbEBI1 z%&)!Hj=%3szp2z_%k(9G$}xK@Up=+OP@ieq=vahzoZoy`CH14jqjOAH=T46|Yn5sC zi{3d9y)CX-lU%cUaG%J8_cKk+vB_N2TOQzNvDRh1GM- z;+bh(`wu)~Zq4kxVEkdffhMZ0>qctcKIy$rdOp#0+^3H8qGy`uuf4LPh1NY&cGi+z zKg*7lvO8OL`tJ$Y1x++cxYiJk;Qpg<1lRku58%E<=ZQ=D2_$~N{vz=U`X1u7FvY^} zVGD6YzAG7`zJM!vlYIbZ=!+(|Z=F$M;$lz5?(?Q6PJP|eUh%eGdb>#$PYOK8bk=$F z(U*GEXI`p4Q?@&IC68D;erAqYku@pz(+0j^Ki@&)OEtbi{c3g5Ib|j9klrse4_Qz5 ziIN^9G?An0>T6y#>02c|k-9Fz8gIwCv>xn0wv%4$x?A(<$5%@Ka!GIo*Kpwo?rFmJ ze&NbK^w2mslPBpfIHq66ACT-nvXu1w9zN;nEIJ5?-f4OUdZ#FSebi22be9!QcPM=M z`h~W0quSyCi|jMdKD6DW6=)yY`l@P|p?#{VEW^U1KfgKN(`v_F7Cw}lZ`GDHkR8xI zv`hb1I|c2tMD5RLwMWoCKd7Cd-GTPWQd_f1Z4Ua|Xb+jZlP$%qhjycbK+B@1l-qJD-t;9I-OPfi4YZM!udi|12vs(GE zn)2={<=JZ9m@x+~_XWR>ctG^|pOhb9pdFkX=ORT4+))R_W`cmt$X(fJ~06ye@?5ZL>E|RP$Nqt4~ zb{AZc`L+JyM}p%xBby9CB^L`#(75L*!T0Z9U!E17i^Vf^- zxvA9=lf7e2YKOYn!Si*r9{K_F*tC+P_1Lu1S6a^#7Cz8QpS1Zod*aJJvt9TM6+Z3D zG(K1^??-)e4oqqPQYoK#R`@(+?N=yTiTxr!LMO3b_! t68%LR;cC-L>=$u`gjOQ1Vsszk%Z{tII)65QbkDm8{4WR`r>7E|{1=t9#R32T literal 37490 zcmeI4cXU-%7RG5I0YV8dC^BLJ2gQL!v!GAU6~WMn*v3JmD@Ywkz|cgBBTNLvhZL1E zjG%*}hylfhASeS!hhzYiB0NJ7DS}H(-UQzFdn`}LGOihh@HqK{%{h0UefR#(+2!tg zlM%%xJvy<_>7QZ#Yg3B!zR90bsZUbxe#x~vbm`l_TkYOSy%KwL>fgUpvOguZL;r4F z`t})+*gvUDqW^h+O4Mmie@b+pt|6A^{VfWO@!wgvv;VLDw*G5d6mE4>kum=6Q7KVf z6O;ROJGG|BBRBiUJ{`JTG^Mc5-@0|{=E0nw)Bm(6l<4m~w$ZzBA)Va6VFWz=4;sGU z{&N}u8y9=*UMSyL%<<$@gU&IR4A4+hW^%h*5yT!=cq#+(xk=cq%#A&`TBKc|C}(@46}-@RD( z=u_>W?j>jjwK>%ciW}SL@$wF-d?Mh951P+q5yyu8&_dA`u?bYu+8?|fc$Pk26$ zd>^g@xDM5g(IdYUeI0cOoKFNi-5iQ8jG@0aAHHt+MPBn^e-eY27e$^;^WxXMBMt$F zz&Sub?Wr5i(w;i~tx`^V%6)y53wQ6$hfb`U+&O4U%<18K_s2Aoi{s*$@ZJ&*C4Y&7 zS6X9m5*&H)M=Dym>WpGe(zXO}}DKL~{Lh#fzPoB5sB%{(aU zZHhdb=0{n_;}G~iA>hqdT@ho|6$|s2UtrZCJYF(iR&uCIWQW)3wfS>)4+|LTm`63{ zI&xSCANs-9MeC>sJ6dw!v*|qPjh^pl94mdDp9(QY{+9^knRhtXoMX=G+j;n1|7C|u zeZKS`J}Bq7l+G^{|4EgcH$;(-U&D{QGLnN_=J99rMUGAJSJrc1;tKwSziGWA4uOb3 zK<%9dXKC-aua9!?Fkt`+2Vq2G&kB#QMLg}QvWD0LyKk0yQA1T0r zi?8gcD+bkH=)+!ch*LVtR{IB>Y<$248Y%cVI~)QIf!qiXC+xh!yq%9Ii=mdpaCy(; ztFGsoxAPkJqfF*;)Q(fggR;)1mR@e?5I8Rg(B>ika!xy|ZYgKATi_?3!;e4V5BLxD zjEwW3$Yq?RyuZYH{QbQ9{D!xFcG%EHDldIlhd$`fyiL)E$7%G05B;Ioh5pE=?qHty z!@T4<;t)8E0C|YImpsP#VV?7bU)l45)1%G9d8SRuI{2U^$4zh9FnNg&ex64vSbhZm zBKco8`9J3#eQ0;V4@Lg#T92P`9lr`w-5(s#AIdxw93{^Yhd@Llp!QV#v$Ur!=C6-P zx1n}Rnprog=TozWF3k0{2xY;q@L#Y;2)lZm4%2C7y>K9WmYI35XpzPNv-{GnS2cZq zXBeNKHeHvr>ANH0G5n>(|Bd(wyP-62M*f8Y82E(JU7Fb-e7xIAHyl^|10UrlX(fbX zz)*(DP7i)?LN54l4v~MM;N$dk2t*PBl+U(5;;+n;FW?V-)sSYsP<)_lrOm_o#DPhN zW_J&oDP>a4)Tu$US8?*l=*LHt`zdHzO`Nmf`<$;9%PNlYcsW4p4vUZSd79+JYaQ_< zU2>x&hx-$!7-LW9F+I&d$paf@4>XeY+S!tS1UOgZaq<-TEL10@8O{m()D`ekPcfgS z^Ho~sbVvEeiY48#)5Mp3vq7=srL4%wGt-~)%N$z0D z8L4$Q>pp3c+fQ<6b5SpSDZRmKv-Eyac5IW~O=agMUB@rM4}V1taRv&m#3TGR|2ykA z3jzYP{irk9Pxe35wmefmQ+L|&nsc)&yWO8YeC8eB-7OY;@2!2+_n7MYfX>|x(F3Br zbe;Q$_yghkUKn|_ad>Tk|1i(%X=qQie`pg`wAw|ywuPVj&<0{2yShuygVGy2u?st~ z``_Zje~<&NjKLTAjLGNt6|^9n8K;9oK<%A+XKC-aua9ycj|`<{+@COepVI$PmDi!N zIA91HHve{om2}bEcPbb&=!d>&MfpJaM%hNWL^<`9@EWD^XO_w+$?-<1oXHd){H;{h zd?S1cYS;uH$}#X^-$@VXPrAZ>@MHcGE4`6cj|m6l?bdkSrrpmiY?$gpKKFOjA>a^* z90c}Qd56cg{hFuzCl2t~OMGFTvY2>7xlViNCBX6gfh?mk#CtNAChl5ADny65A)P5JZ^G+UU{SJu3mG4rgQuY<2P)}HZ4_u)JW_Y zTkrLtS-7|FeI38ZHv6RSh?A42)Vd~UHoh~e+~xs6vq9|_zB;I&_D`(Nc{8i6#A9sX z&A)sN)?4kMAFrSF&s!eeW!5|v(`oF)J%-2ND$+;xd(ffMXS?*FJw;oI*J{zy2OMaZ z@p#QV_+eLD;gKNPMs$zx*k#c*8iOzOA^Ya&xkUi~As%u!{Lc0x+Bx<-ahZ6GAK`cQ zJm<&uuN9s8Z^+CJnzcHgdv!jfmq#0jdX{tiiPmwSke{DvD{)RazqFNRrM7tOi(bpZ ze`__AR|1CDAuA?b|HA&$W)%Gd2#i)o3g2m!0SX zK72h$`;E8)E_QtnE>Qdl|AMm5jyeP)0RgqA>Yk-NbuoN>L>TvF{{6+3zh7vIJbqc? zmp#TC!d{xnj=PkOwdG#5Ya1P1G&pG1>7!-B<2vy-6F=$UcF7M83dt{MM+Y2kk)Fkc zL#D#(km?h>y3LlGVuMahi)yymY*Dzr@vXN@)*GK?It^_%B%x16_-({K;ZnRoY|)8* zf`;&Zr}*26ANfy6eo@J{aVsi)&KGXzV}wsJ>GNjc-t$(*e&pjV%~s*HOljTOmf0?rIM!y&5b>;e0P2MY?B_7^x}`5 z&cxNCZEkO$-YsZ;@T$BYHRiC1m7Kv^r#}_WiI$!HWEYP=;7vRPM{p(|5=V(&(1LuHB1t#)iTztC+|?t zG*`XTRQ&hLzlQ5v?3=Q~SMS?Jrt7fXPtNTz-Bi`N`FcRL8SN`CH(PXluG&P?U-|it zx^Y=%+5YtBhSv?4sBUY=+&?d9E^EAP>f)+d=3u$&C)XboFb!+|9OEAyG#PQJ9~>W% zWrFv#eWGr|40BTLo#QQc_AXpI+pHe&XtTukGEDVRgFDo!cF;Vq_44&6bH3uMawPWT z$b{o&{FZWY3o2%tk6v2Wp;zA%=Jt$x-x_l~XbwqES*;r(z7vv*oF!VveJV(9L&=#f zyMF2#vg0P{JyUj8kX;94$7N& z!*AJd_LpzvGxPqtcUw2skMV(1}1b5nhweQ?9t=PClzA~NWPEDHp znyndO<|vp=qOrVZ3m^Pgz` zQ^|=_J11JSh2(xFzLlaiRL?e-oHn8#NKQ+wt0}&+lDA6hn~ITC2|ru%%Y`W;*g2}f{WD|{amuIxic&EKkd;wscW2b>2e{uIUjiJpry>119! z2&lbN=Pd1=i{a~|W^$wEmOI@gH|6aLVR@s<;?PHi0mH|?l--1v^=fadP&r&#_2R-+xEW1v21P6xNK9yD0Y*aZ#8)uQqpJ!Aay-+{Y_6ew*si$)1 zVYPjZs;pU~GKaQL6YUFKdeUD@>DSuY$8+t-uj+o|zHGm5?etV1&-uB=$H{jHT!;t|PnuZq zq=m-h4Qz|oJbuPJWgX=@Z*Q5WykkDz$~&}`YISd!dgI~@^O@pdImO*IinDS4)+3Kp z4;a4I?5a3gMe!lYYAb!M_RluGj-WkMN9`cqj+WPJ4B6uyQ9J0!C9yv&d_3DcrFgzn z@%(zl@d_&!Je>0KnLpi6`@fp>`Cf7Y`YNfJ+DiNqP53RX+Dd$_(OCEp-v-LA-w6+J zDXTGgg2yBB?azA6LRn3mWZxV;w+Kv7UGcK~CR_V4R`XLePrU*EY|XQ8d?nPJd@D*@ z>7d#f_{CVAmml;Pw_Vq(G>RFR`Q&kPXnX5Nj?N000i}0Mm@znDZk-fAzfR(DQ@X~N z<4Qe{VJ0+an_aO^z+By{h&j;hxVc8>nYPj~wS!Ko{nJFRvns1S)J|anvS*GV-D&2ES!xQGvSF2MKEa_P95*vyGKJ3*!Y83bt)r!KeyMNnkx6Y{E*>x` z!sl7-Bl|^LiMU6bka$NMiN|C7o45$&@t64ASm%Mduc^k+yEHzxr@7!spC`MiOYs-G zUS0R{R=;1U-2s_^#=6~eg+hyhph0|bO&6Ir4x0xr*($U2}x@+R=k@9|O4U&73XgrxI_hpw)=UH5d=y|+*I zb5x)E}3 zzk%^_3EkuUFZz=TpX}*RD%!6{sLPA~=!kLt76rQb|LX7HuN__B-X@Xb{4s@-3ipUl z9MJQ`n8?R(_K$xibh=P0Zv;I#0CpseD-)(&J_sfOm-ZAampx%xSw@-4kwi$@k$p zfb*a(%sg_6%-2JYfeXoiql<&gg<~+g_V9UYU%0i0_mes}Ziqb_-TJHD1CIfZf%CwC z`}il5GdQ%cA1 zBs$#s=T)}O@Hdq97Whtl;g|JiGKKk+pS*hCJ-<8#a>9UH9y#ZBUCaTAyEc-qY|P1I zy*?fTe-Hy_Z(rM(gP({efD=B zNWLkd^$M~tIL~^=4(+8i4tpFg(0-#Gd+ZnN*oa;Ff%pxQ4>jHckAXbFfaIOVXOVZj zuaC~W!-N%(9)t}YJNY$SJf^ayoWeBeptITu(}Z1);f?-;br7B}cBD-7l^Juz;6@95 z_{;jiU&=-|en^|(#m1QM3Y*arFJO(Z(1Jo_ZK~__DV_h>(&SM{>o%z`H zJNpVGZZPkhqORwG$3UKDK=M?hv&d8V`RgOnaOmySHnV>66Hm_`zL2ksoZl422y@x= zzi5{bX{Fc6Ful!e5D$d)60fEgjJzXYzOCG&=G`CdFu#kJ5%?mUahw`*xz{iI=onD^5*F= zkS7?ROt$j^aUiGdcKJ$th#Qo!@Aa6VnOZv8Oq&)o-ziS+A3J1J+24aEdeWQ)`?J1U zj8Yt>+-BT9^~3%o*=enDTxYxN7Ly(P;=Jt-eu3n}b&5wE=ZTx_IdRSmux`jx)CtrN ztefoGsi%6Y-uOoAi*=e(zE7cq+jp3|%e>RLP~x)i{uuqtlsd(W2h3)j=drNd^oX}d zm;c$+tI^)sk{C4YwGR?z)tK*}v)8Qa-1zFK>OnJI`(#9`UBi;TNHeiTZymhxyP!cP z>Y#B`)USf<4wW6o)z@{_%I*N!A@AeI9__!!<>xxh8$Y(m?>pt^7M+Jr_-3ENJ80Ka z@D85y4fY0HPzLP#Ox(5GiNmzp`;Iz%cjhAxe)8;lzV^`z_B(6;;fv8a{6Xt(yD&}2 zdcq$2RHQi6`Q=r&vj50K71h2-vM4eRb{JQ|l54TUKIFA6@+w~eQfITD(lu}V#IGmh zCw_meevE@R8`)Q2Xx`fl^*FqH`|E;k6Yt^{1Cng6zO2B$}%2fkbT;# ztl2G|u+Q~#3v&NJ<^?~b9p+D3L)M{8BkfgEdgFYo($Nv|K-$@(wD-i;4zSQ7;wuwZtc5V@)YGZWixRKc?emFxD`_CgW=a2$V#PFUo&g?vY?rv$6d;9 z;(okjD%y$bd^N#iF0vBWV?LD8JWh~5d20*0_^+T@^j-h^yR6GJ zop%YS*f_@0}seb)N#ma>{Dc=Xz|fT2tU+kHo~io+&2%; zF9Yl^kgIVFX4mevzwCHz?=Rvt_y6;*aU0Vzv-XMB=U&yXy1qj`qMl_Pf2MKlW1g!R z$2vtG!p`jE=r`B(T^jy>t%lb7sUMN9l-(4~2YaPn3cbc3pt`$Af%VUiymXC^`7F~o zH(7~&ZOcl`2bm0h*oVXw8{xx7>V5JbJaa!i^ccts3`m}8a29zgKYV?J`s_=4?d28K z7MjQ*RpNK{o@hwd-{`GR3#Dbs?FzMdht@fQOSan z`UTA zxvb932`>)|m?D+^)uTqO_wm?AS^A;UUTgKQr2g24r!BJ2d{|den$Lye$Gfh_KyGHh zEzZ+^tM+w-_7U-i{Y75*NwQDqM=IM4PH`?$^3L|ZwO#z_XBmd_d1$+9OC-(=nj4;|@&1PLnPy#7!ppyQ zJr!3Y+uqhGrDxFWbE^6~|2#w=0ceCuK%MNmuBR|S&-aPhK zkY8)$2YwfopZG<-vk@L`gjXAL^IFPFbK^eU=Cc0ndVqalw1zCVK6S$D4zOvU%X*4)s>eo3QVdiW|~Da)_3mP z|ERewwd4G8M}h{KsH(<|P`~}M`+)2$(KxPCO?I&}Rr1$9`O!%8o*_T0$gf}I$4dE~ zEqgSrp!zE$VBkJJ+&>wxP-3W+%t6e z+Ho^7%!;&$nV;4QgzNcsYG11M6{=S&O3o=Ryi0b!P(QGSt`jXg;1bz!H0}!Zs~~%; zWT%zJMOp1F)&8;U;0M@JcJb?8^`{iNibr@~uj2>BEBB$B+Toe|Yone2 zIOcf*EFyd7V@-NC^A7`(cj})--pLPNA2rc70=e-mw0(K2LixE#WpPMF4j5z~%5IH! zkbNlY%d1>Q_9?Hlj0%r^uAf`*S(Rh2iYL;|w<^nO$q!^7%B6p+oI>_ltnz1+$|Gc- zA63p!?jZZ5sjOM0G6(y2>Avie%)xwa*8IRfiodzK?u4D4epov**++Qpq6KGo zEis?Unva{T^rhxw%S!w>;Z7_5<5xNHL0&CwKR0V=z?7fd za$f!TBPQyaofAs@In^XJ?vPore!yJOH`4sv^N0yQ9=pg&*XnhaEh~Mg*Ii%fHJI3T zkb{;w$<2y=kR?}z{+x-+M}IzWO@?v)a_N1yG(Bq4cdbcoS3Ny^zMAGkJ%AlsR&q2S zTUPo?^LfhR2U%&!9lxYccr{>lil4#ar){yihe~AqsBiYc$!%XL5-`t;pJ%N5uNl2dhGo$VWaZnDyyl5cET3HgcrGj_`L zn^zs)=tEW-_U@wl#|%7Tswp4Q52RnB`dz7h%m-O%sU<5h4_j9HNd1UA@Pn-Mfcz$3 z!_NfqL%)vd_lEV{0nhv>kotmp0vQRJk^9BI0+E$?PRJ?dz8&DrC@(N@LY_*f_g_JL B*&hG^ literal 37490 zcmeI4d304p62{pRAP@o^Tu?^^Hyk)@D#+8mvP=jd;64s2yMiXbKmvwEk;M@vqT)jq zm1P(~M@10>?h7cw074Wb85RW@$PnI$c)$kcZ zup5fuUId)D@2Eq-A>a`By$E;`07V|x(32f!UL=Y zL!cNCPr}he(%JcgB-ugMy6yi8D2E5mWLzUm+ z;8oK#I0=qO{G}DGLiUEz-a_AzFZ>FA*f-u8ice`>cW0MFpeP7L$|EN|u7`P0;;v1} zS2it*vW~|g@CPB_El^z%ZPgXj8}yq+&Dn74gYz##WNU2`8f z%!3d8;Onh<_*D-}4tzF!S^6OtiXAr9d`BDtzXySEeqnz}SlI8_$BeTtZND45bYImK zAN|`0#owyR-=O$Uy7v34qR6-BA+NgRaG!Df+4e)^LD3t^bFl;ap_=E2L!eY3p!QDl zGaMBp*HDLw_hZP-hgCN}eMQfujh}rlIa7 z58xM!v)}MTdt7jOw0YRiv`Lv4a-(C%r&T7-?9?*#9UuI>kJPgK2>FZUAARNj_!aut z?=8&3pI*^?{EYkfRha7eMfEHAB+n6tK&eJR?WtzRX-}QWUmsD{huSTfX6=**9-T9C zQK7FzC<{)&Ut+fqcJ)6RrZdet;XoRwK6`iha%}@<``LXiYX9z*Fh0M0c}2>XKjwz7 z;jb$G@5N7g3Z;Q#@=FR};1f!BnI=Q{c-K{HIkD1vKFUwha|p+Pp$wIs9{k{hT<~EZ zBEO{Igi_s^q~1paWZ$@ zgQIKw5;QkXp0{vs!B>m1ildanLo{!n_$Z%eNKU-w5l^xtx1!|me4AoV=rujlK(X7V zrD?65Ek#FwePicw#zS?|vHZ@sT~|;~F`g-Zsiu9px#s;9Qf}B{;$vTHULkdv*`a+r z?y;#2s#FP>EZt}S*PK!2)iE`9nZ_61QleT7%lIy6kkeK3@XM-_J6v+cXkH6FCsT3ZPsH8@$#_??+_E zM%mq7c7Cb*;D#SVkz-TpHT;z8qOYvuJRK0AEkT{hda@rvZOg*`MZ9L+*)QAjyZ!b3 z$KUkb)N$cnZ_SIo`&HiuwC^^H?h;MZeV!lU4;Z@*L>uQ#wP)}j#yOsbKA`pwZKB$0 z*C3Z;TllpW;jj7`$F6?TbGP)yPVB-??Eatl@EhdV6nx2tT$9hi7kWDE8K;9oK<%BT z$7%1lua63Ek8Gu7Ue^h;w<-PaS9u)@0|7(Wu=z>f(85nTJw_<6plB8P2s_izj-c&8 zxkNejweT9N@@J09D9Q20s+`FYAN)6}todHy_jC=LtRED7Sa;F{`jf7(AN&|U%Q_7q ztsW2#$lI>#m%rTp)S{N@KIHR!M;!tVfzpG(w-NHrPK{IklaF}qwet?Auf!Y5b=pJE zDLyk#2zt=+dd)8P-^*r=ToyEbE1vUuJ6P=|=IyfduVA&6XoC?q?J^pD!l$o9pf+j=l7rr{E zrS?yZ_IZNUR^m0b?E2pxgC(gQ^z)T-{&Vd;+sx|wqkB%6{H@_NxUTe({T_6r^w}hR zXiw2r;#e(8`hWxNGG4D42S4nh9<%AU!eg67SL+&lp{$#u#fv~BKRVfQ;xh3VKf>?q zar}_`tox@u2WRBu2hAGo&+oK9*uS)asAt*7A8Q`Z3CRx>IyPqi(pH+2-tnc+6PJen z)@mrv1q{a_%O!WC^g)m4XF|vLgVZLfRCeu?57k-aTQ9y9n#b`BcA_80xwMsb%1-pb zzxaBP_8axQO~J#a;A+$2-BnJyL*S1^K<%j}$7xTU311%(u6;TGdUpAFi%hu(&rROi ze^R*ZlBu%eFDh5GI9|QWTkl^yJZRSHqh-S5TJa}{pYHbck{=uvlK;mp@`;}Y4%bT0 zO2Q#W;dMm)C%h&J%S?q~hh|hvSYkFPTwneA8&#T4$}>Gjb{o-kKz2B7><})Mn#Yu% zJRoQY?`_22P5j7zSn|tDzKvUX>2qRBiiZPzjPR)-eO@b@^z!nU4}6@`d?mb=DvjMR zbWo4@ck;}+4O*l;J0f5zo$WtAdh}Y~FQMZV!!N;oApUsqUm$+ubDBaqpQY!Z5Bs)w ztanl!0wo&(;u~?kjQsax$7!EahV#1I`l)G&2j*@ubuPYp%%h1}rjp`V?ewOm;m$0x z@}jG<>o(Y8hAGZPslBuFpIw%G_)(q-wV?whwd)mCQ)dRv)el_w*1DScCM!DS*`IqJ zi>u|kTz_{~-=O))tNU)n@%v1S}n@FYDsU} zR~N{x<+6kL6D2zb$u8m~xDpSqQ``i9J0Fr4p(T5_l^pF^C-Odd%^oNJ+T*NO2h}_6 zwa&N8zaG-Q*fDLhujvnqO`lQQADQ2Ornx}-=9?k)X5D?xGP6PV=c`RL^Mzk;XcCua zmhQ}Y>Y*k9Q?c)w@wdMmH0QS6IDN?ld1iNwE2lOa7BDRv{2J{a7c|*%>F*sJoo9l# z-28Bpmf7a8+B*k3ZAmKIDBr9aa$iF7+u5f6*x}t9*4u6F`0D(1hYP;qtD75hcud!W zX3~ZlaSLnbn-88_+dXmMA#;89t*?(i7&Ln%Csy-Di|>%+BIg~=<2kh?x25FFm0iE| z3E6Rt^qws{Yss!%vSXF(&XJwbw}$NMBAiDF*Sf+H+`;QBYMLy~aM(tZ3Fe}&?g zQ2dtlhnAJxlC#rKji%sAorZsbKjYvEKln3F`>#=(%`+-YT% zJ$j2C5WT~SyPv4eY;D!0LdWYkeB_$(&*hnAIg|1~taq#p)JfwXYaBUoYUf0Wwv*h? z#J61ZV%4)9B&Un$d(x+)=3OejSjk(d`R&EebCApUhmwOGe7%$)x!A??sjtyS}P!6^`J(M%Q-=SJt72##?EeI0?1an{|icPt@_yQ?qN`Q%eg1YVR~YPJ8D}`1?)ZR8G;xS*-Hsag|5UGz+zT0xD;ks@%CpZJ+%rYgVhwq3zR7 z>%!>;>92+KgR+jrO<&Y=CvLra$C{bxz9?Pyke>0Wp>`f&{!^uq>PlauC_oO|4>~Bxd{;wx}_DW7b zUnM1|t;8?UgimSJR^n@owu-C7x1q9&d zL5J1;X{Y0?bJQN{rnVAaVFXnt;$PRP9kj_yzpnDzdEw)-)$fn_bjYea^T6Nh+;VNp zLuSv{tJ1rkzdL*%eRys4BxknTN|n_fyGr`JsJ2oawUs!o8)xAYrMA+_wmPiGQs4aCDP5kg6fkMR=Si*OT&*K*CE^}!aq1P?NW327-^4{|qm9BEIgGUu8`lu{~$tG9P*G){71#-I-c& zd@(#!593Gdmxm;;tN7T*lfAlo8f?z+(N^NCA=*D3)V^siN`1}OPWaWhsej(E;y{Lv zcF?Fd7T-2@$U)Np|J89T{9RR#j}%{h@u3gxu4PtRiR0Wh(kEVWIKJh#2*Bqf)$zB> zZt^wwJS~5s-3E>eFR2mRC(w<6b^ndQp BXwd)w diff --git a/tests/expected_merged_bounds.pickle b/tests/expected_merged_bounds.pickle index fdcecae0e3c0e97fccc1d201254aecc89801a2f4..65e5fc0fca01b3b0c82c261b6b6a5a7703a0dc44 100644 GIT binary patch literal 37490 zcmeI5d301o7RDiLVP6$w89{VJ6l8NjTkDA+Vj`R9xFNd)B*26ONB~6^N0^9;jR8au z9bwcFK@7O-0R)7BfXEVMkwtA|BI^VdlAe(Me(gC0eF#T^*m?YsdtPnt)>rr5y7krT z!`PD3`cEyE|8I=9QF26VGjH;RgA-zhCN}NSXGmPXrm+cw;s^AKi|du>O|H-*u3w)a zgNMb(CG?5+KIcs?mEY2vTxM|JK$GXZQN<>DZ!O;2`wwqdZ=;#Xa7R9XsCP&-<7EC#qPyxA&wrBdY~;3dg_gN7jztIbe#p zLy|e(mih7vE9)*cyQe+!#N5$~O@!(pDpeQN$TiJcRwyxbu-|-NHTA>nu{ow?KV_Zh}-uXWia+|VaVda<783cvY!@8q4`zgssu)l*W>-KA%T#DINmlD=y) zr{sRp&}XVYAF!*78egCA!YH4KP8z*#;;bCgK+nfcIHBKf?o>H;Ra5n=HBN0kZ_)Ve z^xkjv94URk2YIgWao^)2@MlDzgZw61k~y$X*yTJ^Py4Ushd0YF%RE@4U>4{#Mypx^&|D)-B2ypRG&A%^X|{vOwC4Z z4lhaaoBU{$7mxNkbv`IH@u115S>b<+N9LMD#q$dk&l@R@S6R8}q2w2T+YiyFmh{Qe zIJ;`~ZE*VsJB;(}^g%~n*CkUro;{AXTOg~tVwl_c3; zk^3xM!I%A=b#vv}MS#5G75qzbs`8CFdS*XOR{3m|vu+)gPesX($PY%z-m|7}_q6z7 zNzgCGs{DxTV4Lbow}~8{`FOVReckcl!*hLRc*WgM&q?x`YSY>;XcnJsA}`-Lx%}TU zOmeHPxiy;k%*BHu%#Z!DO+DFjv2I--Yu-A;9K3Sw>o+{K+Z>AC5nH@zuDN32+vQqJ z$uYIYC-rF3Fw1n?Qg1`v-fq5lXMB|}hxyHy$`ggQb8|+wg%h*=W|@ z72`e|wl>E&f4}gq8(JSVS-aMzcCMEdtgkD5$Ojl_`kR$^51OCnw|&=5yyp$Saiq_z zSC8M)ylRfA-FEAYrGLpWUrL|HEqwfU1?Kf}+n=%~zvMGJh0h4#)2VEe!{y8Q%-loM zI=xuhXPyy0Pg?7>RXC>$ck&h1i}+}?j>IwIC6aX~k0DNz_uVP`AP&>+>{@f-8~t_R zKwer(&!O^FyBwTZho-k~pILJ95)bjWlh$d2#^s!JSnHdfGHOrekmVlM^??RQV(&}z znWUsZKKzvI#4Guz#_OndXE!jlJ&x`kfk*)e$`RW z_?PWh?NqP%yTV1_4@JP6$3W!~Va>fUC1c8*ptEAiLS|ieg3PKqsm;4p|IjxUs=>aJ zJLA5jsa2H@f30*G8Y(xCr<>|Z+o7SVsr?4Ehlbi3V*Ek{LY^|D4WtKbPd`fk(MnsO zi8d;oxk~8{G?eD$q$%x!hT5d|SF1fV6yuMv#t)^>`ND0Ra3H*sHosLo_KlSlKK5j& z9vbR1ZATrnO+iDIuH~&4IX2A$PS8;At9`WELqjn>IFZ(&4|WR;b-wU%-|He!m=Pdd zBaI`zQ4UR1Oz9T!2Y;vi528~BiaufgXd`;1yy%oj=n`E^a7_W7V$(zq>$=0HiPq~n zMA$l9gFq8i6y204dWm-FqMxo4-Be81Cjr{WXEtlxP&82$>21?QYjxeiyvt~w%!_j3 z6F7o9cp||S$^KiIZ7G~|BR^q%uz%Lk(Q$_Tj^mG=C)hO6eexrlCh9NRqPFa2m#(>> ziDrmSd0F(y3emB@h+e5II%U7EubPPt{Z`jpp=hEtqH%4Si0ib+M6)i@H5&9zrs$?R zqL;>~-2u^0_lR!#Sp9fUU0u^Lj!hGB9cR--KWQGTH190U6PhSWbljuLmud<}=(!(6 z1A*(kS_g0sMH7)%UXcJ-LbcT+x=pSgJa-xf%XP}9Y_4OWT zALJ_Ky^i#Omg4%BYuFBYf27_+c_^9)y=|I^>t6ElRhkd<4>*Hsgm8qOOB25L3Rl*l zx5~kpI7xoNGx<8}j%59j&_rj$CfzZM4g#WgTAYI3$^TXJ5K{%}M z(dx3udU?0l^}9G*68+{I;Wetmu2IRKXPTLnQq7DResf$j*@P!k8kftvz1pXCyCu6Q zje#X+{Bt$$g1(3CnDwL|_+TgKzwp7MX(_$^CP(wHAiG$i@|t?4K32!bY~fp8^=w0L zY})he;?}92jjHGUu5=M_5eOdypp`1i|M7R?rd`fgJ$Ct9iZhEOiBoaI6ZXFn<2Tb4 zpNU(1btEL6$d|H8M6~r8Xr)dazqj;r>g zsU$n^qWmynPW^@6d3%Fj_Z{1GMZ%_^4jHa_V$0k(V&gu)fmZrxU}UdJQ}-HZ*}Bq4 zxH!l->64~ww4oY@enEeeKKR*s=}jI7KKSVt)k7;m6HNJ> zNc;%NOYtAu?|7!3_M!ZXc#Iw4Ux8m`{C3{K4xp7b%RaFiy~n|C$Y-(RPeuRW$DFI^ zhn;f0&Ny>Zqh4D-Xj$+YHK3=@kbgrfrAr^it0?+tsNzKF;%U!3df{3R`Ye||Z)^Xw zap0@n-MJ6k`_WVHc)pX}qmPXb_~AEv9ZNs(uqC*1{>0B%H&?m{gaZPir&^wZo+=7o zAB9>%!ZKlzFrfC%FQOL+m#1kT>NSPF99oUAN*E-aDX(+~8tPZ!a$IQ_G!*F|X%fd3 z!n-|wp$A3^LKl%9gO77l#nzKcf9PR;*azvNx;xNNv=8E{W1Brb+tQx}A6HKofp9_q zIt72nAMsQ63%lHIZ_`A?2b(5>embCQk0>2;p@|kKO{U*|wPXCLqK~50kN0fXH30FL zc@eLe4|+2%B=cLZ^qaUBPHXMXrQir)SJ(^o!M=~(U?-u<$)})6psBFabkQZZi$20X z#_5_X6iw7r*JU$Sx$j89j_EYdDpJE-r-KB{NZqtRSCF{-l5_gH$*b(Ox+T(ZF0e%O4 z(j{uq0cYL6Ju#wVphcjAvlQR4C&t51g^@#^Z_`BVL;SdzYglL^#-Shd4&y>M;fGvz zES6t}qKVL(c`~m@G*9LYeFuK{6S#sGxYI76_cjBr;EtaaCfjj!DY^)V-f49TddL0x z=w!1XY$9nvxYDyzR>yPdi;8?GKYaQOGzG_o9x5kX6ZUxqH`)_kk>P?L>6xeU_=`uh zT^<(|Ok6L}k?$=xSRjBtcKKgJx#*e8-e`KWcSjHv(fDic!_28lLT;U>+j{tT- z-iuv@D#xy9&%Plzow#KO%dh_QUmon4{S3RL9kSfy+0m^N-}hke_&etm=qavS8Ml}G z%GQT^#z#_)K6d_&Um=Mb=*@e;#g^3TeXej32x|mHPhE2gda5XYeRP^N52p8om(#Ql z*J_7wY;PZ=D}-qQ5UPEV!6&dSEed|L3k_97@d?*;cjx6IP!I%y`-mMsiJJv2*UKCv z{cXv3wk(LUZW|YYbA*5ceZsL#?+yNl2K%L)V<-6-$3OBDj(OxI^n-q)UEnK&y#ESB zUr{p-{b z!ETZG1@+jwCK22R@l(cwj=?V|XFU9oacs%^=#OOHNaDut0oCp6B2d^65WRElDd-*d z>!UMXEJr<%XFNcdbYT9-Vx&o$ALU3D7a764kxbkPffMzvbP;e7a1r>u2sj!5$#^^? z&$OI=p(NwE(nX*!BH&0VBs%hpENFSCyeC*5V!sdO0n{V)V)Q7eL|<3B2%Jv@9K9Tg zE<7VoS01dl{329&(4W-62_+fNmZ93K+!ZbYE&^u(0nt;}or0cnzdj0yx1$2ef&PBYv7irdJV>Tf{J}l->ej9MzW45Z z_kFKABa2Rcc2c1W|AxEoPASr-mOG_l-=seMlWTVD+ApDd%|1!J6MJ?}Na&pGPKoH4 z(7kKFz5@~ylDa0kUvsB~UeI)>g!S#_Q+drDS7@xeN#QQ;zq{MI?~W_ns!frx?jE5j zq1_Ub`*%P8Op#|EaF2V*cRM7du*==Lb?X-1ynh$|iz}4q?lP|N;Br2ng7N>@!M7yt z?K#1;kfe_fXZ`by)qh-U4ovR#;_P9IO;O$75}x$N&z&<(k9QZ$-9Ox8nwR;YUeV;0 zCZx?n9WuImjh?qQ_gE0!V9g1$>W=#|qpI&Ql~U`N8izCd`jLaqOb?A+>bGCS_I~`h z(K2QD{uk!;oMC3jKm0s2pvue+l~$U2)DQ1KUw=@N8`2Cz+=Gu=Yda z`(~OYZ*JM(GHm16H*?ILx@I25sh#UT>i*ZVbDOSr$UgXd zXPdXJ{d=U0jg}Q&5JKIIO zv9-tOnsuj;?AGYfGPVBFOtVJ$u&nao2IbjEck59nt9s1$dS5r?(J1AIq^c9G+OaE5 z*ng_c9JbPHzWYa;rC%=3F+0mWS2U^79y28Vo|4Hkyrx#K+dtl1F4vTs+tc;&&*`Sz z-HlHz8FF!5J84CI8anobIeAmWPm2fVnirMNmn)y&t2|zQ^`fU!-n#f)XZce_{v4Mb z&&}QLX#VLQ6Z`7a-lym6FmpQd)|4fN-zopEca`m2TJ1@SXN z^Wx~#-L5)6Eb+(1IJN(z`LIj(qZ)?~%z7c)9NX3UnNzbpW*ZX4h=Ue$?0+iIBk*ut^dUbD+e zZCL!19DhBz;qy_e2CUC9z5Z75(Fdc?m}B3rPmQm9)PEm;Zk0cl?96(<%)Z|9yry{f zjbq!*_nQ0U&pU7b(x`T1j;Yvm=d`7_=9umBXN(P>u*ZD;TG#Z)qvL0JOqTF@Uiicp zuW_nm36Gg~a&nuuig`?m@Oee!$arN4XYyXU@MpXlYWy15;+uSl4AuCPpPTBMy05ve zkq_(o`p)M3XMMi3|CW|Q61i;iF7qP}daaoQ|B^6PPpYsQqJ z`?LD3bg}L|cE{;HPbPaz@(^F%C61UcBV;dLdd%ZVR@AZTyVG1>soq`PIUy}8*K5kC zer_?af2WvbIi|TJ_4O+05m#fT{CWG@(`l~es;7s4xa9HC1F}tZ;`jWfH!m*HkC)ys z=~a~;{fS!RIx_(@N=q;{w}wxr86_l@TDm zh%?5OxI$8IP~Q+|NZN@j;x7M1fa6gxAmGkppz?^Y+U_1jdxV{*vw9SV%)0+PnH4#% z@yC(D(9z5H0os7@rW`7xayVV(^1~{FWW}nYvVD`vi}KRnB>h1BqoKB`{AnV6Xe;bLFZ<9?cHE$$ zZWKcJ!*-JG%m;=0%wHYxh3klLf+kv_<%s198midM?#kgKH@VK*mP3}W7%gio&)7dI z``aa1e(;BRd!z7i-s>QcpAjIxk>@GH*e2QAi7(1I^uHCI(p~fk&o>smQc-kDG3Ch$ zI+j?cV~Rn_bI?SGMHBVbaYwRfqWYqVGISgwY%PvKeiFSCA-V~A3B4V%TTgaEb$sHZ zeLQB9?4FPv`qfPRC@;UCQ9o~0zgDXs;1i~P_E*0skHMRK2#(;)dPp86eX{a*cbZGTowiTVM2AFAbQ5iHtLDvjI_BE0`}0Jn%n*ID zLUe4X=#^VUrySMsRkY~POdWGU6P*)Hlr5TQy=dI3qKTdmO%$c$v(Li*1y6f?HL}GQ zbtd@{+@XJD14t3u8cz`wL`CxSCRI4zr-O}k1TuEYr3giXB%MG zC)`Cz4sJf<);^Ve^j|6CPMQ(R`p@>%-a{b=N|BXO;>o^JxN=C zJn8@Hp-)(YL1_oQ=4a7L8=edAJZ{nf(@XVal>AZuE#xryvrGOY$xatNS5kh2NrJ;} z`OWbp?chhh+NvIlmxN~GYd+`)?JLUM3$C zN5q}Io%v&rFTMuit2@4akzdw_iW}-#=J8jeZFr9_&OQEkVV*L7Ka!oV?NM*zGB*`ctI~*G*y4CJH-lbdSyh#^)8O5*UF#8vSY`gtn$jf-1dL@{N;~a zhlMNtfDd2y`rpU% zjPxyJXD$C%&=2&QeBVyLdMaI$23ZzLtK8YI?FYD=-TLdY3w_F@@1;+Ar@X{|z6VAQ zf)Cp&wvAT9Qqkj!edc1hVteJo(PMIKUmu|lPS^z><{|d~U+{7KbPxz81nlKH@kqQ9 zr>qy~Qqs%W2s%9AH`?9ocw+%n3$l-tltC!|NYKV5cW^c?vl zLv}-Dhxe1GxTZgp(Ma+%Wim3D#@gx2)g!>XVjU+gSkHWQQo3QDpifu@a9GA7T(L@}l%}|_-)G-?L4tC=89C2Avc87>=8l~r; zXVY~|*WX4HQ7`RP{DIdt`Tc_W0lmNNzIm{XYGRG?6f}prxRR zh#%S)X}v>2yFxGVBWUb#jDkMx9JA0)zkGkRB;D|}KK-O$^pk!=-w{990avcU7yDeZ zo`X9wKg=1&2L}PsJM}I@?>Ju{T{=BjUXf^!Zgg!G7Yoh(f+Bs}0gs`qV|xSpPda8> zLiz?b^hv8MYr%paxVXxUzHLz5#V+~_4lIMYqJO~2jt{sXg9RU_9}WTz0+%8{p0KYg zwA%8b^_O|-+a7xi z+hy9(XFaE#cp@H%ALiIQuV|I5XEp09Y} zIrNYv#?NROoxIGi&wd122|eTz+4Ipp9>cuHA80G|k@&+`h{PH9i7UU6e4R)ga(sN9_5Nym__ldYOgfdAr1(U0#i|j*f%C?}UK0K-XZ1?Q#9D6aZ+5XXP-|nd&7PJoQuJ5azi(ge>hk6D*{KF2k5%pjveU;ArK$2}W zenR&`_s}1t>^Qmc_Ml_1 zOFQw5zu2)$;+5xlFa02{h&Mg&2nT_nKtS|P{malh&euoRJi?ZqKwh(fAbEiPBMVU` zseiO1)m&f%{YKJpD*#U1ccg=WgMfp;Z$-e80VMXgMqX<>&jpg$bEJbnenh~MQ%HQ| z8hN$tf$AQAdw})cUk7j>sTbqN)k^$zq=UeXM8ML^k@&(j@=ERg`?k0U)b5WbDOiCd z_Utl{zSfOuQFgAAxxKyTP_v0oQZST5-{xF9c4a2gy0 z90Yn_QaksNJa@H{=N#!E;2@BPz Date: Tue, 1 Sep 2026 16:31:45 -0400 Subject: [PATCH 14/14] fix: gate Davidson convergence on the operator scale --- src/lib.rs | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 5c66cb2..36c8b62 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -17,9 +17,12 @@ //! sparse-times-dense multiplication for complex scalars. //! //! The stopping rule mirrors `pyscf.lib.davidson1`. A cycle converges when the Ritz value has settled -//! (`|Δθ| < tol`) and the residual `A·x - θ·x` of the current Ritz pair `(θ, x)` has norm below -//! `sqrt(tol)`. If the correction vanishes after orthogonalization the subspace is exhausted, and -//! that same residual test then decides convergence, as in `davidson1` (`conv = dx_norm < toloose`). +//! (`|Δθ| < tol`) and the residual `A·x - θ·x` of the current Ritz pair `(θ, x)` has norm below the +//! gate `max(tol, 64·ε·max(‖A‖, 1))`. The `ε‖A‖` term keeps the gate reachable at a tiny `tol`, since +//! a converged eigenvector still leaves a residual on that order, while the `tol` term keeps accuracy +//! tracking the request when `tol` dominates. `‖A‖` is estimated by the maximum absolute row sum. If +//! the correction vanishes after orthogonalization the subspace is exhausted, and that same residual +//! test then decides convergence, as in `davidson1` (`conv = dx_norm < toloose`). //! //! The Jacobi preconditioner divides each correction entry by the shift `diag[i] - θ`. A shift whose //! magnitude falls below `floor = 1e-12 * max(diag_scale, 1)`, where `diag_scale` is the largest @@ -55,6 +58,18 @@ impl CsrOp { } y } + + /// Maximum absolute row sum, the induced infinity-norm, which for a Hermitian operator upper + /// bounds the spectral radius. Used as a cheap operator-scale estimate. + fn norm_bound(&self) -> f64 { + (0..self.dim) + .map(|row| { + (self.indptr[row] as usize..self.indptr[row + 1] as usize) + .map(|k| self.data[k].norm()) + .sum::() + }) + .fold(0.0_f64, f64::max) + } } /// Diagonalizes the small `k x k` Hermitian Rayleigh-Ritz matrix, returning the smallest eigenvalue @@ -90,7 +105,10 @@ fn davidson( let diag_scale = diag.iter().map(|z| z.norm()).fold(0.0_f64, f64::max); let precondition = diag_scale > 0.0; let floor = 1e-12 * diag_scale.max(1.0); - let residual_tol = tol.sqrt(); + + // Residual gate relative to the operator scale (see the module documentation). + let anorm = op.norm_bound(); + let residual_tol = tol.max(64.0 * f64::EPSILON * anorm.max(1.0)); // Subspace basis vectors `s` and their images `A @ s`, grown one vector per cycle. let mut images: Vec> = vec![op.apply(&seed)];