diff --git a/Cargo.lock b/Cargo.lock index 84ea97efc..110a5f3c1 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -39,6 +39,25 @@ version = "0.2.21" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "683d7910e743518b0e34f1186f92494becacb047c7b6bf616c96772180fef923" +[[package]] +name = "alloy-rlp" +version = "0.3.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "24671b1f62edcf0f9b62994c7bf72cd621a04a4b99f5020ece1a647b40e2f103" +dependencies = [ + "arrayvec", + "bytes", +] + +[[package]] +name = "android_system_properties" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae221649c9976a6f6c56ae1facf410f3ddb33cc661c4b7b61020a912d4237fbc" +dependencies = [ + "libc", +] + [[package]] name = "anes" version = "0.1.6" @@ -129,12 +148,298 @@ dependencies = [ "yansi", ] +[[package]] +name = "ark-ff" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b3235cc41ee7a12aaaf2c575a2ad7b46713a8a50bda2fc3b003a04845c05dd6" +dependencies = [ + "ark-ff-asm 0.3.0", + "ark-ff-macros 0.3.0", + "ark-serialize 0.3.0", + "ark-std 0.3.0", + "derivative", + "num-bigint", + "num-traits", + "paste", + "rustc_version 0.3.3", + "zeroize", +] + +[[package]] +name = "ark-ff" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec847af850f44ad29048935519032c33da8aa03340876d351dfab5660d2966ba" +dependencies = [ + "ark-ff-asm 0.4.2", + "ark-ff-macros 0.4.2", + "ark-serialize 0.4.2", + "ark-std 0.4.0", + "derivative", + "digest 0.10.7", + "itertools 0.10.5", + "num-bigint", + "num-traits", + "paste", + "rustc_version 0.4.1", + "zeroize", +] + +[[package]] +name = "ark-ff" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a177aba0ed1e0fbb62aa9f6d0502e9b46dad8c2eab04c14258a1212d2557ea70" +dependencies = [ + "ark-ff-asm 0.5.0", + "ark-ff-macros 0.5.0", + "ark-serialize 0.5.0", + "ark-std 0.5.0", + "arrayvec", + "digest 0.10.7", + "educe", + "itertools 0.13.0", + "num-bigint", + "num-traits", + "paste", + "zeroize", +] + +[[package]] +name = "ark-ff" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f7a806ac6c8307b929df4645776290a50ee2aac754ad09d8bdf73391309e43af" +dependencies = [ + "ark-ff-asm 0.6.0", + "ark-ff-macros 0.6.0", + "ark-serialize 0.6.0", + "ark-std 0.6.0", + "digest 0.10.7", + "educe", + "num-bigint", + "num-traits", + "zeroize", +] + +[[package]] +name = "ark-ff-asm" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "db02d390bf6643fb404d3d22d31aee1c4bc4459600aef9113833d17e786c6e44" +dependencies = [ + "quote", + "syn 1.0.109", +] + +[[package]] +name = "ark-ff-asm" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ed4aa4fe255d0bc6d79373f7e31d2ea147bcf486cba1be5ba7ea85abdb92348" +dependencies = [ + "quote", + "syn 1.0.109", +] + +[[package]] +name = "ark-ff-asm" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62945a2f7e6de02a31fe400aa489f0e0f5b2502e69f95f853adb82a96c7a6b60" +dependencies = [ + "quote", + "syn 2.0.118", +] + +[[package]] +name = "ark-ff-asm" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1479009684adc073dff49a1025d3a7065b317a9ead25aaaca38cdc70058ba8a2" +dependencies = [ + "quote", + "syn 2.0.118", +] + +[[package]] +name = "ark-ff-macros" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "db2fd794a08ccb318058009eefdf15bcaaaaf6f8161eb3345f907222bac38b20" +dependencies = [ + "num-bigint", + "num-traits", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "ark-ff-macros" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7abe79b0e4288889c4574159ab790824d0033b9fdcb2a112a3182fac2e514565" +dependencies = [ + "num-bigint", + "num-traits", + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "ark-ff-macros" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09be120733ee33f7693ceaa202ca41accd5653b779563608f1234f78ae07c4b3" +dependencies = [ + "num-bigint", + "num-traits", + "proc-macro2", + "quote", + "syn 2.0.118", +] + +[[package]] +name = "ark-ff-macros" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4a0691ed21ef00ef89c1e9bda832eba493dda3ec2f8d892fb25b705f73f06bb8" +dependencies = [ + "num-bigint", + "num-traits", + "proc-macro2", + "quote", + "syn 2.0.118", +] + +[[package]] +name = "ark-serialize" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d6c2b318ee6e10f8c2853e73a83adc0ccb88995aa978d8a3408d492ab2ee671" +dependencies = [ + "ark-std 0.3.0", + "digest 0.9.0", +] + +[[package]] +name = "ark-serialize" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb7b85a02b83d2f22f89bd5cac66c9c89474240cb6207cb1efc16d098e822a5" +dependencies = [ + "ark-std 0.4.0", + "digest 0.10.7", + "num-bigint", +] + +[[package]] +name = "ark-serialize" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f4d068aaf107ebcd7dfb52bc748f8030e0fc930ac8e360146ca54c1203088f7" +dependencies = [ + "ark-std 0.5.0", + "arrayvec", + "digest 0.10.7", + "num-bigint", +] + +[[package]] +name = "ark-serialize" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a74dd304fd536fb95d0a328e72be759209cc496a9da094c5bc56e5fea4f9e86b" +dependencies = [ + "ark-serialize-derive", + "ark-std 0.6.0", + "digest 0.10.7", + "num-bigint", + "serde_with", +] + +[[package]] +name = "ark-serialize-derive" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4f153690697a2b91e5e1251ff98411ee5371500a111a0fd317a70e588eb300f9" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", +] + +[[package]] +name = "ark-std" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1df2c09229cbc5a028b1d70e00fdb2acee28b1055dfb5ca73eea49c5a25c4e7c" +dependencies = [ + "num-traits", + "rand 0.8.8", +] + +[[package]] +name = "ark-std" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94893f1e0c6eeab764ade8dc4c0db24caf4fe7cbbaafc0eba0a9030f447b5185" +dependencies = [ + "num-traits", + "rand 0.8.8", +] + +[[package]] +name = "ark-std" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "246a225cc6131e9ee4f24619af0f19d67761fff15d7ccc22e42b80846e69449a" +dependencies = [ + "num-traits", + "rand 0.8.8", +] + +[[package]] +name = "ark-std" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "367c9c827ed431bff6868b7aa926e05b16eb46603cc8b6e768e4a5553fa1d155" +dependencies = [ + "num-traits", + "rand 0.8.8", +] + +[[package]] +name = "arrayvec" +version = "0.7.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3fb67a6e08acf24fdeccbac2cb6ac4305825bd1f117462e0e6f2f193345ad56" + +[[package]] +name = "auto_impl" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ffdcb70bdbc4d478427380519163274ac86e52916e10f0a8889adf0f96d3fee7" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", +] + [[package]] name = "autocfg" version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" +[[package]] +name = "base64" +version = "0.22.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" + [[package]] name = "bincode" version = "2.0.1" @@ -226,7 +531,16 @@ dependencies = [ "proc-macro2", "quote", "rustversion", - "syn", + "syn 2.0.118", +] + +[[package]] +name = "bs58" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bf88ba1141d185c399bee5288d850d63b8369520c1eafc32a0430b5b6c287bf4" +dependencies = [ + "tinyvec", ] [[package]] @@ -235,6 +549,12 @@ version = "3.19.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "46c5e41b57b8bba42a04676d81cb89e9ee8e859a1a66f80a5a72e1cb76b34d43" +[[package]] +name = "byte-slice-cast" +version = "1.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7575182f7272186991736b70173b0ea045398f984bf5ebbb3804736ce1330c9d" + [[package]] name = "bytemuck" version = "1.25.0" @@ -247,6 +567,12 @@ version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" +[[package]] +name = "bytes" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc652a48c352aef3ea3aed32080501cf3ef6ed5da78602a020c991775b0aff04" + [[package]] name = "cassowary" version = "0.3.0" @@ -295,6 +621,18 @@ dependencies = [ "rand_core 0.10.0", ] +[[package]] +name = "chrono" +version = "0.4.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1aa79e62e7697b8e29b513a68abacf485adcd1fe8284a4316c5ae868e6633327" +dependencies = [ + "iana-time-zone", + "num-traits", + "serde", + "windows-link", +] + [[package]] name = "chumsky" version = "0.10.1" @@ -381,7 +719,7 @@ dependencies = [ "heck", "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -452,6 +790,27 @@ dependencies = [ "windows-sys 0.59.0", ] +[[package]] +name = "const_format" +version = "0.2.36" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4481a617ad9a412be3b97c5d403fef8ed023103368908b9c50af598ff467cc1e" +dependencies = [ + "const_format_proc_macros", + "konst", +] + +[[package]] +name = "const_format_proc_macros" +version = "0.2.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d57c2eccfb16dbac1f4e61e206105db5820c9d26c3c472bc17c774259ef7744" +dependencies = [ + "proc-macro2", + "quote", + "unicode-xid", +] + [[package]] name = "convert_case" version = "0.6.0" @@ -461,6 +820,12 @@ dependencies = [ "unicode-segmentation", ] +[[package]] +name = "core-foundation-sys" +version = "0.8.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" + [[package]] name = "cpufeatures" version = "0.3.0" @@ -594,6 +959,16 @@ version = "0.2.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" +[[package]] +name = "crypto-common" +version = "0.1.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78c8292055d1c1df0cce5d180393dc8cce0abec0a7102adb6c7b1eef6016d60a" +dependencies = [ + "generic-array", + "typenum", +] + [[package]] name = "darling" version = "0.21.3" @@ -625,7 +1000,7 @@ dependencies = [ "proc-macro2", "quote", "strsim", - "syn", + "syn 2.0.118", ] [[package]] @@ -638,7 +1013,7 @@ dependencies = [ "proc-macro2", "quote", "strsim", - "syn", + "syn 2.0.118", ] [[package]] @@ -649,7 +1024,7 @@ checksum = "d38308df82d1080de0afee5d069fa14b0326a88c14f15c5ccda35b4a6c414c81" dependencies = [ "darling_core 0.21.3", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -660,7 +1035,7 @@ checksum = "ac3984ec7bd6cfa798e62b4a642426a5be0e68f9401cfc2a01e3fa9ea2fcdb8d" dependencies = [ "darling_core 0.23.0", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -697,7 +1072,7 @@ dependencies = [ "defmt-parser", "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -709,6 +1084,62 @@ dependencies = [ "thiserror", ] +[[package]] +name = "deranged" +version = "0.5.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7cd812cc2bc1d69d4764bd80df88b4317eaef9e773c75226407d9bc0876b211c" +dependencies = [ + "serde_core", +] + +[[package]] +name = "derivative" +version = "2.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fcc3dd5e9e9c0b295d6e1e4d811fb6f157d5ffd784b8d202fc62eac8035a770b" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "digest" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3dd60d1080a57a05ab032377049e0591415d2b31afd7028356dbf3cc6dcb066" +dependencies = [ + "generic-array", +] + +[[package]] +name = "digest" +version = "0.10.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ed9a281f7bc9b7576e61468ba615a66a5c8cfdff42420a70aa82701a3b1e292" +dependencies = [ + "crypto-common", +] + +[[package]] +name = "dyn-clone" +version = "1.0.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d0881ea181b1df73ff77ffaaf9c7544ecc11e82fba9b5f27b262a3c73a332555" + +[[package]] +name = "educe" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d7bc049e1bd8cdeb31b68bbd586a9464ecf9f3944af3958a7a9d0f8b9799417" +dependencies = [ + "enum-ordinalize", + "proc-macro2", + "quote", + "syn 2.0.118", +] + [[package]] name = "either" version = "1.15.0" @@ -721,6 +1152,26 @@ version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "34aa73646ffb006b8f5147f3dc182bd4bcb190227ce861fc4a4844bf8e3cb2c0" +[[package]] +name = "enum-ordinalize" +version = "4.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "89dd01549b09589510cf0647475075d12071456586d70f5c75c98ae2a5537677" +dependencies = [ + "enum-ordinalize-derive", +] + +[[package]] +name = "enum-ordinalize-derive" +version = "4.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a65863d15a4ce2888bd2f0f543cc963d3879c3a022c8ee43f6141d479a3ac815" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.4", +] + [[package]] name = "env_filter" version = "2.0.0" @@ -776,12 +1227,46 @@ version = "2.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9f1f227452a390804cdb637b74a86990f2a7d7ba4b7d5693aac9b4dd6defd8d6" +[[package]] +name = "fastrlp" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "139834ddba373bbdd213dffe02c8d110508dcf1726c2be27e8d1f7d7e1856418" +dependencies = [ + "arrayvec", + "auto_impl", + "bytes", +] + +[[package]] +name = "fastrlp" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ce8dba4714ef14b8274c371879b175aa55b16b30f269663f19d576f380018dc4" +dependencies = [ + "arrayvec", + "auto_impl", + "bytes", +] + [[package]] name = "find-msvc-tools" version = "0.1.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5baebc0774151f905a1a2cc41989300b1e6fbb29aff0ceffa1064fdd3088d582" +[[package]] +name = "fixed-hash" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "835c052cb0c08c1acf6ffd71c022172e18723949c8282f2b9f27efbc51e64534" +dependencies = [ + "byteorder", + "rand 0.8.8", + "rustc-hex", + "static_assertions", +] + [[package]] name = "fnv" version = "1.0.7" @@ -809,6 +1294,27 @@ dependencies = [ "byteorder", ] +[[package]] +name = "generic-array" +version = "0.14.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85649ca51fd72272d7821adaf274ad91c288277713d9c18820d8499a7ff69e9a" +dependencies = [ + "typenum", + "version_check", +] + +[[package]] +name = "getrandom" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0" +dependencies = [ + "cfg-if", + "libc", + "wasi 0.11.1+wasi-snapshot-preview1", +] + [[package]] name = "getrandom" version = "0.3.3" @@ -856,6 +1362,12 @@ dependencies = [ "crunchy", ] +[[package]] +name = "hashbrown" +version = "0.12.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8a9ee70c43aaf417c914396645a0fa852624801b24ebb7ae78fe8272889ac888" + [[package]] name = "hashbrown" version = "0.14.5" @@ -885,6 +1397,36 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "hex" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f24254aa9a54b5c858eaee2f5bccdb46aaf0e486a595ed5fd8f86ba55232a70" + +[[package]] +name = "iana-time-zone" +version = "0.1.65" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e31bc9ad994ba00e440a8aa5c9ef0ec67d5cb5e5cb0cc7f8b744a35b389cc470" +dependencies = [ + "android_system_properties", + "core-foundation-sys", + "iana-time-zone-haiku", + "js-sys", + "log", + "wasm-bindgen", + "windows-core", +] + +[[package]] +name = "iana-time-zone-haiku" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f31827a206f56af32e590ba56d5d2d085f558508192593743f16b2306495269f" +dependencies = [ + "cc", +] + [[package]] name = "id-arena" version = "2.3.0" @@ -897,12 +1439,43 @@ version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" +[[package]] +name = "impl-codec" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba6a270039626615617f3f36d15fc827041df3b78c439da2cadfa47455a77f2f" +dependencies = [ + "parity-scale-codec", +] + +[[package]] +name = "impl-trait-for-tuples" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a0eb5a3343abf848c0984fe4604b2b105da9539376e24fc0a3b0007411ae4fd9" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", +] + [[package]] name = "indenter" version = "0.3.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "964de6e86d545b246d84badc0fef527924ace5134f30641c203ef52ba83f58d5" +[[package]] +name = "indexmap" +version = "1.9.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bd070e393353796e801d209ad339e89596eb4c8d430d18ede6a1cced8fafbd99" +dependencies = [ + "autocfg", + "hashbrown 0.12.3", + "serde", +] + [[package]] name = "indexmap" version = "2.11.4" @@ -946,7 +1519,7 @@ dependencies = [ "indoc", "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -964,6 +1537,15 @@ version = "1.70.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a6cb138bb79a146c1bd460005623e142ef0181e3d0219cb493e02f7d08a35695" +[[package]] +name = "itertools" +version = "0.10.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b0fd2260e829bddf4cb6ea802289de2f86d6a7a690192fbe91b3f46e0f2c8473" +dependencies = [ + "either", +] + [[package]] name = "itertools" version = "0.13.0" @@ -996,10 +1578,12 @@ checksum = "ccfe6121cbe750cf81efa362d85c0bde7ea298ec43092d3a193baca59cdbd634" dependencies = [ "defmt", "jiff-static", + "jiff-tzdb-platform", "log", "portable-atomic", "portable-atomic-util", "serde_core", + "windows-link", ] [[package]] @@ -1010,7 +1594,22 @@ checksum = "e165e897f662d428f3cd3828a919dbe067c2d42bb1031eede74ef9d27ecdedd2" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.118", +] + +[[package]] +name = "jiff-tzdb" +version = "0.1.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "142bd39932ad231f10513df9ab62661fead8719872150b7ad02a2df79f4e141e" + +[[package]] +name = "jiff-tzdb-platform" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "875a5a69ac2bab1a891711cf5eccbec1ce0341ea805560dcd90b7a2e925132e8" +dependencies = [ + "jiff-tzdb", ] [[package]] @@ -1023,6 +1622,21 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "konst" +version = "0.2.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "128133ed7824fcd73d6e7b17957c5eb7bacb885649bd8c69708b2331a10bcefb" +dependencies = [ + "konst_macro_rules", +] + +[[package]] +name = "konst_macro_rules" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a4933f3f57a8e9d9da04db23fb153356ecaf00cbd14aee46279c33dc80925c37" + [[package]] name = "leb128fmt" version = "0.1.0" @@ -1035,6 +1649,12 @@ version = "0.2.186" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "68ab91017fe16c622486840e4c83c9a37afeff978bd239b5293d61ece587de66" +[[package]] +name = "libm" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" + [[package]] name = "libmimalloc-sys" version = "0.1.49" @@ -1081,6 +1701,16 @@ dependencies = [ "hashbrown 0.15.5", ] +[[package]] +name = "matrixmultiply" +version = "0.3.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f607c237553f086e7043417a51df26b2eb899d3caff94e6a67592ff992fedc7" +dependencies = [ + "autocfg", + "rawpointer", +] + [[package]] name = "memchr" version = "2.7.6" @@ -1108,6 +1738,21 @@ dependencies = [ "windows-sys 0.61.2", ] +[[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" version = "0.4.3" @@ -1138,9 +1783,16 @@ version = "0.4.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" dependencies = [ + "bytemuck", "num-traits", ] +[[package]] +name = "num-conv" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "521739c6d2bac4aa25192232afe6841231376b2b26d4d9fae5ecf8ca5772e441" + [[package]] name = "num-integer" version = "0.1.46" @@ -1179,6 +1831,23 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", + "libm", +] + +[[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]] @@ -1218,6 +1887,34 @@ dependencies = [ "winapi", ] +[[package]] +name = "parity-scale-codec" +version = "3.7.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "799781ae679d79a948e13d4824a40970bfa500058d245760dd857301059810fa" +dependencies = [ + "arrayvec", + "bitvec", + "byte-slice-cast", + "const_format", + "impl-trait-for-tuples", + "parity-scale-codec-derive", + "rustversion", + "serde", +] + +[[package]] +name = "parity-scale-codec-derive" +version = "3.7.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34b4653168b563151153c9e4c08ebed57fb8262bebfa79711552fa983c623e7a" +dependencies = [ + "proc-macro-crate", + "proc-macro2", + "quote", + "syn 2.0.118", +] + [[package]] name = "parking_lot" version = "0.12.4" @@ -1247,6 +1944,16 @@ version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" +[[package]] +name = "pest" +version = "2.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a07a60cc7a4d00c91f95c685609d1d2f79050e6804b70ebedd7650f0b839bcf" +dependencies = [ + "memchr", + "ucd-trie", +] + [[package]] name = "plotters" version = "0.3.7" @@ -1290,6 +1997,12 @@ dependencies = [ "portable-atomic", ] +[[package]] +name = "powerfmt" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "439ee305def115ba05938db6eb1644ff94165c5ab5e9420d1c1bcedbba909391" + [[package]] name = "ppv-lite86" version = "0.2.21" @@ -1319,6 +2032,21 @@ dependencies = [ "ratatui", ] +[[package]] +name = "ppvm-lindblad" +version = "0.1.0" +dependencies = [ + "fxhash", + "ndarray", + "num", + "ppvm-pauli-sum", + "ppvm-pauli-word", + "ppvm-traits", + "quspin-expm", + "quspin-types", + "rayon", +] + [[package]] name = "ppvm-pauli-sum" version = "0.1.0" @@ -1329,7 +2057,7 @@ dependencies = [ "dashmap", "fxhash", "gxhash", - "indexmap", + "indexmap 2.11.4", "insta", "itertools 0.14.0", "num", @@ -1359,7 +2087,11 @@ name = "ppvm-python-native" version = "0.1.0" dependencies = [ "bnum", + "mimalloc", + "num", + "numpy", "paste", + "ppvm-lindblad", "ppvm-pauli-sum", "ppvm-stim", "ppvm-tableau", @@ -1455,7 +2187,7 @@ dependencies = [ "dashmap", "fxhash", "gxhash", - "indexmap", + "indexmap 2.11.4", "insta", "num", "rayon", @@ -1498,7 +2230,18 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "479ca8adacdd7ce8f1fb39ce9ecccbfe93a3f1344b3d0d97f20bc0196208f62b" dependencies = [ "proc-macro2", - "syn", + "syn 2.0.118", +] + +[[package]] +name = "primitive-types" +version = "0.12.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b34d9fd68ae0b74a41b21c03c2f62847aa0ffea044eee893b4c140b37e244e2" +dependencies = [ + "fixed-hash", + "impl-codec", + "uint", ] [[package]] @@ -1530,7 +2273,7 @@ dependencies = [ "bitflags 2.11.1", "num-traits", "rand 0.9.4", - "rand_chacha", + "rand_chacha 0.9.0", "rand_xorshift", "regex-syntax 0.8.11", "rusty-fork", @@ -1591,7 +2334,7 @@ dependencies = [ "proc-macro2", "pyo3-macros-backend", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -1603,7 +2346,7 @@ dependencies = [ "heck", "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -1621,6 +2364,28 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "quspin-expm" +version = "0.1.0" +source = "git+https://github.com/QuSpin/QuSpin-rust?rev=a0ad6c9fe2e8063208f9ba1c6677150c993bb554#a0ad6c9fe2e8063208f9ba1c6677150c993bb554" +dependencies = [ + "ndarray", + "num-complex", + "quspin-types", + "rayon", +] + +[[package]] +name = "quspin-types" +version = "0.1.0" +source = "git+https://github.com/QuSpin/QuSpin-rust?rev=a0ad6c9fe2e8063208f9ba1c6677150c993bb554#a0ad6c9fe2e8063208f9ba1c6677150c993bb554" +dependencies = [ + "ndarray", + "num-complex", + "ruint", + "thiserror", +] + [[package]] name = "r-efi" version = "5.3.0" @@ -1633,13 +2398,24 @@ version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "dc33ff2d4973d518d823d61aa239014831e521c75da58e3df4840d3f47749d09" +[[package]] +name = "rand" +version = "0.8.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e058c7de0b26af77780c769414d6257830bb240f3c38477dbc2c16e5f54d6d4c" +dependencies = [ + "libc", + "rand_chacha 0.3.1", + "rand_core 0.6.4", +] + [[package]] name = "rand" version = "0.9.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "44c5af06bb1b7d3216d91932aed5265164bf384dc89cd6ba05cf59a35f5f76ea" dependencies = [ - "rand_chacha", + "rand_chacha 0.9.0", "rand_core 0.9.5", ] @@ -1654,6 +2430,16 @@ dependencies = [ "rand_core 0.10.0", ] +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core 0.6.4", +] + [[package]] name = "rand_chacha" version = "0.9.0" @@ -1664,6 +2450,15 @@ dependencies = [ "rand_core 0.9.5", ] +[[package]] +name = "rand_core" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom 0.2.17", +] + [[package]] name = "rand_core" version = "0.9.5" @@ -1709,6 +2504,12 @@ dependencies = [ "unicode-width 0.2.0", ] +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + [[package]] name = "rayon" version = "1.11.0" @@ -1738,6 +2539,26 @@ dependencies = [ "bitflags 2.11.1", ] +[[package]] +name = "ref-cast" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7e440fb4e4b4147295338efb76001ab9e4efc0e5839df2c47fc5ac2381d365c3" +dependencies = [ + "ref-cast-impl", +] + +[[package]] +name = "ref-cast-impl" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "92ecd8964f8453721699a1ed72037b0db49ce2f5a5138486ee89bed6f67cdf3a" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.4", +] + [[package]] name = "regex" version = "1.12.4" @@ -1784,6 +2605,81 @@ version = "0.8.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" +[[package]] +name = "rlp" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb919243f34364b6bd2fc10ef797edbfa75f33c252e7998527479c6d6b47e1ec" +dependencies = [ + "bytes", + "rustc-hex", +] + +[[package]] +name = "ruint" +version = "1.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f5e99bff0393163bb25029a6af25d3d8d202ba5b5438a74d1bd8789f5c822970" +dependencies = [ + "alloy-rlp", + "ark-ff 0.3.0", + "ark-ff 0.4.2", + "ark-ff 0.5.0", + "ark-ff 0.6.0", + "bytes", + "fastrlp 0.3.1", + "fastrlp 0.4.0", + "num-bigint", + "num-integer", + "num-traits", + "parity-scale-codec", + "primitive-types", + "proptest", + "rand 0.8.8", + "rand 0.9.4", + "rlp", + "ruint-macro", + "serde_core", + "valuable", + "zeroize", +] + +[[package]] +name = "ruint-macro" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "48fd7bd8a6377e15ad9d42a8ec25371b94ddc67abe7c8b9127bec79bebaaae18" + +[[package]] +name = "rustc-hash" +version = "2.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b1e7f9a428571be2dc5bc0505c13fb6bf936822b894ec87abf8a08a4e51742d" + +[[package]] +name = "rustc-hex" +version = "2.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3e75f6a532d0fd9f7f13144f392b6ad56a32696bfcd9c78f797f16bbb6f072d6" + +[[package]] +name = "rustc_version" +version = "0.3.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0dfe2087c51c460008730de8b57e6a320782fbfb312e1f4d520e6c6fae155ee" +dependencies = [ + "semver 0.11.0", +] + +[[package]] +name = "rustc_version" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cfcb3a22ef46e85b45de6ee7e79d063319ebb6594faafcf1c225ea92ab6e9b92" +dependencies = [ + "semver 1.0.27", +] + [[package]] name = "rustix" version = "0.38.44" @@ -1843,18 +2739,60 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "schemars" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4cd191f9397d57d581cddd31014772520aa448f65ef991055d7f61582c65165f" +dependencies = [ + "dyn-clone", + "ref-cast", + "serde", + "serde_json", +] + +[[package]] +name = "schemars" +version = "1.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "687274d293b6cdc6e73e0fee520bf2049650090d7164f87672d212a3c530cf4a" +dependencies = [ + "dyn-clone", + "ref-cast", + "serde", + "serde_json", +] + [[package]] name = "scopeguard" version = "1.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" +[[package]] +name = "semver" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f301af10236f6df4160f7c3f04eec6dbc70ace82d23326abad5edee88801c6b6" +dependencies = [ + "semver-parser", +] + [[package]] name = "semver" version = "1.0.27" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d767eb0aabc880b29956c35734170f26ed551a859dbd361d140cdbeca61ab1e2" +[[package]] +name = "semver-parser" +version = "0.10.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9900206b54a3527fdc7b8a938bffd94a568bac4f4aa8113b209df75a09c0dec2" +dependencies = [ + "pest", +] + [[package]] name = "serde" version = "1.0.228" @@ -1882,7 +2820,7 @@ checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -1898,6 +2836,26 @@ dependencies = [ "serde_core", ] +[[package]] +name = "serde_with" +version = "3.22.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ee78f1fbe43ac4a0e47aadb3dbd357b69eb0d3793e948624cd03dd2750ab1c0a" +dependencies = [ + "base64", + "bs58", + "chrono", + "hex", + "indexmap 1.9.3", + "indexmap 2.11.4", + "jiff", + "schemars 0.9.0", + "schemars 1.2.2", + "serde_core", + "serde_json", + "time", +] + [[package]] name = "shlex" version = "1.3.0" @@ -2000,7 +2958,18 @@ dependencies = [ "proc-macro2", "quote", "rustversion", - "syn", + "syn 2.0.118", +] + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", ] [[package]] @@ -2014,6 +2983,17 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "syn" +version = "3.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6275cddf4610d1775e6d1fe9469b2e77d0f39fd98fb7450901b821e0c53649f" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + [[package]] name = "tap" version = "1.0.1" @@ -2065,7 +3045,37 @@ checksum = "ebc4ee7f67670e9b64d05fa4253e753e016c6c95ff35b89b7941d6b856dec1d5" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.118", +] + +[[package]] +name = "time" +version = "0.3.55" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cdb87b95ec50ddfa440816d227a17b2ccbdda963a316a727fda0fc4334f7d134" +dependencies = [ + "deranged", + "num-conv", + "powerfmt", + "serde_core", + "time-core", + "time-macros", +] + +[[package]] +name = "time-core" +version = "0.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e1c906769ad99c88eaa54e728060edef082f8e358ff32030cb7c7d315e81109" + +[[package]] +name = "time-macros" +version = "0.2.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7e689342a48d2ea927c87ea50cabf8594854bf940e9310208848d680d668ed85" +dependencies = [ + "num-conv", + "time-core", ] [[package]] @@ -2078,6 +3088,21 @@ dependencies = [ "serde_json", ] +[[package]] +name = "tinyvec" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb4ebadaa0af04fab11ae01eb5f9fdb5f9c5b875506e210e71c07873528baa7f" +dependencies = [ + "tinyvec_macros", +] + +[[package]] +name = "tinyvec_macros" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20" + [[package]] name = "toml_datetime" version = "1.1.1+spec-1.1.0" @@ -2093,7 +3118,7 @@ version = "0.25.6+spec-1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0db3bae107c9522f86d361697dee1d7386a2ddcf659d5aea5159819a21a3c4a7" dependencies = [ - "indexmap", + "indexmap 2.11.4", "toml_datetime", "toml_parser", "winnow", @@ -2108,6 +3133,30 @@ dependencies = [ "winnow", ] +[[package]] +name = "typenum" +version = "1.20.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6f5e870be6c3b371b77fe0ee0bafb859fa4964b4404c27de1d380043c4dda20" + +[[package]] +name = "ucd-trie" +version = "0.1.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2896d95c02a80c6d6a5d6e953d479f5ddf2dfdb6a244441010e373ac0fb88971" + +[[package]] +name = "uint" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76f64bba2c53b04fcab63c01a7d7427eadc821e3bc48c34dc9ba29c501164b52" +dependencies = [ + "byteorder", + "crunchy", + "hex", + "static_assertions", +] + [[package]] name = "unarray" version = "0.1.4" @@ -2167,6 +3216,12 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" +[[package]] +name = "valuable" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba73ea9cf16a25df0c8caa16c51acb937d5712a8429db78a3ee29d5dcacd3a65" + [[package]] name = "version_check" version = "0.9.5" @@ -2219,7 +3274,7 @@ dependencies = [ "proc-macro-crate", "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -2293,7 +3348,7 @@ dependencies = [ "eyre", "proc-macro2", "quote", - "syn", + "syn 2.0.118", "vihaco-parser", ] @@ -2322,7 +3377,7 @@ dependencies = [ "proc-macro-crate", "proc-macro2", "quote", - "syn", + "syn 2.0.118", ] [[package]] @@ -2428,7 +3483,7 @@ dependencies = [ "log", "proc-macro2", "quote", - "syn", + "syn 2.0.118", "wasm-bindgen-shared", ] @@ -2450,7 +3505,7 @@ checksum = "9f07d2f20d4da7b26400c9f4a0511e6e0345b040694e8a75bd41d578fa4421d7" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.118", "wasm-bindgen-backend", "wasm-bindgen-shared", ] @@ -2481,7 +3536,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bb0e353e6a2fbdc176932bbaab493762eb1255a7900fe0fea1a2f96c296cc909" dependencies = [ "anyhow", - "indexmap", + "indexmap 2.11.4", "wasm-encoder", "wasmparser", ] @@ -2494,8 +3549,8 @@ checksum = "47b807c72e1bac69382b3a6fb3dbe8ea4c0ed87ff5629b8685ae6b9a611028fe" dependencies = [ "bitflags 2.11.1", "hashbrown 0.15.5", - "indexmap", - "semver", + "indexmap 2.11.4", + "semver 1.0.27", ] [[package]] @@ -2539,12 +3594,65 @@ version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" +[[package]] +name = "windows-core" +version = "0.62.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b8e83a14d34d0623b51dce9581199302a221863196a1dde71a7663a4c2be9deb" +dependencies = [ + "windows-implement", + "windows-interface", + "windows-link", + "windows-result", + "windows-strings", +] + +[[package]] +name = "windows-implement" +version = "0.60.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "053e2e040ab57b9dc951b72c264860db7eb3b0200ba345b4e4c3b14f67855ddf" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", +] + +[[package]] +name = "windows-interface" +version = "0.59.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f316c4a2570ba26bbec722032c4099d8c8bc095efccdc15688708623367e358" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", +] + [[package]] name = "windows-link" version = "0.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" +[[package]] +name = "windows-result" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7781fa89eaf60850ac3d2da7af8e5242a5ea78d1a11c49bf2910bb5a73853eb5" +dependencies = [ + "windows-link", +] + +[[package]] +name = "windows-strings" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7837d08f69c77cf6b07689544538e017c1bfcf57e34b4c0ff58e6c2cd3b37091" +dependencies = [ + "windows-link", +] + [[package]] name = "windows-sys" version = "0.59.0" @@ -2670,9 +3778,9 @@ checksum = "b7c566e0f4b284dd6561c786d9cb0142da491f46a9fbed79ea69cdad5db17f21" dependencies = [ "anyhow", "heck", - "indexmap", + "indexmap 2.11.4", "prettyplease", - "syn", + "syn 2.0.118", "wasm-metadata", "wit-bindgen-core", "wit-component", @@ -2688,7 +3796,7 @@ dependencies = [ "prettyplease", "proc-macro2", "quote", - "syn", + "syn 2.0.118", "wit-bindgen-core", "wit-bindgen-rust", ] @@ -2701,7 +3809,7 @@ checksum = "9d66ea20e9553b30172b5e831994e35fbde2d165325bec84fc43dbf6f4eb9cb2" dependencies = [ "anyhow", "bitflags 2.11.1", - "indexmap", + "indexmap 2.11.4", "log", "serde", "serde_derive", @@ -2720,9 +3828,9 @@ checksum = "ecc8ac4bc1dc3381b7f59c34f00b67e18f910c2c0f50015669dde7def656a736" dependencies = [ "anyhow", "id-arena", - "indexmap", + "indexmap 2.11.4", "log", - "semver", + "semver 1.0.27", "serde", "serde_derive", "serde_json", @@ -2762,5 +3870,25 @@ checksum = "88d2b8d9c68ad2b9e4340d7832716a4d21a22a1154777ad56ea55c51a9cf3831" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.118", +] + +[[package]] +name = "zeroize" +version = "1.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e13c156562582aa81c60cb29407084cdb54c4164760106ab78e6c5b0858cf64e" +dependencies = [ + "zeroize_derive", +] + +[[package]] +name = "zeroize_derive" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c50655cbb0fe3fc43170059e702f1ce5e19b84cec58dc87b037a09935c2f328" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", ] diff --git a/Cargo.toml b/Cargo.toml index 35cb8102a..b957023b0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -14,6 +14,7 @@ members = [ "crates/ppvm-pauli-word", "crates/ppvm-pauli-sum", "crates/ppvm-sym", + "crates/ppvm-lindblad", "crates/ppvm-python-native", "crates/ppvm-tableau", "crates/ppvm-stim", "crates/stim-parser", diff --git a/crates/ppvm-lindblad/Cargo.toml b/crates/ppvm-lindblad/Cargo.toml new file mode 100644 index 000000000..69f425341 --- /dev/null +++ b/crates/ppvm-lindblad/Cargo.toml @@ -0,0 +1,23 @@ +[package] +name = "ppvm-lindblad" +version = "0.1.0" +edition = "2024" +description = "Direct Heisenberg-picture Lindbladian evolution on an adaptive Pauli-string basis." + +[dependencies] +fxhash = "0.2.1" +ndarray = "0.17" +num = "0.4.3" +ppvm-traits = { version = "0.1.0", path = "../ppvm-traits" } +ppvm-pauli-word = { version = "0.1.0", path = "../ppvm-pauli-word" } +rayon = "1.11" +# Matrix-exponential action (Al-Mohy & Higham). QuSpin-rust is MIT-licensed; +# the pinned rev is the commit that added the LICENSE file. +quspin-expm = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2e8063208f9ba1c6677150c993bb554" } +# `QuSpinError` (the error type returned by the `LinearOperator` trait methods +# we implement in `mf_expm.rs`) is not re-exported from `quspin-expm`'s root, +# so we depend on `quspin-types` directly. Same git rev as `quspin-expm`. +quspin-types = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2e8063208f9ba1c6677150c993bb554" } + +[dev-dependencies] +ppvm-pauli-sum = { version = "0.1.0", path = "../ppvm-pauli-sum" } diff --git a/crates/ppvm-lindblad/src/algebra.rs b/crates/ppvm-lindblad/src/algebra.rs new file mode 100644 index 000000000..da44c3042 --- /dev/null +++ b/crates/ppvm-lindblad/src/algebra.rs @@ -0,0 +1,78 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Hot-path Pauli product / commutator on raw word chunks. +//! +//! Phase encoding: Pauli product `P·Q = ω · R` where `ω = i^phase` and +//! `phase ∈ {0,1,2,3}` ↔ `ω ∈ {1, i, -1, -i}`. The per-byte XOR/AND +//! formulas are the same ones used by +//! [`ppvm_pauli_word::phase::PhasedPauliWord`]'s `MulAssign`. This module +//! keeps a copy that returns the unpacked `(word, phase)` pair without +//! constructing a phased wrapper. + +use crate::word::{W_CHUNKS, Word}; +use num::Complex; +use ppvm_traits::PauliWordTrait; + +#[inline(always)] +pub(crate) fn phase_factor(phase: u8) -> Complex { + match phase & 3 { + 0 => Complex::new(1.0, 0.0), + 1 => Complex::new(0.0, 1.0), + 2 => Complex::new(-1.0, 0.0), + _ => Complex::new(0.0, -1.0), + } +} + +/// `true` if Pauli words `a` and `b` anti-commute. +/// +/// Two Pauli strings anti-commute iff +/// `popcount(a.x & b.z) + popcount(a.z & b.x)` is odd. +#[inline(always)] +pub(crate) fn anti_commutes(a: &Word, b: &Word) -> bool { + let mut bits: u32 = 0; + for i in 0..W_CHUNKS { + bits += (a.xbits.data[i] & b.zbits.data[i]).count_ones(); + bits += (a.zbits.data[i] & b.xbits.data[i]).count_ones(); + } + bits & 1 == 1 +} + +/// Commutator product `h · p`: returns `(out, eps)` where `out = h ⊕ p` and +/// +/// - `eps = 0` if `h` and `p` commute (caller should skip — `[h,p] = 0`), +/// - `eps = -2.0` if `h·p` has phase `+i` (so `i·[h,p] = -2·out`), +/// - `eps = +2.0` if `h·p` has phase `-i` (so `i·[h,p] = +2·out`). +#[inline(always)] +pub(crate) fn comm_product(h: &Word, p: &Word) -> (Word, f64) { + let (out, phase) = pauli_mul(h, p); + let eps = match phase { + 1 => -2.0, + 3 => 2.0, + _ => 0.0, + }; + (out, eps) +} + +/// Full Pauli product `p · q`: returns `(out, phase)` where the product +/// is `ω · out` with `ω = i^phase`. +#[inline(always)] +pub(crate) fn pauli_mul(p: &Word, q: &Word) -> (Word, u8) { + let mut out = Word::new(p.n_qubits()); + let mut sign_count: u32 = 0; + let mut imag_count: u32 = 0; + for i in 0..W_CHUNKS { + let a = p.xbits.data[i]; + let b = p.zbits.data[i]; + let c = q.xbits.data[i]; + let d = q.zbits.data[i]; + let sign = (a & b & c & !d) | (a & !b & !c & d) | (!a & b & c & d); + let imag = (a & !b & d) | (a & !c & d) | (!a & b & c) | (b & c & !d); + sign_count += sign.count_ones(); + imag_count += imag.count_ones(); + out.xbits.data[i] = a ^ c; + out.zbits.data[i] = b ^ d; + } + out.rehash(); + (out, ((2 * sign_count + imag_count) & 3) as u8) +} diff --git a/crates/ppvm-lindblad/src/basis.rs b/crates/ppvm-lindblad/src/basis.rs new file mode 100644 index 000000000..9f352d972 --- /dev/null +++ b/crates/ppvm-lindblad/src/basis.rs @@ -0,0 +1,256 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Basis-level `L*` operators: in-basis generator and off-basis leakage. + +use crate::Error; +use crate::spec::LindbladSpec; +use crate::word::{Word, word_hash}; +use fxhash::{FxBuildHasher, FxHashMap}; +use num::Complex; +use rayon::prelude::*; + +/// Build a `word → row` map for a basis assumed to contain unique Pauli +/// words; debug-asserts the uniqueness invariant. +pub fn build_basis_index(basis: &[Word]) -> FxHashMap { + let mut index: FxHashMap = FxHashMap::default(); + for (i, w) in basis.iter().enumerate() { + let prev = index.insert(*w, i as u32); + debug_assert!( + prev.is_none(), + "basis contains duplicate Pauli word at positions {} and {}", + prev.unwrap(), + i, + ); + } + index +} + +impl LindbladSpec { + /// Off-basis component of `L*( Σ_j coeffs[j] · basis[j] )`. Output + /// strings that lie in `basis` or in `protected` are dropped. + pub fn leakage( + &self, + basis: &[Word], + coeffs: &[f64], + protected: &[Word], + ) -> Result, Error> { + self.leakage_with_prune(basis, coeffs, protected, usize::MAX, 0.0) + } + + /// Like [`Self::leakage`], but caps the live off-basis leakage map to + /// the *available room* `room = max_basis − basis.len()` — only the + /// strings we could actually add to the basis are worth keeping. The + /// cap is applied during accumulation (after each chunk), keeping the + /// `room` largest-magnitude entries. + /// + /// Basis indices are processed in descending-`|c|` order so the + /// running cap keeps the entries that are most likely to be the true + /// largest contributors. When `max_basis` is large enough that + /// `room ≥ all candidates`, nothing is dropped — the near-exact case. + pub fn leakage_with_prune( + &self, + basis: &[Word], + coeffs: &[f64], + protected: &[Word], + max_basis: usize, + tau_add: f64, + ) -> Result, Error> { + if basis.len() != coeffs.len() { + return Err(Error::LengthMismatch { + what: "basis and coeffs", + a: basis.len(), + b: coeffs.len(), + }); + } + // Hash-only membership tables: storing 8-byte `u64` keys instead + // of 48-byte Words shrinks the in-basis structure ~6×, keeping it + // in L3 (and often L2) at basis sizes where the full-Word version + // would spill to DRAM. + let in_basis: FxHashMap = basis.iter().map(|w| (word_hash(w), ())).collect(); + let protected_set: FxHashMap = + protected.iter().map(|w| (word_hash(w), ())).collect(); + + // Descending sort by |c|: process largest-magnitude contributors + // first so the running room-cap keeps the right entries. + let mut order: Vec = (0..basis.len()).collect(); + order.sort_by(|&a, &b| { + coeffs[b] + .abs() + .partial_cmp(&coeffs[a].abs()) + .unwrap_or(std::cmp::Ordering::Equal) + }); + + const CHUNK_SIZE: usize = 4096; + let room = max_basis.saturating_sub(basis.len()); + let n_qubits = self.n_qubits(); + let mut merged: FxHashMap = FxHashMap::default(); + for chunk_indices in order.chunks(CHUNK_SIZE) { + let local: Vec> = chunk_indices + .par_iter() + .map_init( + || { + ( + Vec::::with_capacity(n_qubits), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), &i| { + let p = &basis[i]; + let c = coeffs[i]; + let terms = self.compute_action_terms(p, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + for (w, v) in terms.iter() { + let h = word_hash(w); + if !in_basis.contains_key(&h) && !protected_set.contains_key(&h) { + out.push((*w, c * *v)); + } + } + out + }, + ) + .collect(); + for v in local { + for (k, val) in v { + *merged.entry(k).or_insert(0.0) += val; + } + } + + // Room-cap: keep only the `room` largest-magnitude entries. + if merged.len() > room { + if room == 0 { + merged.clear(); + } else { + let mut mags: Vec = merged.values().map(|v| v.abs()).collect(); + let k = room.min(mags.len() - 1); + mags.select_nth_unstable_by(k, |a, b| { + b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal) + }); + let cutoff = mags[k]; + merged.retain(|_, &mut v| v.abs() >= cutoff); + } + } + } + // Rate-based admission: keep only candidates whose leakage rate + // exceeds `tau_add`. `tau_add = 0` admits everything except exact + // zeros. + Ok(merged + .into_iter() + .filter(|(_, c)| c.abs() > tau_add) + .collect()) + } + + /// Sparse generator matrix in COO form: returns `(row, col, val)` + /// triplets. Row = output Pauli's position in `basis`; col = input + /// Pauli's position. Output Paulis not in `basis` are silently dropped. + /// + /// Precondition: `basis` must not contain duplicate Pauli words + /// (asserted in debug builds). + pub fn generator(&self, basis: &[Word]) -> Vec<(usize, usize, f64)> { + let index = build_basis_index(basis); + let n_qubits = self.n_qubits(); + + // `compute_action_terms` returns a deduplicated `Vec<(Word, f64)>`, + // so it can be scattered directly into COO triplets. + let local: Vec> = basis + .par_iter() + .enumerate() + .map_init( + || { + ( + Vec::::with_capacity(n_qubits), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), (col, p)| { + let terms = self.compute_action_terms(p, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + for (w, v) in terms.iter() { + if let Some(&row) = index.get(w) { + out.push((row as usize, col, *v)); + } + } + out + }, + ) + .collect(); + + // Pre-allocate the flat output to avoid sequential push reallocation. + let total: usize = local.iter().map(|v| v.len()).sum(); + let mut flat = Vec::with_capacity(total); + for v in local { + flat.extend(v); + } + flat + } + + /// Complex-coefficient variant of [`Self::leakage`]: off-basis + /// component of `L*( Σ_j coeffs[j] · basis[j] )` with complex `coeffs`. + pub fn leakage_complex( + &self, + basis: &[Word], + coeffs: &[Complex], + protected: &[Word], + ) -> Result)>, Error> { + if basis.len() != coeffs.len() { + return Err(Error::LengthMismatch { + what: "basis and coeffs", + a: basis.len(), + b: coeffs.len(), + }); + } + let in_basis: FxHashMap = basis.iter().map(|w| (word_hash(w), ())).collect(); + let protected_set: FxHashMap = + protected.iter().map(|w| (word_hash(w), ())).collect(); + + const CHUNK_SIZE: usize = 4096; + let n_qubits = self.n_qubits(); + let mut merged: FxHashMap> = FxHashMap::default(); + for chunk_start in (0..basis.len()).step_by(CHUNK_SIZE) { + let chunk_end = (chunk_start + CHUNK_SIZE).min(basis.len()); + let chunk_basis = &basis[chunk_start..chunk_end]; + let chunk_coeffs = &coeffs[chunk_start..chunk_end]; + let local: Vec)>> = chunk_basis + .par_iter() + .zip(chunk_coeffs.par_iter()) + .map_init( + || { + ( + Vec::::with_capacity(n_qubits), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), (p, &c)| { + let terms = self.compute_action_terms(p, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + for (w, v) in terms.iter() { + let h = word_hash(w); + if !in_basis.contains_key(&h) && !protected_set.contains_key(&h) { + out.push((*w, c * *v)); + } + } + out + }, + ) + .collect(); + for v in local { + for (k, val) in v { + *merged.entry(k).or_insert(Complex::new(0.0, 0.0)) += val; + } + } + } + Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect()) + } +} diff --git a/crates/ppvm-lindblad/src/config.rs b/crates/ppvm-lindblad/src/config.rs new file mode 100644 index 000000000..3fe32a71d --- /dev/null +++ b/crates/ppvm-lindblad/src/config.rs @@ -0,0 +1,59 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Configuration objects for the predictor-corrector stepper. + +/// Truncation and execution policy for a single predictor-corrector step +/// ([`crate::LindbladSpec::pc_step`], [`crate::LindbladSpec::pc_step_timed`]). +/// +/// These are the per-run *tuning knobs*, kept separate from the per-call data +/// (`basis`, `coeffs`, `dt`, `protected`). +/// +/// `max_basis` is the primary accuracy/cost dial; `admit_basis` selects the +/// displacement scheme; `drop_tol` is the churn valve of the admission-bound +/// scheme; `tau_add` is a wall optimization at most. +#[derive(Debug, Clone, Copy)] +pub struct PcStepConfig { + /// Hard rank cap on the retained basis: after the corrector, only the + /// top-`max_basis` strings by `|coeff|` are kept (protected words always + /// survive). The primary convergence dial — verify by re-running at 2×. + pub max_basis: usize, + /// Working-set (admission) bound. When `Some(a)` with `a > max_basis`, + /// enrichment may grow the live basis to `a` and the final cap performs a + /// genuine top-`max_basis`-of-union rank displacement (the analog of + /// two-site TDVP truncation at `χ_max`); `drop_tol` is then not needed + /// for membership turnover. `None` bounds admission by `max_basis` + /// itself — the valve scheme, which requires `drop_tol > 0` to keep the + /// basis adapting once it fills. + pub admit_basis: Option, + /// Magnitude prune applied after the corrector: basis entries whose + /// `|coeff|` is below `drop_tol` are discarded (protected words are + /// always kept). `<= 0.0` disables pruning — valid only with + /// `admit_basis` set, otherwise the basis freezes once it fills the cap. + pub drop_tol: f64, + /// Optional absolute rate threshold on leakage admission: a candidate is + /// admitted only if its inflow rate exceeds `tau_add`. This is the + /// natural (dt- and drop_tol-independent) parameterization — the + /// admission accuracy cliff sits at a fixed `tau_add`. `None` = no + /// filter, the recommended default with cap-based truncation. + pub tau_add: Option, + /// When `Some(n)`, run the entire step inside a freshly built rayon + /// thread pool of `n` threads (useful for benchmarking parallel + /// scaling). When `None`, the global rayon pool is used. + pub num_threads: Option, +} + +impl Default for PcStepConfig { + /// Uncapped, unfiltered, no pruning: the near-exact reference + /// configuration. Production runs should set `max_basis` (and usually + /// `admit_basis ≈ 2-3×` it). + fn default() -> Self { + Self { + max_basis: usize::MAX, + admit_basis: None, + drop_tol: 0.0, + tau_add: None, + num_threads: None, + } + } +} diff --git a/crates/ppvm-lindblad/src/error.rs b/crates/ppvm-lindblad/src/error.rs new file mode 100644 index 000000000..e0c1eb5f0 --- /dev/null +++ b/crates/ppvm-lindblad/src/error.rs @@ -0,0 +1,76 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Error type for [`crate::LindbladSpec`] construction and stepping. + +use crate::MAX_QUBITS; +use std::fmt; + +/// Errors raised when constructing a [`crate::LindbladSpec`]. +#[derive(Debug, Clone)] +pub enum Error { + TooManyQubits { + got: usize, + }, + LengthMismatch { + what: &'static str, + a: usize, + b: usize, + }, + InvalidPauliCode { + code: u8, + }, + InvalidPauliChar { + c: char, + }, + WrongLength { + expected: usize, + got: usize, + }, + NegativeRate { + index: usize, + rate: f64, + }, + EmptyLincomb { + index: usize, + }, + Internal(String), +} + +impl fmt::Display for Error { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Error::TooManyQubits { got } => { + write!( + f, + "LindbladSpec supports n_qubits ≤ {MAX_QUBITS}; got {got}" + ) + } + Error::LengthMismatch { what, a, b } => { + write!(f, "{what}: expected matching lengths, got {a} and {b}") + } + Error::InvalidPauliCode { code } => write!( + f, + "Pauli code must be 0 (I), 1 (X), 2 (Z), or 3 (Y); got {code}" + ), + Error::InvalidPauliChar { c } => { + write!(f, "invalid Pauli character '{c}'; expected I, X, Y, or Z") + } + Error::WrongLength { expected, got } => { + write!(f, "Pauli string has length {got} but n_qubits = {expected}") + } + Error::NegativeRate { index, rate } => { + write!(f, "jump rate must be non-negative; got γ_{index} = {rate}") + } + Error::EmptyLincomb { index } => { + write!( + f, + "jump {index}: lincomb must contain at least one Pauli term" + ) + } + Error::Internal(msg) => write!(f, "internal error: {msg}"), + } + } +} + +impl std::error::Error for Error {} diff --git a/crates/ppvm-lindblad/src/expm.rs b/crates/ppvm-lindblad/src/expm.rs new file mode 100644 index 000000000..f203f4f56 --- /dev/null +++ b/crates/ppvm-lindblad/src/expm.rs @@ -0,0 +1,145 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Taylor-partition selection for the `quspin-expm`-backed `exp(t·A)·b` +//! engine (driven from [`crate::mf_expm`]): the `(m, s)` selection tables +//! [`THETA`] / [`THETA_LOOSE`] from Al-Mohy & Higham (2011), used to pick +//! the partition handed to `quspin-expm`'s `from_parts`. + +/// `θ_m` table from Al-Mohy & Higham (2011), Table A.3, for double +/// precision (unit roundoff `u = 2^{-53}`). +/// +/// `θ_m` bounds `‖A‖₁` such that the degree-`m` Taylor polynomial +/// approximates `exp(A)` to within `u`. We pick `(m, s)` with +/// `s ≥ ⌈‖tA‖₁ / θ_m⌉` and minimise `m·s` (total SpMV count). +pub(crate) const THETA: &[(u32, f64)] = &[ + (1, 2.29e-16), + (2, 2.58e-8), + (3, 1.39e-5), + (4, 3.40e-4), + (5, 2.40e-3), + (6, 9.07e-3), + (7, 2.38e-2), + (8, 5.00e-2), + (9, 8.96e-2), + (10, 1.44e-1), + (11, 2.14e-1), + (12, 3.00e-1), + (13, 4.00e-1), + (14, 5.14e-1), + (15, 6.41e-1), + (16, 7.81e-1), + (17, 9.31e-1), + (18, 1.09), + (19, 1.26), + (20, 1.44), + (21, 1.62), + (22, 1.82), + (23, 2.01), + (24, 2.22), + (25, 2.43), + (26, 2.64), + (27, 2.86), + (28, 3.08), + (29, 3.31), + (30, 3.54), +]; + +/// `θ_m` table for a relaxed backward-error tolerance `tol = 1e-6`, computed +/// with the same Al-Mohy & Higham (2011) construction as [`THETA`] (the +/// backward-error series `h_{m+1}(x) = log(e^{-x} T_m(x))`; validated by +/// reproducing the `u = 2^{-53}` table above to ~2 significant figures). +/// +/// The predictor-corrector truncates the Pauli basis at `drop_tol` (typically +/// 1e-3), so computing `exp` to double-precision backward error (~1e-16) is +/// ~10 orders more accurate than the state it acts on. Using `tol = 1e-6` +/// (still ~1000x tighter than the truncation) admits a lower-degree Taylor +/// polynomial for the same `‖tA‖`, cutting the SpMV count (e.g. 23 -> 13 at +/// `‖tA‖ ≈ 2`) with no measurable effect on the truncated result. +pub(crate) const THETA_LOOSE: &[(u32, f64)] = &[ + (1, 2.000e-06), + (2, 2.447e-03), + (3, 2.863e-02), + (4, 1.025e-01), + (5, 2.262e-01), + (6, 3.911e-01), + (7, 5.866e-01), + (8, 8.045e-01), + (9, 1.039), + (10, 1.285), + (11, 1.539), + (12, 1.801), + (13, 2.067), + (14, 2.337), + (15, 2.610), + (16, 2.885), + (17, 3.162), + (18, 3.441), + (19, 3.721), + (20, 4.001), + (21, 4.282), + (22, 4.564), + (23, 4.847), + (24, 5.129), + (25, 5.412), + (26, 5.696), + (27, 5.979), + (28, 6.263), + (29, 6.546), + (30, 6.830), +]; + +/// Pick `(m, s)` minimising `s·m` subject to `s ≥ ⌈t_norm / θ_m⌉, s ≥ 1`, +/// using the `θ_m` table `theta`. Restricted to the table's `m` range; for +/// larger norms `s` simply grows linearly. +fn select_ms_with(t_norm: f64, theta: &[(u32, f64)]) -> (u32, u32) { + if t_norm <= 0.0 { + return (1, 1); + } + let mut best_m = 1u32; + let mut best_s = 1u32; + let mut best_cost = u64::MAX; + for &(m, th) in theta { + let s_f = (t_norm / th).ceil(); + let s = if s_f >= 1.0 { s_f as u32 } else { 1 }; + let cost = (m as u64) * (s as u64); + if cost < best_cost { + best_cost = cost; + best_m = m; + best_s = s; + } + } + (best_m, best_s) +} + +/// `(m, s)` selection at double-precision backward error ([`THETA`]). +pub(crate) fn select_ms(t_norm: f64) -> (u32, u32) { + select_ms_with(t_norm, THETA) +} + +/// `(m, s)` selection at the relaxed `tol = 1e-6` backward error +/// ([`THETA_LOOSE`]) — fewer SpMVs, used on the truncated PC expm path. +pub(crate) fn select_ms_loose(t_norm: f64) -> (u32, u32) { + select_ms_with(t_norm, THETA_LOOSE) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn ms_selection_sane() { + // tiny norm → small m, s = 1 + let (m, s) = select_ms(1e-9); + assert!(m <= 5, "expected small m for tiny norm, got m={m}"); + assert_eq!(s, 1); + + // moderate norm → m·s should be ~10-50 + let (m, s) = select_ms(1.0); + assert!((m * s) <= 50, "moderate norm cost too high: m={m} s={s}"); + + // large norm → s grows + let (_m, s) = select_ms(100.0); + assert!(s >= 20, "large norm should require many steps, got s={s}"); + } +} diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs new file mode 100644 index 000000000..da5414e66 --- /dev/null +++ b/crates/ppvm-lindblad/src/lib.rs @@ -0,0 +1,56 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Direct Heisenberg-picture Lindbladian evolution on an adaptive +//! Pauli-string basis. +//! +//! For a Hermitian Pauli Hamiltonian `H = Σ c_i P_i` and jump operators +//! `L_k = Σ_a λ_{k,a} P_{k,a}` (each a Hermitian-Pauli linear combination +//! with possibly complex coefficients) with rates `γ_k ≥ 0`, the adjoint +//! Lindbladian acts on a single Pauli string `p` as +//! +//! ```text +//! L*(p) = i [H, p] + Σ_k γ_k ( L_k† p L_k − 1/2 {L_k† L_k, p} ). +//! ``` +//! +//! Two jump shapes are supported with separate code paths: +//! +//! - **Hermitian Pauli** (`L = P`, `λ ∈ ℝ`): the dissipator collapses to a +//! diagonal `-2γ` on Pauli strings that anti-commute with `P`. Same fast +//! path used by every dephasing-style model. +//! +//! - **General** (complex `λ_a`, e.g. `σ± = (X ± iY)/2`): the dissipator +//! becomes a double sum `Σ_{a,b} λ_a* λ_b P_a p P_b` plus a Pauli- +//! linear-combination anti-commutator with `L†L`, which is precomputed +//! once at construction. Intermediate coefficients are complex; the +//! result is real because `L*` preserves Hermiticity, so we cast back +//! to `f64` at the boundary (with a debug-only check that `|Im|` is at +//! FP noise). +//! +//! Pauli strings are stored as [`ppvm_pauli_word::word::PauliWord`] backed by +//! two 64-bit chunks (≤128 qubits; four 32-bit chunks on 32-bit targets) +//! with cached hashes for fast HashMap lookup. The hot-path commutator/ +//! product loops bypass the higher-level word API and operate directly on +//! the raw chunks for speed. + +mod algebra; +mod basis; +pub mod config; +pub mod error; +pub(crate) mod expm; +mod spec; +mod step; +mod word; + +/// Matrix-free / quspin-expm-backed `exp(dt·L*)·b` engine. See module docs. +pub(crate) mod mf_expm; + +pub use basis::build_basis_index; +pub use config::PcStepConfig; +pub use error::Error; +pub use spec::{JumpInput, LindbladSpec}; +pub use step::PcStepTimings; +pub use word::{MAX_QUBITS, Word, codes_from_word, parse_pauli_string, word_from_codes}; + +#[cfg(test)] +mod tests; diff --git a/crates/ppvm-lindblad/src/mf_expm.rs b/crates/ppvm-lindblad/src/mf_expm.rs new file mode 100644 index 000000000..f2e5b73f3 --- /dev/null +++ b/crates/ppvm-lindblad/src/mf_expm.rs @@ -0,0 +1,306 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Matrix-free `exp(dt · L*) · b` for the real (`f64`) path, driven by the +//! external `quspin-expm` crate. +//! +//! Instead of materialising the in-basis-restricted generator as a CSR, the +//! per-column generator action is computed ONCE per expm call (via +//! [`build_mf_cols`]) and reused, CSC-style, across every Krylov/Taylor matvec +//! by [`CscOp`] (a [`quspin_types::LinearOperator`]) fed to +//! [`quspin_expm::ExpmOp::from_parts`]. Each matvec is then a cheap CSC +//! scatter; the Pauli-commutator action is never recomputed per matvec. +//! `from_parts` (rather than `ExpmOp::new`) supplies the diagonal shift `μ`, +//! the partition count `s`, and the truncation order `m*` directly, bypassing +//! quspin's adaptive parameter selection — so the 1-norm *estimator* and +//! `dot_transpose` are never invoked on the single-vector `apply` path; only +//! [`LinearOperator::dot`] runs. +//! +//! `μ`, the trace, and the exact column 1-norm of `A − μ·I` are computed in +//! the same single action pass as the cache. The `(m, s)` Taylor partition is +//! picked with the tolerance-matched tables in [`crate::expm`]: a relaxed +//! `tol=1e-6` table when the PC prunes coarsely (`drop_tol ≥ 1e-4`), else the +//! double-precision table (keeping the exact-reference test paths bit-exact). + +use crate::{LindbladSpec, Word, build_basis_index, expm}; +use fxhash::{FxBuildHasher, FxHashMap}; +use num::Complex; +use quspin_types::{ExpmComputation, LinearOperator, QuSpinError}; +use rayon::prelude::*; + +/// Per-column in-basis action of the real generator `M`, plus the data the +/// `(m, s)`/`μ` selection needs — all from ONE action pass over the basis. +/// +/// Returns `(cols, per_col)` where `cols[c]` holds `(row, coeff)` for every +/// action output of `L*(basis[c])` that lands back in `basis` (CSC column +/// `c`), and `per_col[c] = (raw, diag)` with `raw = Σ|coeff|` over ALL action +/// outputs (in- and out-of-basis, an upper bound on the column 1-norm) and +/// `diag` the coefficient of the output Word equal to the input Word. The +/// cache is reused by [`CscOp`] across every Krylov/Taylor matvec. +/// CSC columns of the cached in-basis action: `cols[c]` = `(row, coeff)`. +type MfCols = Vec>; +/// Per-column `(raw, diag)` for the `μ`/1-norm selection. +type MfPerCol = Vec<(f64, f64)>; + +fn build_mf_cols( + spec: &LindbladSpec, + basis: &[Word], + index: &FxHashMap, +) -> (MfCols, MfPerCol) { + basis + .par_iter() + .map_init( + || { + ( + Vec::::with_capacity(spec.n_qubits()), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), p| { + let terms = spec.compute_action_terms(p, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + let mut raw = 0.0; + let mut diag = 0.0; + for (w, c) in terms.iter() { + raw += c.abs(); + if w == p { + diag = *c; + } + if let Some(&row) = index.get(w) { + out.push((row, *c)); + } + } + (out, (raw, diag)) + }, + ) + .unzip() +} + +/// Borrowed CSC-style view of an in-basis-restricted generator `M`, backed +/// by a cached per-column action computed once per expm call +/// ([`build_mf_cols`]). `dot` performs the CSC matvec `y = M·x` against the cache; the +/// remaining `LinearOperator` entry points are unused on the `from_parts` + +/// single-vector `apply` path. +/// +/// Borrowed, not owned: `quspin-types` provides a blanket `LinearOperator` +/// impl for `&T`, so `ExpmOp::from_parts(op, ...)` accepts a `CscOp` by +/// value while it keeps borrowing `cols`. +pub(crate) struct CscOp<'a, T> { + pub(crate) cols: &'a [Vec<(u32, T)>], + pub(crate) dim: usize, +} + +impl LinearOperator for CscOp<'_, T> +where + T: ExpmComputation + + Copy + + PartialEq + + num::Zero + + std::ops::AddAssign + + std::ops::Mul + + Send + + Sync, +{ + fn dim(&self) -> usize { + self.dim + } + + fn parallel_hint(&self) -> bool { + // `dot` parallelises internally over column chunks, and we drive the + // sequential single-vector `apply` path; never let quspin run its + // persistent-thread pool on top of our rayon parallelism. + false + } + + fn dot(&self, overwrite: bool, input: &[T], output: &mut [T]) -> Result<(), QuSpinError> { + let n = self.dim; + if n == 0 { + return Ok(()); + } + let num_threads = rayon::current_num_threads().max(1); + let chunk_size = n.div_ceil(num_threads); + + // Parallelise over column chunks; each thread accumulates into a dense + // local `y` of length `dim`, reading the cached action; the partials + // are reduced into `output` sequentially at the end. + let partial_ys: Vec> = self + .cols + .par_chunks(chunk_size) + .enumerate() + .map(|(chunk_idx, chunk)| { + let c_offset = chunk_idx * chunk_size; + let mut y_local = vec![T::zero(); n]; + for (c_local, col) in chunk.iter().enumerate() { + let xc = input[c_offset + c_local]; + if xc == T::zero() { + continue; + } + for &(row, val) in col.iter() { + y_local[row as usize] += val * xc; + } + } + y_local + }) + .collect(); + + if overwrite { + output.fill(T::zero()); + } + for partial in &partial_ys { + for (oi, &pi) in output.iter_mut().zip(partial.iter()) { + *oi += pi; + } + } + Ok(()) + } + + fn trace(&self) -> T { + // Computed eagerly by the callers; never reached on the + // `from_parts` + single-vector `apply` path. + unreachable!("CscOp::trace not used on the from_parts apply path") + } + + fn onenorm(&self, _shift: T) -> ::Real { + unreachable!("CscOp::onenorm not used on the from_parts apply path") + } + + fn dot_transpose( + &self, + _overwrite: bool, + _input: &[T], + _output: &mut [T], + ) -> Result<(), QuSpinError> { + Err(QuSpinError::RuntimeError( + "CscOp: dot_transpose not used on the from_parts apply path".into(), + )) + } + + fn dot_many( + &self, + _overwrite: bool, + _input: ndarray::ArrayView2<'_, T>, + _output: ndarray::ArrayViewMut2<'_, T>, + ) -> Result<(), QuSpinError> { + Err(QuSpinError::RuntimeError( + "CscOp: dot_many not used on the from_parts apply path".into(), + )) + } + + fn dot_chunk( + &self, + _overwrite: bool, + _input: &[T], + _output_chunk: &mut [T], + _row_start: usize, + ) -> Result<(), QuSpinError> { + Err(QuSpinError::RuntimeError( + "CscOp: dot_chunk not used on the from_parts apply path".into(), + )) + } + + fn dot_transpose_chunk( + &self, + _input: &[T], + _output: &[::Atomic], + _rows: std::ops::Range, + ) -> Result<(), QuSpinError> { + Err(QuSpinError::RuntimeError( + "CscOp: dot_transpose_chunk not used on the from_parts apply path".into(), + )) + } +} + +/// Compute `exp(dt · M) · coeffs` for the in-basis-restricted generator +/// `M`, matrix-free, via `quspin-expm`. Returns a fresh `Vec` of length +/// `basis.len()`. +/// +/// One matrix-free pass extracts the diagonal shift `μ = tr(M)/n` and the +/// exact column 1-norm of `M − μ·I`; from `‖dt·(M−μI)‖₁` we pick the Taylor +/// partition `(m*, s)` and hand everything to +/// [`quspin_expm::ExpmOp::from_parts`]. +pub(crate) fn expm_apply_mf( + spec: &LindbladSpec, + basis: &[Word], + dt: f64, + coeffs: &[f64], + drop_tol: f64, +) -> Vec { + let n = basis.len(); + if n == 0 { + return Vec::new(); + } + + // ONE action pass: build the CSC cache `cols` (reused across every matvec) + // and, in the same pass, `per_col = (raw, diag)` for the `μ`/1-norm + // selection. `raw = Σ|coeff|` (all outputs), `diag` = coeff of the + // output == input term. From these: `trace = Σ diag`, `μ = trace/n`, and + // the column 1-norm of `M − μ·I` is `raw − |diag| + |diag − μ|`. + let index = build_basis_index(basis); + let (cols, per_col) = build_mf_cols(spec, basis, &index); + + let trace: f64 = per_col.iter().map(|(_, d)| *d).sum(); + let mu = trace / n as f64; + let onenorm = per_col + .iter() + .map(|(raw, diag)| raw - diag.abs() + (diag - mu).abs()) + .fold(0.0_f64, f64::max); + + // Pick the Taylor backward-error tolerance to match the basis truncation: + // when the PC prunes coarsely (drop_tol >= 1e-4) a double-precision exp is + // ~10 orders more accurate than the state it acts on, so the relaxed + // (tol=1e-6, still >=100x tighter than the cut) table is used — it admits a + // lower-degree Taylor polynomial and cuts the SpMV count with no effect on + // the truncated result. At tight/zero drop_tol we keep double precision so + // the exact-reference paths (orbit-rep / merged) still agree bit-for-bit. + let t_norm = dt.abs() * onenorm; + let (m_star, s, expm_tol) = if drop_tol >= 1e-4 { + let (m, s) = expm::select_ms_loose(t_norm); + (m, s, 1e-6_f64) + } else { + let (m, s) = expm::select_ms(t_norm); + (m, s, 1e-12_f64) + }; + + let mut v = coeffs.to_vec(); + let op = CscOp { + cols: &cols, + dim: n, + }; + let expm = quspin_expm::ExpmOp::from_parts(op, dt, mu, s as usize, m_star as usize, expm_tol); + expm.apply(ndarray::ArrayViewMut1::from(v.as_mut_slice())) + .expect("expm apply"); + v +} + +/// `exp(dt · M) · b` where `M` is the REAL in-basis-restricted generator but +/// the input vector `b` is complex. Because `M` is real, +/// `exp(dt·M)·(re + i·im) = exp(dt·M)·re + i·exp(dt·M)·im`, so we split the +/// complex vector into its real and imaginary parts, run two real +/// matrix-free applies, and recombine. Used by the test-only full-space +/// complex reference step. +#[cfg(test)] +pub(crate) fn expm_apply_mf_cxvec( + spec: &LindbladSpec, + basis: &[Word], + dt: f64, + b: &[Complex], + drop_tol: f64, +) -> Vec> { + let n = basis.len(); + if n == 0 { + return Vec::new(); + } + let re: Vec = b.iter().map(|z| z.re).collect(); + let im: Vec = b.iter().map(|z| z.im).collect(); + let re_out = expm_apply_mf(spec, basis, dt, &re, drop_tol); + let im_out = expm_apply_mf(spec, basis, dt, &im, drop_tol); + re_out + .into_iter() + .zip(im_out) + .map(|(r, i)| Complex::new(r, i)) + .collect() +} diff --git a/crates/ppvm-lindblad/src/spec.rs b/crates/ppvm-lindblad/src/spec.rs new file mode 100644 index 000000000..4a8046eea --- /dev/null +++ b/crates/ppvm-lindblad/src/spec.rs @@ -0,0 +1,300 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Precompiled Lindbladian: construction and the single-Pauli `L*` kernel. + +use crate::Error; +use crate::algebra::{anti_commutes, comm_product, pauli_mul, phase_factor}; +use crate::word::{MAX_QUBITS, Word, parse_pauli_string, word_support}; +use fxhash::FxHashMap; +use num::Complex; + +/// Parsed Hamiltonian term. +#[derive(Clone)] +struct HTerm { + word: Word, + coeff: f64, +} + +/// One Pauli term in a complex linear combination (a single summand of +/// `L = Σ_a λ_a P_a` or of the precomputed `L†L`). +#[derive(Clone)] +struct PauliTerm { + word: Word, + coeff: Complex, +} + +/// One jump operator `L_k` with rate `γ_k`. The `HermitianPauli` variant +/// is a fast path; `General` handles arbitrary complex Pauli sums. +#[derive(Clone)] +enum JumpKind { + HermitianPauli { + word: Word, + rate: f64, + }, + General { + terms: Vec, // L = Σ_a λ_a P_a + dagger_dagger: Vec, // L†L = Σ_c μ_c P_c (μ_c ∈ ℝ) + rate: f64, + }, +} + +/// Expand `L†L = (Σ_a λ_a P_a)† (Σ_b λ_b P_b) = Σ_{a,b} λ_a* λ_b P_a P_b` +/// as a Pauli linear combination, dropping FP-noise zeros. Coefficients are +/// real because `L†L` is Hermitian; we keep them complex for arithmetic +/// uniformity. +fn precompute_ldagger_l(terms: &[PauliTerm]) -> Vec { + let zero = Complex::new(0.0, 0.0); + let mut acc: FxHashMap> = FxHashMap::default(); + for a in terms { + for b in terms { + let (word, phase) = pauli_mul(&a.word, &b.word); + let coeff = a.coeff.conj() * b.coeff * phase_factor(phase); + *acc.entry(word).or_insert(zero) += coeff; + } + } + acc.into_iter() + .filter(|(_, c)| c.norm() > 1e-14) + .map(|(word, coeff)| PauliTerm { word, coeff }) + .collect() +} + +/// Union of `index[q]` for each `q ∈ p_support`, deduped. +#[inline] +fn candidate_terms(p_support: &[u32], index: &[Vec], scratch: &mut Vec) { + scratch.clear(); + for &q in p_support { + scratch.extend_from_slice(&index[q as usize]); + } + scratch.sort_unstable(); + scratch.dedup(); +} + +/// Precompiled Lindbladian. Constructed once from string-form Hamiltonian +/// terms + jump operators; reused across many calls to [`Self::action`], +/// [`Self::leakage`], [`Self::generator`]. `L*(p)` is recomputed on every +/// call rather than cached: for sparse-local Hamiltonians a per-word cache +/// costs more than the recompute (hash lookup ≳ recompute) and its several +/// KB per cached word dominate memory at large basis sizes. +pub struct LindbladSpec { + n_qubits: usize, + h_terms: Vec, + j_kinds: Vec, + /// `h_support[q]` = indices of Hamiltonian terms acting on qubit `q`. + h_support: Vec>, + /// `j_support[q]` = indices of jumps whose support contains qubit `q`. + j_support: Vec>, +} + +/// User-facing description of one jump operator: a complex Pauli linear +/// combination together with its rate. +#[derive(Clone, Debug)] +pub struct JumpInput { + /// `(pauli_string, λ)` pairs forming `L_k = Σ_a λ_a P_a`. + pub lincomb: Vec<(String, Complex)>, + /// Non-negative GKSL rate `γ_k`. + pub rate: f64, +} + +impl LindbladSpec { + /// Construct a Lindbladian spec from Hamiltonian terms and jump operators. + /// + /// `h_terms` are `(pauli_string, coefficient)` pairs forming the Hermitian + /// Hamiltonian. Each jump operator is a complex Pauli linear combination; + /// a length-1 jump with imaginary part `0` is routed to the Hermitian-Pauli + /// fast path (with rate scaled by the squared real coefficient). + pub fn new( + n_qubits: usize, + h_terms: &[(String, f64)], + jumps: &[JumpInput], + ) -> Result { + if n_qubits > MAX_QUBITS { + return Err(Error::TooManyQubits { got: n_qubits }); + } + + let mut h_parsed: Vec = Vec::with_capacity(h_terms.len()); + let mut h_support_idx: Vec> = vec![Vec::new(); n_qubits]; + for (i, (s, c)) in h_terms.iter().enumerate() { + let (word, support) = parse_pauli_string(s, n_qubits)?; + for q in support { + h_support_idx[q as usize].push(i as u32); + } + h_parsed.push(HTerm { word, coeff: *c }); + } + + let mut j_kinds: Vec = Vec::with_capacity(jumps.len()); + let mut j_support_idx: Vec> = vec![Vec::new(); n_qubits]; + for (k, jump) in jumps.iter().enumerate() { + if jump.rate < 0.0 { + return Err(Error::NegativeRate { + index: k, + rate: jump.rate, + }); + } + if jump.lincomb.is_empty() { + return Err(Error::EmptyLincomb { index: k }); + } + + // Fast path: single-term, purely real → Hermitian Pauli. + if jump.lincomb.len() == 1 && jump.lincomb[0].1.im == 0.0 { + let (s, c) = &jump.lincomb[0]; + let (word, support) = parse_pauli_string(s, n_qubits)?; + for q in support { + j_support_idx[q as usize].push(k as u32); + } + j_kinds.push(JumpKind::HermitianPauli { + word, + rate: jump.rate * c.re * c.re, + }); + continue; + } + + // General path: parse all terms, precompute L†L, record union support. + let mut terms: Vec = Vec::with_capacity(jump.lincomb.len()); + let mut union_support: std::collections::BTreeSet = + std::collections::BTreeSet::new(); + for (s, c) in &jump.lincomb { + let (word, support) = parse_pauli_string(s, n_qubits)?; + for q in &support { + union_support.insert(*q); + } + terms.push(PauliTerm { word, coeff: *c }); + } + for q in union_support { + j_support_idx[q as usize].push(k as u32); + } + let dagger_dagger = precompute_ldagger_l(&terms); + j_kinds.push(JumpKind::General { + terms, + dagger_dagger, + rate: jump.rate, + }); + } + + Ok(Self { + n_qubits, + h_terms: h_parsed, + j_kinds, + h_support: h_support_idx, + j_support: j_support_idx, + }) + } + + pub fn n_qubits(&self) -> usize { + self.n_qubits + } + + pub fn num_h_terms(&self) -> usize { + self.h_terms.len() + } + + pub fn num_jump_terms(&self) -> usize { + self.j_kinds.len() + } + + /// Apply `L*` to a single Pauli string `p`. Returns the output Pauli + /// strings and their real coefficients (zero entries omitted). + pub fn action(&self, p: &Word) -> Vec<(Word, f64)> { + let mut out: FxHashMap = FxHashMap::default(); + let mut s1 = Vec::new(); + let mut s2 = Vec::new(); + self.accumulate_action(p, 1.0, &mut out, &mut s1, &mut s2); + out.into_iter().filter(|(_, c)| *c != 0.0).collect() + } + + /// Compute the unscaled list of `(output, coefficient)` pairs that + /// `L*(p)` contributes (without the input coefficient). + pub(crate) fn compute_action_terms( + &self, + p: &Word, + scratch_support: &mut Vec, + scratch_cands: &mut Vec, + scratch_local: &mut FxHashMap>, + ) -> Vec<(Word, f64)> { + word_support(p, scratch_support); + let zero = Complex::new(0.0, 0.0); + scratch_local.clear(); + let local = scratch_local; + + // ── i [H, p] ───────────────────────────────────────────────── + candidate_terms(scratch_support, &self.h_support, scratch_cands); + for &i in scratch_cands.iter() { + let h = &self.h_terms[i as usize]; + let (r, eps) = comm_product(&h.word, p); + if eps != 0.0 { + *local.entry(r).or_insert(zero) += Complex::new(h.coeff * eps, 0.0); + } + } + + // ── dissipator ─────────────────────────────────────────────── + candidate_terms(scratch_support, &self.j_support, scratch_cands); + for &k in scratch_cands.iter() { + match &self.j_kinds[k as usize] { + JumpKind::HermitianPauli { word, rate } => { + if anti_commutes(word, p) { + *local.entry(*p).or_insert(zero) += Complex::new(-2.0 * *rate, 0.0); + } + } + JumpKind::General { + terms, + dagger_dagger, + rate, + } => { + let rate_c = Complex::new(*rate, 0.0); + // Sandwich: γ Σ_{a,b} λ_a* λ_b P_a p P_b. + for a in terms { + let (r_ap, phi1) = pauli_mul(&a.word, p); + for b in terms { + let (s, phi2) = pauli_mul(&r_ap, &b.word); + let coeff = + a.coeff.conj() * b.coeff * phase_factor(phi1 + phi2) * rate_c; + *local.entry(s).or_insert(zero) += coeff; + } + } + // -1/2 γ {L†L, p}. For Hermitian Pauli P_c and Pauli p, + // {P_c, p} = 2·sign·R if they commute (P_c·p = sign·R), + // = 0 if they anti-commute. + for c_term in dagger_dagger { + let (r, phase) = pauli_mul(&c_term.word, p); + if phase & 1 == 0 { + let sign = if phase == 0 { 1.0 } else { -1.0 }; + let coeff = -c_term.coeff * rate_c * Complex::new(sign, 0.0); + *local.entry(r).or_insert(zero) += coeff; + } + } + } + } + } + + // L* preserves Hermiticity; imaginary parts must cancel to FP noise. + // `drain()` empties `scratch_local` so its allocation can be reused + // by the next call on the same thread (`Vec` keeps capacity). + local + .drain() + .filter_map(|(w, c)| { + debug_assert!( + c.im.abs() < 1e-9, + "L*(p) produced non-real coefficient {c}; bug in dissipator" + ); + if c.re == 0.0 { None } else { Some((w, c.re)) } + }) + .collect() + } + + /// Accumulate `scale · L*(p)` into `out`. + fn accumulate_action( + &self, + p: &Word, + scale: f64, + out: &mut FxHashMap, + scratch_support: &mut Vec, + scratch_cands: &mut Vec, + ) { + let mut scratch_local = FxHashMap::default(); + let terms = + self.compute_action_terms(p, scratch_support, scratch_cands, &mut scratch_local); + for (w, c) in terms.iter() { + *out.entry(*w).or_insert(0.0) += scale * c; + } + } +} diff --git a/crates/ppvm-lindblad/src/step.rs b/crates/ppvm-lindblad/src/step.rs new file mode 100644 index 000000000..fe77ce83f --- /dev/null +++ b/crates/ppvm-lindblad/src/step.rs @@ -0,0 +1,272 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Predictor-corrector adaptive step `O ← exp(dt·L*) O`. + +use crate::spec::LindbladSpec; +use crate::word::Word; +use crate::{Error, PcStepConfig, mf_expm}; +use fxhash::FxHashSet; +use std::time::Instant; + +/// Per-phase timing breakdown (microseconds) returned by +/// [`LindbladSpec::pc_step_timed`]. +#[derive(Default, Clone, Copy, Debug)] +pub struct PcStepTimings { + pub leakage1_us: u64, + pub expand1_us: u64, + pub expm1_us: u64, + pub leakage2_us: u64, + pub expand2_us: u64, + pub expm2_us: u64, +} + +impl PcStepTimings { + pub fn total_us(&self) -> u64 { + self.leakage1_us + + self.expand1_us + + self.expm1_us + + self.leakage2_us + + self.expand2_us + + self.expm2_us + } +} + +/// Clock for one `pc_step` phase. Disarmed (`None`) on the untimed path, so +/// [`LindbladSpec::pc_step`] pays no `Instant` syscalls. +struct Phase(Option); + +impl Phase { + fn start(timed: bool) -> Self { + Self(timed.then(Instant::now)) + } + + fn stop(self, slot: &mut u64) { + if let Some(t0) = self.0 { + *slot = t0.elapsed().as_micros() as u64; + } + } +} + +/// Compact `basis` / `coeffs` in place: drop entries whose absolute +/// coefficient is below `drop_tol` unless the word appears in `protected`. +/// No-op when `drop_tol ≤ 0`. +fn prune_basis(basis: &mut Vec, coeffs: &mut Vec, drop_tol: f64, protected: &[Word]) { + if drop_tol <= 0.0 { + return; + } + debug_assert_eq!(basis.len(), coeffs.len()); + let protected_set: FxHashSet<&Word> = protected.iter().collect(); + let mut write = 0; + for read in 0..basis.len() { + if coeffs[read].abs() >= drop_tol || protected_set.contains(&basis[read]) { + if write != read { + basis.swap(write, read); + coeffs.swap(write, read); + } + write += 1; + } + } + basis.truncate(write); + coeffs.truncate(write); +} + +/// Global max-basis cap (PauliStrings.jl-style top-M trim): keep only the +/// `max_basis` largest-|coeff| terms (protected strings always kept), +/// dropping the rest. Rank-based total-basis bound; dual of `drop_tol`. +/// A `max_basis` large enough to cover the whole basis is a no-op. +fn cap_basis(basis: &mut Vec, coeffs: &mut Vec, max_basis: usize, protected: &[Word]) { + if basis.len() <= max_basis { + return; + } + let protected_set: FxHashSet<&Word> = protected.iter().collect(); + let n_prot = basis.iter().filter(|w| protected_set.contains(w)).count(); + let slots = max_basis.saturating_sub(n_prot); + let mut mags: Vec = basis + .iter() + .zip(coeffs.iter()) + .filter(|(w, _)| !protected_set.contains(w)) + .map(|(_, c)| c.abs()) + .collect(); + let cutoff = if slots == 0 { + f64::INFINITY + } else if slots >= mags.len() { + return; + } else { + let k = slots - 1; + mags.select_nth_unstable_by(k, |a, b| { + b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal) + }); + mags[k] + }; + let mut write = 0; + for read in 0..basis.len() { + if protected_set.contains(&basis[read]) || coeffs[read].abs() >= cutoff { + if write != read { + basis.swap(write, read); + coeffs.swap(write, read); + } + write += 1; + } + } + basis.truncate(write); + coeffs.truncate(write); +} + +/// Add the largest leakage strings to the basis, up to the available room +/// `room = max_basis − basis.len()` — so the in-step basis (hence the +/// expm/leakage peak memory) never exceeds `max_basis`. New strings get +/// coefficient 0; the surrounding expm fills them. No magnitude filter: the +/// top-`room` by `|leakage|` are added (a large `max_basis` adds them all). +fn add_leakage_capped( + basis: &mut Vec, + coeffs: &mut Vec, + mut leak: Vec<(Word, f64)>, + max_basis: usize, +) { + let room = max_basis.saturating_sub(basis.len()); + if leak.len() > room { + if room > 0 { + leak.select_nth_unstable_by(room - 1, |a, b| { + b.1.abs() + .partial_cmp(&a.1.abs()) + .unwrap_or(std::cmp::Ordering::Equal) + }); + } + leak.truncate(room); + } + for (w, _) in leak { + basis.push(w); + coeffs.push(0.0); + } +} + +impl LindbladSpec { + /// One predictor-corrector step `O ← exp(dt·L*) O` in the adaptive + /// real-coefficient Pauli basis: first-hop leakage admission, predictor + /// exponential, second-hop admission from the predicted state, corrector + /// exponential from the saved pre-step state, then truncation (prune + + /// rank cap) per [`PcStepConfig`]. Exact in `dt` within the working + /// basis — the only error is basis truncation. + /// + /// `protected` words are never dropped. All tuning knobs live in `cfg`. + pub fn pc_step( + &self, + basis: &mut Vec, + coeffs: &mut Vec, + dt: f64, + protected: &[Word], + cfg: &PcStepConfig, + ) -> Result<(), Error> { + self.run_in_pool(cfg, |this| { + this.pc_step_inner(basis, coeffs, dt, protected, cfg, false) + .map(|_| ()) + }) + } + + /// Same as [`Self::pc_step`] but also returns a per-phase timing + /// breakdown (microseconds), for profiling parallel scaling and hot + /// spots. + pub fn pc_step_timed( + &self, + basis: &mut Vec, + coeffs: &mut Vec, + dt: f64, + protected: &[Word], + cfg: &PcStepConfig, + ) -> Result { + self.run_in_pool(cfg, |this| { + this.pc_step_inner(basis, coeffs, dt, protected, cfg, true) + }) + } + + fn run_in_pool( + &self, + cfg: &PcStepConfig, + f: impl FnOnce(&Self) -> Result + Send, + ) -> Result { + if let Some(n) = cfg.num_threads { + let pool = rayon::ThreadPoolBuilder::new() + .num_threads(n) + .build() + .map_err(|e| Error::Internal(format!("rayon pool build: {e}")))?; + pool.install(|| f(self)) + } else { + f(self) + } + } + + fn pc_step_inner( + &self, + basis: &mut Vec, + coeffs: &mut Vec, + dt: f64, + protected: &[Word], + cfg: &PcStepConfig, + timed: bool, + ) -> Result { + let PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + .. + } = *cfg; + // Admission bound: enrichment may grow the live basis to `admit` + // >= `max_basis`; the final `cap_basis` then keeps the top- + // `max_basis` strings by evolved |coeff| over the whole union + // (retained + admitted) — rank displacement. With `admit_basis = + // None` admission is bounded by `max_basis` itself, `cap_basis` is + // a no-op, and membership turnover requires `drop_tol > 0`. + let admit = admit_basis.unwrap_or(max_basis).max(max_basis); + let tau_add = tau_add.unwrap_or(0.0); + let mut t = PcStepTimings::default(); + + // 1. First-hop expansion. After this, `coeffs` contains the pre-step + // coefficients followed by zeros for the newly-added leakage strings. + // We rely on `coeffs` itself as the pre-step buffer for the corrector + // — no `.clone()` is needed because `expm_step` only borrows it. + let p = Phase::start(timed); + let leak = self.leakage_with_prune(basis, coeffs, protected, admit, tau_add)?; + p.stop(&mut t.leakage1_us); + + let p = Phase::start(timed); + add_leakage_capped(basis, coeffs, leak, admit); + p.stop(&mut t.expand1_us); + + // 2. Predictor: `expm_step` reads `coeffs` immutably and returns a + // new owned vector with the predicted state. + let p = Phase::start(timed); + let coeffs_predict = self.expm_step(basis, dt, coeffs, drop_tol); + p.stop(&mut t.expm1_us); + + // 3. Second-hop expansion from the predicted state. After leakage2 + // we no longer need `coeffs_predict`. Extend `coeffs` with zeros for + // any newly-added second-hop strings so it remains a valid input + // (pre-step state) for the corrector. + let p = Phase::start(timed); + let leak2 = self.leakage_with_prune(basis, &coeffs_predict, protected, admit, tau_add)?; + p.stop(&mut t.leakage2_us); + drop(coeffs_predict); + + let p = Phase::start(timed); + add_leakage_capped(basis, coeffs, leak2, admit); + p.stop(&mut t.expand2_us); + + // 4. Corrector: redo from pre-step state on the doubly-enlarged basis. + let p = Phase::start(timed); + *coeffs = self.expm_step(basis, dt, coeffs, drop_tol); + p.stop(&mut t.expm2_us); + + // 5. Prune basis entries below `drop_tol` (protected words never dropped). + prune_basis(basis, coeffs, drop_tol, protected); + cap_basis(basis, coeffs, max_basis, protected); + Ok(t) + } + + /// Compute `exp(dt · M) · b` for the in-basis-restricted generator + /// `M`, matrix-free, via `quspin-expm` (see [`crate::mf_expm`]). + fn expm_step(&self, basis: &[Word], dt: f64, b: &[f64], drop_tol: f64) -> Vec { + mf_expm::expm_apply_mf(self, basis, dt, b, drop_tol) + } +} diff --git a/crates/ppvm-lindblad/src/tests.rs b/crates/ppvm-lindblad/src/tests.rs new file mode 100644 index 000000000..4217b3e13 --- /dev/null +++ b/crates/ppvm-lindblad/src/tests.rs @@ -0,0 +1,275 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +use super::*; +use fxhash::FxHashMap; +use num::Complex; + +/// Test-only full-space complex predictor-corrector step, UNTRUNCATED: +/// adds every nonzero leakage string (two hops) and applies the exact +/// in-basis exponential to the complex coefficient vector. Reference +/// bridge between the real `pc_step` and the orbit-rep path. +fn pc_step_complex_full( + spec: &LindbladSpec, + basis: &mut Vec, + coeffs: &mut Vec>, + dt: f64, +) { + let protected: Vec = Vec::new(); + let leak = spec.leakage_complex(basis, coeffs, &protected).unwrap(); + for (w, v) in leak { + if v.norm() > 0.0 { + basis.push(w); + coeffs.push(Complex::new(0.0, 0.0)); + } + } + let predict = mf_expm::expm_apply_mf_cxvec(spec, basis, dt, coeffs, 0.0); + let leak2 = spec.leakage_complex(basis, &predict, &protected).unwrap(); + drop(predict); + for (w, v) in leak2 { + if v.norm() > 0.0 { + basis.push(w); + coeffs.push(Complex::new(0.0, 0.0)); + } + } + *coeffs = mf_expm::expm_apply_mf_cxvec(spec, basis, dt, coeffs, 0.0); +} + +fn jump_hpauli(s: &str, rate: f64) -> JumpInput { + JumpInput { + lincomb: vec![(s.to_string(), Complex::new(1.0, 0.0))], + rate, + } +} + +#[test] +fn z_dephasing_action_on_x() { + // L = Z on a single qubit; L*(X) = γ(ZXZ - X) = γ(-X - X) = -2γ X. + let spec = LindbladSpec::new( + 1, + &[("X".to_string(), 0.0)], // no Hamiltonian + &[jump_hpauli("Z", 0.5)], + ) + .unwrap(); + let (x, _) = parse_pauli_string("X", 1).unwrap(); + let terms = spec.action(&x); + assert_eq!(terms.len(), 1); + assert!((terms[0].1 - (-1.0)).abs() < 1e-12); // -2·0.5 = -1 +} + +#[test] +fn amplitude_damping_action_on_z() { + // Single-qubit σ⁻ jump: L*(Z) = -γ(I + Z). With γ=1 we expect + // I coefficient = -1, Z coefficient = -1. + let sigma_minus = JumpInput { + lincomb: vec![ + ("X".to_string(), Complex::new(0.5, 0.0)), + ("Y".to_string(), Complex::new(0.0, -0.5)), + ], + rate: 1.0, + }; + let spec = LindbladSpec::new(1, &[], &[sigma_minus]).unwrap(); + let (z, _) = parse_pauli_string("Z", 1).unwrap(); + let terms = spec.action(&z); + let (i_word, _) = parse_pauli_string("I", 1).unwrap(); + let mut i_coeff = 0.0; + let mut z_coeff = 0.0; + for (w, c) in &terms { + if w == &i_word { + i_coeff = *c; + } else if w == &z { + z_coeff = *c; + } + } + assert!((i_coeff - (-1.0)).abs() < 1e-10, "I coeff = {i_coeff}"); + assert!((z_coeff - (-1.0)).abs() < 1e-10, "Z coeff = {z_coeff}"); +} + +#[test] +fn word_codec_roundtrip() { + let codes = [0u8, 1, 2, 3, 1, 0, 3, 2]; + let w = word_from_codes(&codes).unwrap(); + let mut out = vec![0u8; codes.len()]; + codes_from_word(&w, &mut out); + assert_eq!(out.as_slice(), &codes); +} + +/// The full-space complex step at momentum k=0 must reproduce the real +/// pc_step on the same trajectory exactly. +#[test] +fn complex_full_matches_real_at_kzero() { + let n = 4usize; + let dt = 0.01f64; + let n_steps = 5usize; + let mut h_terms: Vec<(String, f64)> = Vec::new(); + for j in 0..n { + let nxt = (j + 1) % n; + for op in ["X", "Y"] { + let mut s = vec!['I'; n]; + s[j] = op.chars().next().unwrap(); + s[nxt] = op.chars().next().unwrap(); + h_terms.push((s.into_iter().collect(), 1.0)); + } + } + let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap(); + + let mut basis_r: Vec = (0..n) + .map(|j| { + let mut s = vec!['I'; n]; + s[j] = 'Z'; + let st: String = s.into_iter().collect(); + let (w, _) = parse_pauli_string(&st, n).unwrap(); + w + }) + .collect(); + let mut coeffs_r: Vec = vec![1.0; n]; + + let mut basis_c = basis_r.clone(); + let mut coeffs_c: Vec> = coeffs_r.iter().map(|&v| Complex::new(v, 0.0)).collect(); + + let protected: Vec = Vec::new(); + for _ in 0..n_steps { + // Large max_basis: rank cap never binds, so the real path + // enriches fully (adds every leakage string). Match the + // complex path by setting its tau_add=0.0 (also full + // enrichment) so the two stay in lock-step at k=0. + spec.pc_step( + &mut basis_r, + &mut coeffs_r, + dt, + &protected, + &PcStepConfig { + max_basis: 10_000_000, + ..Default::default() + }, + ) + .unwrap(); + pc_step_complex_full(&spec, &mut basis_c, &mut coeffs_c, dt); + } + // Match as (word → coeff) maps. + let map_r: FxHashMap = basis_r.into_iter().zip(coeffs_r).collect(); + let map_c: FxHashMap> = basis_c.into_iter().zip(coeffs_c).collect(); + assert_eq!( + map_r.len(), + map_c.len(), + "real and complex pc_step produced different basis sizes ({} vs {})", + map_r.len(), + map_c.len() + ); + let mut max_diff = 0.0_f64; + for (w, cr) in &map_r { + let cc = map_c + .get(w) + .copied() + .unwrap_or_else(|| panic!("word {:?} in real but not complex", w)); + assert!(cc.im.abs() < 1e-10, "expected zero imag at k=0, got {cc:?}"); + max_diff = max_diff.max((cr - cc.re).abs()); + } + assert!( + max_diff < 1e-10, + "real vs complex pc_step diverged: max |Δc| = {max_diff:e}" + ); +} + +/// Small-system end-to-end check that orbit-rep merging gives the +/// same physics as standard evolution, when no truncation is applied. +/// +/// Setup: n=4 qubit chain, PBC, translation-invariant XY Hamiltonian +/// `H = Σ_j (X_j X_{j+1} + Y_j Y_{j+1})`, no dissipation. Initial +/// operator `O(0) = Σ_j Z_j` is translation-invariant (k=0 sector). +/// +/// Run 10 pc_step iterations with `drop_tol = 0` (no truncation): +/// once without merging, once applying `canonicalize_pauli_sum` +/// after each step. Canonicalize the un-merged final state once at +/// the end. The two orbit-rep representations should be +/// bit-identical up to FP noise. +#[test] +fn pc_step_matches_symmetry_merged_on_small_chain() { + use ppvm_pauli_sum::symmetry::{TranslationGroup, canonicalize_pauli_sum}; + + let n = 4usize; + let dt = 0.05f64; + let n_steps = 10usize; + + // Build XY-chain Hamiltonian with PBC. 8 terms (4 bonds × {XX, YY}). + let mut h_terms: Vec<(String, f64)> = Vec::new(); + for j in 0..n { + let nxt = (j + 1) % n; + for op in ["X", "Y"] { + let mut s = vec!['I'; n]; + s[j] = op.chars().next().unwrap(); + s[nxt] = op.chars().next().unwrap(); + h_terms.push((s.into_iter().collect(), 1.0)); + } + } + // No dissipation. + let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap(); + let group = TranslationGroup::chain_1d(n); + + // Initial: O(0) = Σ_j Z_j (translation-invariant). + let mut basis_u: Vec = (0..n) + .map(|j| { + let mut s = vec!['I'; n]; + s[j] = 'Z'; + let st: String = s.into_iter().collect(); + let (w, _) = parse_pauli_string(&st, n).unwrap(); + w + }) + .collect(); + let mut coeffs_u: Vec = vec![1.0; n]; + + // Mirror state for the "with merging" run. + let mut basis_m = basis_u.clone(); + let mut coeffs_m = coeffs_u.clone(); + + let protected: Vec = Vec::new(); + for _ in 0..n_steps { + // max_basis == current basis size → room = 0: no leakage + // enrichment, only the expm step (the regime where merging + // commutes with evolution). drop_tol = 0 → no truncation. + let cfg_u = PcStepConfig { + max_basis: basis_u.len(), + ..Default::default() + }; + spec.pc_step(&mut basis_u, &mut coeffs_u, dt, &protected, &cfg_u) + .unwrap(); + + let cfg_m = PcStepConfig { + max_basis: basis_m.len(), + ..Default::default() + }; + spec.pc_step(&mut basis_m, &mut coeffs_m, dt, &protected, &cfg_m) + .unwrap(); + // Apply symmetry merging on the "with merging" run only. + canonicalize_pauli_sum(&mut basis_m, &mut coeffs_m, &group); + } + + // Canonicalize the un-merged final state once. + canonicalize_pauli_sum(&mut basis_u, &mut coeffs_u, &group); + + // Both representations should now be in orbit-rep form; compare + // as (word → coeff) maps with FP tolerance. + let map_u: FxHashMap = basis_u.into_iter().zip(coeffs_u).collect(); + let map_m: FxHashMap = basis_m.into_iter().zip(coeffs_m).collect(); + assert_eq!( + map_u.len(), + map_m.len(), + "merged basis size {} != post-merged-unmerged basis size {}", + map_m.len(), + map_u.len() + ); + let mut max_diff = 0.0f64; + for (w, c_u) in &map_u { + let c_m = map_m.get(w).copied().unwrap_or_else(|| { + panic!( + "rep {:?} present in un-merged-then-canonicalized but not in merged", + w + ); + }); + max_diff = max_diff.max((c_u - c_m).abs()); + } + assert!( + max_diff < 1e-9, + "with-merging vs without-merging diverged: max |Δc| = {max_diff:e}" + ); +} diff --git a/crates/ppvm-lindblad/src/word.rs b/crates/ppvm-lindblad/src/word.rs new file mode 100644 index 000000000..15c762d9b --- /dev/null +++ b/crates/ppvm-lindblad/src/word.rs @@ -0,0 +1,128 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Packed Pauli-word type and the string / `u8`-label codecs. + +use crate::Error; +use fxhash::FxBuildHasher; +use ppvm_pauli_word::word::PauliWord; +use ppvm_traits::PauliWordTrait; + +/// Pauli-word storage chunk: `u64` on 64-bit targets, `u32` elsewhere +/// (`bitvec` implements `BitStore` for `u64` only on 64-bit targets, so +/// e.g. wasm32 builds use four 32-bit chunks instead of two 64-bit ones). +#[cfg(target_pointer_width = "64")] +pub(crate) type Chunk = u64; +#[cfg(not(target_pointer_width = "64"))] +pub(crate) type Chunk = u32; + +/// Chunks per word; words pack up to 128 qubits on every target. +#[cfg(target_pointer_width = "64")] +pub(crate) const W_CHUNKS: usize = 2; +#[cfg(not(target_pointer_width = "64"))] +pub(crate) const W_CHUNKS: usize = 4; + +/// Maximum number of qubits supported by [`Word`]. +pub const MAX_QUBITS: usize = 128; + +/// The Pauli-word storage type used throughout this crate. +/// +/// `[Chunk; W_CHUNKS]` covers up to 128 qubits; the `FxBuildHasher` +/// matches the hash used by the `FxHashMap` keys we wrap with; +/// `REHASH=true` means `set()` keeps the cached hash in sync. +pub type Word = PauliWord<[Chunk; W_CHUNKS], FxBuildHasher, true>; + +/// Build a [`Word`] from a length-`n_qubits` slice of Pauli labels +/// (`0=I, 1=X, 2=Z, 3=Y` — the [`ppvm_traits::char::Pauli`] discriminants). +/// Sets all bits and rehashes once. +pub fn word_from_codes(codes: &[u8]) -> Result { + let n_qubits = codes.len(); + if n_qubits > MAX_QUBITS { + return Err(Error::TooManyQubits { got: n_qubits }); + } + let mut w = Word::new(n_qubits); + for (q, &b) in codes.iter().enumerate() { + if b > 3 { + return Err(Error::InvalidPauliCode { code: b }); + } + if b & 1 != 0 { + w.xbits.set(q, true); + } + if b & 2 != 0 { + w.zbits.set(q, true); + } + } + w.rehash(); + Ok(w) +} + +/// Inverse of [`word_from_codes`]: write `n_qubits` Pauli labels into `out`. +pub fn codes_from_word(w: &Word, out: &mut [u8]) { + debug_assert_eq!(out.len(), w.n_qubits()); + for (q, slot) in out.iter_mut().enumerate() { + let xb = w.xbits[q] as u8; + let zb = w.zbits[q] as u8; + *slot = xb | (zb << 1); + } +} + +/// Parse a `"IXYZ..."` string into a [`Word`] together with the list of +/// qubits where the Pauli is non-identity (the term's support). +pub fn parse_pauli_string(s: &str, n_qubits: usize) -> Result<(Word, Vec), Error> { + if n_qubits > MAX_QUBITS { + return Err(Error::TooManyQubits { got: n_qubits }); + } + let chars: Vec = s.chars().filter(|c| *c != '_').collect(); + if chars.len() != n_qubits { + return Err(Error::WrongLength { + expected: n_qubits, + got: chars.len(), + }); + } + let mut w = Word::new(n_qubits); + let mut support = Vec::new(); + for (q, c) in chars.into_iter().enumerate() { + match c { + 'I' => {} + 'X' => { + w.xbits.set(q, true); + support.push(q as u32); + } + 'Z' => { + w.zbits.set(q, true); + support.push(q as u32); + } + 'Y' => { + w.xbits.set(q, true); + w.zbits.set(q, true); + support.push(q as u32); + } + other => return Err(Error::InvalidPauliChar { c: other }), + } + } + w.rehash(); + Ok((w, support)) +} + +/// Compute the support (non-identity qubits) of `w`. +pub(crate) fn word_support(w: &Word, out: &mut Vec) { + out.clear(); + for q in 0..w.n_qubits() { + if w.xbits[q] || w.zbits[q] { + out.push(q as u32); + } + } +} + +/// Compact 64-bit hash of a [`Word`], used as the key in cache-friendly +/// membership tables: an `FxHashMap` over the basis has a working +/// set ~6× smaller than `FxHashMap`. The hash mixes the word's +/// cached hash once through `FxHasher` and never touches the 32-byte +/// payload. +#[inline(always)] +pub(crate) fn word_hash(w: &Word) -> u64 { + use std::hash::{Hash, Hasher}; + let mut h = fxhash::FxHasher::default(); + w.hash(&mut h); + h.finish() +} diff --git a/crates/ppvm-python-native/Cargo.toml b/crates/ppvm-python-native/Cargo.toml index 7d948575b..e926b570e 100644 --- a/crates/ppvm-python-native/Cargo.toml +++ b/crates/ppvm-python-native/Cargo.toml @@ -13,7 +13,14 @@ test = false [dependencies] bnum = "0.13.0" +# mimalloc reduces peak RSS for the allocation-heavy adaptive-Pauli (lindblad) +# paths; installed as the global allocator in lib.rs. +mimalloc = { version = "0.1", default-features = false } +num = "0.4.3" +# numpy 0.29 pairs with pyo3 0.29 (used by the lindblad + symmetry array bindings). +numpy = "0.29" paste = "1.0.15" +ppvm-lindblad = { version = "0.1.0", path = "../ppvm-lindblad" } ppvm-pauli-sum = { version = "0.1.0", path = "../ppvm-pauli-sum" } ppvm-stim = { version = "0.1.0", path = "../ppvm-stim", features = ["rayon"] } ppvm-tableau = { version = "0.1.0", path = "../ppvm-tableau" } diff --git a/crates/ppvm-python-native/src/lib.rs b/crates/ppvm-python-native/src/lib.rs index b323c1adc..13f4c22ec 100644 --- a/crates/ppvm-python-native/src/lib.rs +++ b/crates/ppvm-python-native/src/lib.rs @@ -1,12 +1,20 @@ // SPDX-FileCopyrightText: 2026 The PPVM Authors // SPDX-License-Identifier: Apache-2.0 +// mimalloc returns freed pages to the kernel more aggressively than the +// default system allocator. This materially reduces peak RSS for the +// allocation-heavy adaptive-Pauli paths (leakage + generator each +// allocate hundreds of MB of transient Vec data per pc_step). +#[global_allocator] +static GLOBAL: mimalloc::MiMalloc = mimalloc::MiMalloc; + use pyo3::exceptions::PyValueError; use pyo3::prelude::*; pub mod interface; pub mod interface_tableau; pub mod interface_tableau_sum; +pub mod lindblad; pub mod stim_program; pub(crate) fn flat_pairs(targets: &[usize]) -> PyResult> { @@ -298,4 +306,8 @@ pub mod _core { // Stim #[pymodule_export] pub use crate::stim_program::PyStimProgram; + + // Lindbladian time evolution + #[pymodule_export] + pub use crate::lindblad::LindbladSpec; } diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs new file mode 100644 index 000000000..25d4e0b0d --- /dev/null +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -0,0 +1,370 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! PyO3 wrapper around [`ppvm_lindblad::LindbladSpec`]. +//! +//! All algorithmic work — Pauli arithmetic, active-site iteration, and the +//! dissipator branches (Hermitian Pauli fast path vs general complex Pauli +//! sum) — lives in the [`ppvm_lindblad`] crate. This module is responsible +//! only for the Python boundary: decoding the `(N, n_qubits)` numpy uint8 +//! arrays into [`ppvm_lindblad::Word`] vectors, and re-encoding outputs +//! back into numpy. + +use std::collections::HashMap; + +use num::Complex; +use numpy::{IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2}; +use ppvm_lindblad::{JumpInput, LindbladSpec as CoreSpec, Word, codes_from_word, word_from_codes}; +use pyo3::{exceptions::PyValueError, prelude::*}; + +type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); +type PyCoo<'py> = ( + Bound<'py, PyArray1>, + Bound<'py, PyArray1>, + Bound<'py, PyArray1>, +); + +fn map_err(e: ppvm_lindblad::Error) -> PyErr { + PyValueError::new_err(e.to_string()) +} + +/// Reject a basis that contains the same Pauli word at two distinct rows. +/// Duplicate rows would silently overwrite each other in the generator's +/// row-index map and produce an incorrect sparse matrix. +fn assert_basis_unique(basis: &[Word]) -> PyResult<()> { + let mut seen: HashMap<&Word, usize> = HashMap::with_capacity(basis.len()); + for (i, w) in basis.iter().enumerate() { + if let Some(prev) = seen.insert(w, i) { + return Err(PyValueError::new_err(format!( + "basis contains duplicate Pauli word at row {prev} and row {i}" + ))); + } + } + Ok(()) +} + +/// Decode a `(N, n_qubits)` uint8 ndarray view into `N` packed [`Word`]s. +pub(crate) fn decode_basis( + view: &numpy::ndarray::ArrayView2, + n_qubits: usize, +) -> PyResult> { + let n_basis = view.shape()[0]; + let n_cols = view.shape()[1]; + if n_cols != n_qubits { + return Err(PyValueError::new_err(format!( + "basis has {n_cols} columns but spec.n_qubits = {n_qubits}" + ))); + } + let mut out = Vec::with_capacity(n_basis); + let mut row_buf = vec![0u8; n_qubits]; + for i in 0..n_basis { + let row = view.row(i); + for (q, slot) in row_buf.iter_mut().enumerate() { + *slot = row[q]; + } + out.push(word_from_codes(&row_buf).map_err(map_err)?); + } + Ok(out) +} + +/// Pack `Vec<(Word, f64)>` into the standard PyO3 return shape. +fn pack_pauli_map<'py>( + py: Python<'py>, + pairs: Vec<(Word, f64)>, + n_qubits: usize, +) -> PyResult> { + let m = pairs.len(); + let mut basis = vec![0u8; m * n_qubits]; + let mut coeffs = vec![0f64; m]; + for (i, (w, c)) in pairs.into_iter().enumerate() { + codes_from_word(&w, &mut basis[i * n_qubits..(i + 1) * n_qubits]); + coeffs[i] = c; + } + let basis_arr = basis + .into_pyarray(py) + .reshape([m, n_qubits]) + .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + Ok((basis_arr, coeffs.into_pyarray(py))) +} + +/// PyO3 facade exposing [`ppvm_lindblad::LindbladSpec`] to Python. +#[pyclass] +pub struct LindbladSpec { + inner: CoreSpec, +} + +#[pymethods] +impl LindbladSpec { + /// Construct a Lindbladian spec from Hamiltonian terms and jump operators. + /// + /// `jump_lincombs[k]` is a list of `(pauli_string, real, imag)` triples + /// encoding `L_k = Σ_a (re + i·im) P_a`. A length-1 jump with `im == 0` + /// is routed to the Hermitian-Pauli fast path (with rate scaled by `re²`). + #[new] + #[pyo3(signature = (n_qubits, h_terms, h_coeffs, jump_lincombs, jump_rates))] + fn new( + n_qubits: usize, + h_terms: Vec, + h_coeffs: Vec, + jump_lincombs: Vec>, + jump_rates: Vec, + ) -> PyResult { + if h_terms.len() != h_coeffs.len() { + return Err(PyValueError::new_err( + "h_terms and h_coeffs must have the same length", + )); + } + if jump_lincombs.len() != jump_rates.len() { + return Err(PyValueError::new_err( + "jump_lincombs and jump_rates must have the same length", + )); + } + let h: Vec<(String, f64)> = h_terms.into_iter().zip(h_coeffs).collect(); + let jumps: Vec = jump_lincombs + .into_iter() + .zip(jump_rates) + .map(|(lincomb, rate)| JumpInput { + lincomb: lincomb + .into_iter() + .map(|(s, re, im)| (s, Complex::new(re, im))) + .collect(), + rate, + }) + .collect(); + let inner = CoreSpec::new(n_qubits, &h, &jumps).map_err(map_err)?; + Ok(Self { inner }) + } + + #[getter] + fn n_qubits(&self) -> usize { + self.inner.n_qubits() + } + + #[getter] + fn num_h_terms(&self) -> usize { + self.inner.num_h_terms() + } + + #[getter] + fn num_jump_terms(&self) -> usize { + self.inner.num_jump_terms() + } + + /// Apply `L*` to a single Pauli string `p`. + fn action<'py>( + &self, + py: Python<'py>, + p: PyReadonlyArray1<'py, u8>, + ) -> PyResult> { + let p_slice = p.as_slice()?; + let p_word = word_from_codes(p_slice).map_err(map_err)?; + let pairs = self.inner.action(&p_word); + pack_pauli_map(py, pairs, self.inner.n_qubits()) + } + + /// Off-basis component of `L*( Σ_j coeffs[j] · basis[j] )`. + #[pyo3(signature = (basis, coeffs, protected = None))] + fn leakage<'py>( + &self, + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, f64>, + protected: Option>, + ) -> PyResult> { + let n_q = self.inner.n_qubits(); + let basis_view = basis.as_array(); + let basis_words = decode_basis(&basis_view, n_q)?; + let coeffs_slice = coeffs.as_slice()?; + if coeffs_slice.len() != basis_words.len() { + return Err(PyValueError::new_err(format!( + "coeffs has length {} but basis has {} rows", + coeffs_slice.len(), + basis_words.len() + ))); + } + let protected_words: Vec = if let Some(ref prot) = protected { + let pv = prot.as_array(); + decode_basis(&pv, n_q)? + } else { + Vec::new() + }; + let pairs = self + .inner + .leakage(&basis_words, coeffs_slice, &protected_words) + .map_err(map_err)?; + pack_pauli_map(py, pairs, n_q) + } + + /// One predictor-corrector adaptive step. + /// + /// Internally: expand basis with first-hop leakage, predictor step + /// (`exp(dt·M)`), expand again with second-hop leakage from the + /// predicted state, then redo the step from the pre-step coefficients + /// on the doubly-enlarged basis. The matrix exponential is computed in + /// Rust via `quspin-expm`; no scipy required. + /// + /// Returns `(new_basis, new_coeffs)`. + /// + /// `max_basis` is a hard rank cap on the retained basis; `admit_basis` + /// (when `> max_basis`) bounds in-step enrichment instead, so the final + /// cap selects the top-`max_basis` strings over the whole union + /// (displacement truncation). `drop_tol` prunes by magnitude after the + /// step; `tau_add` filters leakage admission by inflow rate. Protected + /// words are never dropped. + #[pyo3(signature = ( + basis, coeffs, dt, max_basis, + drop_tol = 0.0, + protected = None, + num_threads = None, + admit_basis = None, + tau_add = None, + ))] + #[allow(clippy::too_many_arguments)] + fn pc_step<'py>( + &self, + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, f64>, + dt: f64, + max_basis: usize, + drop_tol: f64, + protected: Option>, + num_threads: Option, + admit_basis: Option, + tau_add: Option, + ) -> PyResult> { + let n_q = self.inner.n_qubits(); + let basis_view = basis.as_array(); + let mut basis_words = decode_basis(&basis_view, n_q)?; + assert_basis_unique(&basis_words)?; + let mut coeffs_vec = coeffs.as_slice()?.to_vec(); + if coeffs_vec.len() != basis_words.len() { + return Err(PyValueError::new_err(format!( + "coeffs has length {} but basis has {} rows", + coeffs_vec.len(), + basis_words.len() + ))); + } + let protected_words: Vec = if let Some(ref p) = protected { + decode_basis(&p.as_array(), n_q)? + } else { + Vec::new() + }; + self.inner + .pc_step( + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads, + }, + ) + .map_err(map_err)?; + + // Pack output. Basis may have grown; coeffs has the same new length. + let pairs: Vec<(Word, f64)> = basis_words.into_iter().zip(coeffs_vec).collect(); + pack_pauli_map(py, pairs, n_q) + } + + /// Same as [`Self::pc_step`] but also returns a dict mapping phase + /// name → microseconds spent in that phase, for profiling. + #[pyo3(signature = ( + basis, coeffs, dt, max_basis, + drop_tol = 0.0, + protected = None, + num_threads = None, + admit_basis = None, + tau_add = None, + ))] + #[allow(clippy::too_many_arguments)] + fn pc_step_timed<'py>( + &self, + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, f64>, + dt: f64, + max_basis: usize, + drop_tol: f64, + protected: Option>, + num_threads: Option, + admit_basis: Option, + tau_add: Option, + ) -> PyResult<(PyPauliMap<'py>, Bound<'py, pyo3::types::PyDict>)> { + let n_q = self.inner.n_qubits(); + let basis_view = basis.as_array(); + let mut basis_words = decode_basis(&basis_view, n_q)?; + assert_basis_unique(&basis_words)?; + let mut coeffs_vec = coeffs.as_slice()?.to_vec(); + if coeffs_vec.len() != basis_words.len() { + return Err(PyValueError::new_err(format!( + "coeffs has length {} but basis has {} rows", + coeffs_vec.len(), + basis_words.len() + ))); + } + let protected_words: Vec = if let Some(ref p) = protected { + decode_basis(&p.as_array(), n_q)? + } else { + Vec::new() + }; + let timings = self + .inner + .pc_step_timed( + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads, + }, + ) + .map_err(map_err)?; + + let pairs: Vec<(Word, f64)> = basis_words.into_iter().zip(coeffs_vec).collect(); + let map = pack_pauli_map(py, pairs, n_q)?; + let d = pyo3::types::PyDict::new(py); + d.set_item("leakage1_us", timings.leakage1_us)?; + d.set_item("expand1_us", timings.expand1_us)?; + d.set_item("expm1_us", timings.expm1_us)?; + d.set_item("leakage2_us", timings.leakage2_us)?; + d.set_item("expand2_us", timings.expand2_us)?; + d.set_item("expm2_us", timings.expm2_us)?; + Ok((map, d)) + } + + /// Sparse generator matrix in COO form: `(rows, cols, vals)`. + fn generator<'py>( + &self, + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + ) -> PyResult> { + let n_q = self.inner.n_qubits(); + let basis_view = basis.as_array(); + let basis_words = decode_basis(&basis_view, n_q)?; + assert_basis_unique(&basis_words)?; + let triplets = self.inner.generator(&basis_words); + let total = triplets.len(); + let mut rows = Vec::with_capacity(total); + let mut cols = Vec::with_capacity(total); + let mut vals = Vec::with_capacity(total); + for (r, c, v) in triplets { + rows.push(r as u64); + cols.push(c as u64); + vals.push(v); + } + Ok(( + rows.into_pyarray(py), + cols.into_pyarray(py), + vals.into_pyarray(py), + )) + } +} diff --git a/ppvm-python/demo/lindblad_adaptive.py b/ppvm-python/demo/lindblad_adaptive.py new file mode 100644 index 000000000..e6ee9dc1a --- /dev/null +++ b/ppvm-python/demo/lindblad_adaptive.py @@ -0,0 +1,313 @@ +# --- +# jupyter: +# jupytext: +# cell_metadata_filter: -all +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.19.1 +# kernelspec: +# display_name: ppvm (3.12.12) +# language: python +# name: python3 +# --- + +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +# %% [markdown] +# # Adaptive Pauli-Lindbladian time evolution +# +# Direct Heisenberg-picture evolution of a transport observable on a growing +# Pauli-string basis, without Trotterisation. Composes the three +# `ppvm.Lindbladian` primitives — `leakage`, `generator`, and the matrix +# exponential — into a predictor-corrector integrator on an all-to-all XY +# model (Kac-normalised $1/r^\alpha$ couplings) with single-site Z dephasing. +# +# The observable is kept as a finite sum over Pauli strings (the *basis*), +# and each step: +# +# 1. measures the leakage $(\mathbf 1 - P_B)\,\mathcal L^\dagger(\mathcal O)$ +# out of the current basis $B$; +# 2. adds the leakage strings above `add_tol` to $B$; +# 3. advances the coefficient vector by $\exp(dt\,M)$, where +# $M = P_B \mathcal L^\dagger P_B$ is the Lindbladian restricted to $B$; +# 4. prunes coefficients below `drop_tol`; the protected set (the target +# observable's own support) is never dropped. +# +# The matrix exponential is exact in `dt` within the basis, so there is no +# Trotter splitting error; the only approximation is the finite basis, +# monitored by the cumulative discarded weight. +# +# With `predictor_corrector=True`, the predicted state is fed back as a +# second leakage probe — enriching $B$ with the strings the predictor flows +# into before re-running the step from the pre-step state. This lifts the +# per-step adaptive-integration error from $\mathcal O(dt^2)$ to $\mathcal O(dt^3)$. +# +# The observable is the single-Z Fourier mode +# $\mathcal O_k = \sum_j \cos(k x_j) Z_j$; its decay +# $C_k(t)/C_k(0)$ is the (infinite-temperature) spin transport coefficient. + +# %% +import matplotlib.pyplot as plt +import numpy as np +import scipy.sparse as sp +from scipy.sparse.linalg import expm_multiply + +from ppvm import Lindbladian + +# %% [markdown] +# ## Parameters + +# %% +L = 8 +alpha = 3.0 +gamma = 0.1 +dt = 0.05 +steps = 20 +kmax = 3 +add_tol = 1e-8 +drop_tol = 1e-10 +max_pauli_weight = L +max_basis = 0 # 0 = no cap +predictor_corrector = True + +times = np.arange(steps + 1) * dt +k_indices = np.arange(1, kmax + 1) +k_modes = 2 * np.pi * k_indices / L +x = (np.arange(L) - L // 2 + L // 2) % L - L // 2 + + +# %% [markdown] +# ## Model +# +# All-to-all XY Hamiltonian $H = \sum_{a drop_tol: + rows.append(zterm_codes(j)) + coeffs.append(c) + basis_arr = np.array(rows, dtype=np.uint8) + coeff = np.array(coeffs, dtype=float) + index = {row.tobytes(): i for i, row in enumerate(basis_arr)} + protected_keys = set(index) + return basis_arr, coeff, index, protected_keys + + +def add_leakage_to_basis(basis_arr, probe_coeff, extend_coeffs, protected_arr, index): + """Compute leakage from `probe_coeff`, add above-threshold strings to + `basis_arr`, and pad each vector in `extend_coeffs` with zeros for the + new rows. Returns ``(basis_arr, [extended...], index, rate_below)`` where + `rate_below` is the l2 norm of the leakage too small to add.""" + leak_basis, leak_coeffs = L_op.leakage_arr(basis_arr, probe_coeff, protected_arr) + if not len(leak_coeffs): + return basis_arr, list(extend_coeffs), index, 0.0 + weights = weights_of(leak_basis) + add_mask = (np.abs(leak_coeffs) > add_tol) & (weights <= max_pauli_weight) + rate_below = float(np.linalg.norm(leak_coeffs[~add_mask])) + if not add_mask.any(): + return basis_arr, list(extend_coeffs), index, rate_below + cand = leak_basis[add_mask] + cand = cand[np.argsort(np.abs(leak_coeffs[add_mask]))[::-1]] + cand = cand[np.array([row.tobytes() not in index for row in cand])] + if max_basis: + cand = cand[: max(max_basis - len(basis_arr), 0)] + if not len(cand): + return basis_arr, list(extend_coeffs), index, rate_below + n0 = len(basis_arr) + basis_arr = np.vstack([basis_arr, cand]) + extended = [np.r_[c, np.zeros(len(cand))] for c in extend_coeffs] + for i, row in enumerate(cand): + index[row.tobytes()] = n0 + i + return basis_arr, extended, index, rate_below + + +def prune_and_cap(basis_arr, coeff, protected_keys): + """Drop below-`drop_tol` coefficients; if `max_basis` is set, cap by + keeping protected rows + largest-magnitude others.""" + keep = np.array( + [(row.tobytes() in protected_keys) or (abs(v) >= drop_tol) + for row, v in zip(basis_arr, coeff)] + ) + basis_arr = basis_arr[keep] + coeff = coeff[keep] + keys = [row.tobytes() for row in basis_arr] + index = {kk: i for i, kk in enumerate(keys)} + protected_keys = {pk for pk in protected_keys if pk in index} + if max_basis and len(basis_arr) > max_basis: + slots = max(max_basis - len(protected_keys), 0) + is_protected = np.array([kk in protected_keys for kk in keys]) + order = sorted(np.where(~is_protected)[0].tolist(), + key=lambda i: abs(coeff[i]), reverse=True) + keep2 = is_protected.copy() + keep2[order[:slots]] = True + basis_arr = basis_arr[keep2] + coeff = coeff[keep2] + keys = [row.tobytes() for row in basis_arr] + index = {kk: i for i, kk in enumerate(keys)} + protected_keys = {pk for pk in protected_keys if pk in index} + return basis_arr, coeff, index, protected_keys + + +# %% [markdown] +# ## Run one $k$-mode +# +# Returns per-step $C_k(t)$, basis size, max weight, and cumulative discarded +# weight (trapezoidal a-posteriori error estimate normalised by +# $\|\mathcal O_k\|$). + +# %% +def run_mode(kk): + basis_arr, coeff, index, protected_keys = init_mode(kk) + # `basis_arr` rebinds to a new ndarray on every vstack/slice; nothing + # mutates the initial array in place, so the target/protected views can + # simply alias it without a copy. + target_arr = basis_arr + target_coeff = coeff + protected_arr = basis_arr + norm0 = float(np.dot(coeff, coeff)) + norm_target = np.sqrt(norm0) + + ck = np.empty(steps + 1) + n_basis_t = np.empty(steps + 1, dtype=np.int64) + max_w_t = np.empty(steps + 1, dtype=np.int64) + discarded_cum = np.zeros(steps + 1) + + for nt in range(steps + 1): + # Overlap with the (fixed) initial target. + c_t = sum(coeff[index[trow.tobytes()]] * tc + for trow, tc in zip(target_arr, target_coeff) + if trow.tobytes() in index) + ck[nt] = c_t / norm0 + n_basis_t[nt] = len(basis_arr) + max_w_t[nt] = int(weights_of(basis_arr).max()) + if nt == steps: + break + + # Predictor: enrich basis with leakage from current state, then + # advance by exp(dt · M). + basis_arr, [coeff], index, rate_before = add_leakage_to_basis( + basis_arr, coeff, [coeff], protected_arr, index + ) + coeff_pre = coeff.copy() + coeff = expm_multiply(dt * generator_sparse(L_op, basis_arr), coeff) + + if predictor_corrector: + # Probe leakage with the predicted state; extend the pre-step + # vector with zeros for the new rows and re-run on the enlarged + # basis. Lifts O(dt²) -> O(dt³). + basis_arr, [coeff_pre], index, _ = add_leakage_to_basis( + basis_arr, coeff, [coeff_pre], protected_arr, index + ) + coeff = expm_multiply(dt * generator_sparse(L_op, basis_arr), coeff_pre) + + # Post-step leakage rate, for the trapezoidal error estimate. + norm2_pre_prune = float(np.dot(coeff, coeff)) + _, leak_after = L_op.leakage_arr(basis_arr, coeff, protected_arr) + rate_after = float(np.linalg.norm(leak_after)) if len(leak_after) else 0.0 + + basis_arr, coeff, index, protected_keys = prune_and_cap( + basis_arr, coeff, protected_keys + ) + # Pruning only removes entries, so the discarded l2 weight is + # sqrt(‖c‖²_pre − ‖c‖²_post). + dropped_w = np.sqrt(max(norm2_pre_prune - float(np.dot(coeff, coeff)), 0.0)) + d_total = 0.5 * dt * (rate_before + rate_after) + dropped_w + discarded_cum[nt + 1] = discarded_cum[nt] + d_total / norm_target + + return ck, n_basis_t, max_w_t, discarded_cum + + +# %% [markdown] +# ## Run all $k$-modes and plot + +# %% +Ck = np.empty((steps + 1, kmax)) +n_basis = np.empty((steps + 1, kmax), dtype=np.int64) +max_weight = np.empty((steps + 1, kmax), dtype=np.int64) +discarded_cum = np.empty((steps + 1, kmax)) +for m, kk in enumerate(k_modes): + Ck[:, m], n_basis[:, m], max_weight[:, m], discarded_cum[:, m] = run_mode(kk) + +# %% +fig, ax = plt.subplots() +for m in range(kmax): + ax.plot(times, Ck[:, m], "o-", ms=3, label=rf"$k = 2\pi\cdot{k_indices[m]}/L$") +ax.set_xlabel("$t$") +ax.set_ylabel(r"$C_k(t)/C_k(0)$") +ax.set_title(f"Adaptive Lindbladian evolution L={L} γ={gamma} α={alpha}") +ax.legend() +plt.tight_layout() +plt.show() + +# %% +fig, ax = plt.subplots(1, 2, figsize=(10, 4)) +for m in range(kmax): + ax[0].plot(times, n_basis[:, m], "o-", ms=3, label=rf"$k_{k_indices[m]}$") + ax[1].semilogy(times, discarded_cum[:, m] + 1e-16, "o-", ms=3, + label=rf"$k_{k_indices[m]}$") +ax[0].set(xlabel="$t$", ylabel="|basis|", title="Basis-size growth") +ax[1].set(xlabel="$t$", ylabel=r"cum. discarded / $||O_k||$", + title="A-posteriori error estimate") +for a in ax: + a.legend() +plt.tight_layout() +plt.show() diff --git a/ppvm-python/demo/lindblad_pc_scaling.py b/ppvm-python/demo/lindblad_pc_scaling.py new file mode 100644 index 000000000..528442155 --- /dev/null +++ b/ppvm-python/demo/lindblad_pc_scaling.py @@ -0,0 +1,185 @@ +# --- +# jupyter: +# jupytext: +# cell_metadata_filter: -all +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.19.1 +# kernelspec: +# display_name: ppvm (3.12.12) +# language: python +# name: python3 +# --- + +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +# %% [markdown] +# # `pc_step` parallel scaling +# +# End-to-end wall-time scaling of the pure-Rust predictor-corrector step +# (`ppvm.Lindbladian.pc_step`) with rayon thread count. The entire `pc_step` +# body — both leakage calls, the action-cache build, and both matrix +# exponentials — runs inside a rayon pool of the requested size: +# +# * leakage and the action cache parallelise over basis elements; +# * the matrix exponential parallelises over SpMV columns. +# +# So the speedup numbers reflect overall PC throughput, not just SpMV. + +# %% +from statistics import median +import time + +import matplotlib.pyplot as plt +import numpy as np + +from ppvm import Lindbladian + + +# %% [markdown] +# ## Parameters + +# %% +L = 51 +J = 1.0 +gamma = 1.0 +alpha = 1.0 +dt = 0.05 +n_steps = 20 +max_basis = 10_000_000 # large: rank cap never binds (full enrichment) +max_cores = 4 +warmup_steps = 4 +model = "long-range" # "nn" or "long-range" + + +# %% [markdown] +# ## Model +# +# All-to-all XY with $1/r^\alpha$ couplings (Kac-normalised) and per-site Z +# dephasing. Long-range activates every bond every step, giving a basis +# size that meaningfully exercises parallel scaling. + +# %% +def build_nn_xy_dephasing(L, J, gamma): + h_terms = [] + for i in range(L - 1): + a, b = i, i + 1 + xs = ["I"] * L + xs[a] = xs[b] = "X" + ys = ["I"] * L + ys[a] = ys[b] = "Y" + h_terms += [("".join(xs), J), ("".join(ys), J)] + jump_terms = [("I" * j + "Z" + "I" * (L - j - 1), gamma) for j in range(L)] + return h_terms, jump_terms + + +def build_long_range_xy_dephasing(L, J, alpha, gamma): + pairs = [ + (a, b, 1.0 / min(b - a, L - b + a) ** alpha) + for a in range(L) + for b in range(a + 1, L) + ] + kac = sum(j for _, _, j in pairs) / L + pairs = [(a, b, j / kac) for a, b, j in pairs] + h_terms = [] + for a, b, j in pairs: + for q in "XY": + term = ["I"] * L + term[a] = term[b] = q + h_terms.append(("".join(term), J * j)) + jump_terms = [("I" * j + "Z" + "I" * (L - j - 1), gamma) for j in range(L)] + return h_terms, jump_terms + + +if model == "nn": + h_terms, jump_terms = build_nn_xy_dephasing(L, J, gamma) +else: + h_terms, jump_terms = build_long_range_xy_dephasing(L, J, alpha, gamma) +L_op = Lindbladian(L, h_terms, jump_terms) + + +# %% [markdown] +# ## Timing harness +# +# Each call to `run_pc_steps` runs `n_steps` consecutive PC steps from +# $Z_{L/2}$, returning the per-step wall times and the final basis size. +# The `num_threads` kwarg pins this call to a freshly-built rayon pool of +# that size, isolating thread-count effects from JIT cache state. + +# %% +def run_pc_steps(L_op, L, site0, dt, n_steps, max_basis, num_threads): + z_strings = ["I" * j + "Z" + "I" * (L - j - 1) for j in range(L)] + basis = [z_strings[site0]] + coeffs = np.array([1.0]) + protected = [z_strings[site0]] + times = [] + for _ in range(n_steps): + t0 = time.perf_counter() + basis, coeffs = L_op.pc_step( + basis, + coeffs, + dt, + max_basis, + protected=protected, + num_threads=num_threads, + ) + times.append(time.perf_counter() - t0) + return times, len(basis) + + +# %% [markdown] +# ## Warmup +# +# Each thread count pre-builds its rayon pool and amortises one-time setup +# before the real timing pass. + +# %% +site0 = L // 2 +for n in range(1, max_cores + 1): + run_pc_steps(L_op, L, site0, dt, warmup_steps, max_basis, n) + + +# %% [markdown] +# ## Scaling sweep + +# %% +results = [] +for n in range(1, max_cores + 1): + times, basis_size = run_pc_steps(L_op, L, site0, dt, n_steps, max_basis, n) + first = times[0] * 1000.0 + steady = median(times[1:]) * 1000.0 + results.append({"threads": n, "first_ms": first, "steady_ms": steady, + "basis": basis_size}) + +baseline = results[0]["steady_ms"] +print(f"{'threads':>8s} {'first-step (ms)':>16s} {'steady (ms)':>12s} " + f"{'speedup':>9s} {'|basis|':>8s}") +for r in results: + speedup = baseline / r["steady_ms"] + print(f"{r['threads']:>8d} {r['first_ms']:>16.1f} {r['steady_ms']:>12.2f} " + f"{speedup:>8.2f}x {r['basis']:>8d}") + + +# %% [markdown] +# ## Plot +# +# Steady-state speedup vs thread count, with the linear-scaling reference. + +# %% +threads = np.array([r["threads"] for r in results]) +steady_ms = np.array([r["steady_ms"] for r in results]) +speedup = steady_ms[0] / steady_ms + +fig, ax = plt.subplots() +ax.plot(threads, speedup, "o-", label="measured") +ax.plot(threads, threads, "k--", alpha=0.4, label="linear") +ax.set_xlabel("threads") +ax.set_ylabel("speedup (vs 1 thread)") +ax.set_title(f"pc_step parallel scaling L={L} model={model} " + f"|basis|={results[-1]['basis']}") +ax.legend() +plt.tight_layout() +plt.show() diff --git a/ppvm-python/pyproject.toml b/ppvm-python/pyproject.toml index 5189852df..9afd25b59 100644 --- a/ppvm-python/pyproject.toml +++ b/ppvm-python/pyproject.toml @@ -10,6 +10,7 @@ requires-python = ">=3.10" dependencies = [ "bloqade-circuit>=0.14.1", "kirin-toolchain~=0.22.2", + "numpy>=1.26", ] [build-system] @@ -58,3 +59,8 @@ dev = [ "pytest>=9.0.2", "pytest-benchmark>=5.2.3", ] +# Optional: only the `demo/` scripts use it (`expm_multiply` for the reference +# matrix exponential). Runtime ppvm and the test suite have no scipy dep. +demo = [ + "scipy>=1.13", +] diff --git a/ppvm-python/src/ppvm/__init__.py b/ppvm-python/src/ppvm/__init__.py index 8eaff9490..286250622 100644 --- a/ppvm-python/src/ppvm/__init__.py +++ b/ppvm-python/src/ppvm/__init__.py @@ -5,6 +5,9 @@ from .generalized_tableau import GeneralizedTableau as GeneralizedTableau from .generalized_tableau import MeasurementResult as MeasurementResult from .generalized_tableau import sample_stim as sample_stim +from .lindblad import Lindbladian as Lindbladian +from .lindblad import sigma_minus as sigma_minus +from .lindblad import sigma_plus as sigma_plus from .paulisum import LossyPauliSum as LossyPauliSum from .paulisum import PauliSum as PauliSum from .squin_interpreter.device import ( diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index bd890adc7..4b2cfac7f 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -1,5 +1,7 @@ from collections.abc import Sequence +import numpy as np + class _PauliSumBase: def __init__( self, @@ -343,3 +345,51 @@ class TableauSumSampler29(_TableauSumSamplerBase): ... class TableauSumSampler30(_TableauSumSamplerBase): ... class TableauSumSampler31(_TableauSumSamplerBase): ... class TableauSumSampler32(_TableauSumSamplerBase): ... + +class LindbladSpec: + def __init__( + self, + n_qubits: int, + h_terms: list[str], + h_coeffs: list[float], + jump_lincombs: list[list[tuple[str, float, float]]], + jump_rates: list[float], + ) -> None: ... + @property + def n_qubits(self) -> int: ... + @property + def num_h_terms(self) -> int: ... + @property + def num_jump_terms(self) -> int: ... + def action(self, p: np.ndarray) -> tuple[np.ndarray, np.ndarray]: ... + def leakage( + self, + basis: np.ndarray, + coeffs: np.ndarray, + protected: np.ndarray | None = None, + ) -> tuple[np.ndarray, np.ndarray]: ... + def pc_step( + self, + basis: np.ndarray, + coeffs: np.ndarray, + dt: float, + max_basis: int, + drop_tol: float = 0.0, + protected: np.ndarray | None = None, + num_threads: int | None = None, + admit_basis: int | None = None, + tau_add: float | None = None, + ) -> tuple[np.ndarray, np.ndarray]: ... + def pc_step_timed( + self, + basis: np.ndarray, + coeffs: np.ndarray, + dt: float, + max_basis: int, + drop_tol: float = 0.0, + protected: np.ndarray | None = None, + num_threads: int | None = None, + admit_basis: int | None = None, + tau_add: float | None = None, + ) -> tuple[tuple[np.ndarray, np.ndarray], dict[str, int]]: ... + def generator(self, basis: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: ... diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py new file mode 100644 index 000000000..a31ed4e5f --- /dev/null +++ b/ppvm-python/src/ppvm/lindblad.py @@ -0,0 +1,372 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +"""Direct Pauli-Lindbladian time evolution on an adaptive Pauli-string basis. + +Given a Hermitian Pauli Hamiltonian H = Σ c_i P_i and jump operators +L_k = Σ_a λ_{k,a} P_{k,a} (each a complex linear combination of Pauli +strings) with rates γ_k ≥ 0, this module exposes the primitives for +adaptive Heisenberg-picture evolution: + +- ``pc_step(...)`` / ``pc_step_arr(...)``: one adaptive predictor-corrector + step ``O ← exp(dt·L*) O`` +- ``action(p)`` / ``action_arr(p)``: L*(p) for one Pauli string p +- ``leakage(basis, coeffs)`` / ``leakage_arr(...)``: off-basis component of + L*(Σ c_j p_j), driving basis expansion +- ``generator(basis)``: COO triples ``(rows, cols, vals)`` for the generator + matrix M such that L* restricted to ``basis`` is ``M @ coeffs``. Users + wanting a sparse matrix can wrap them — e.g. + ``scipy.sparse.coo_matrix((vals, (rows, cols)), shape=(N, N)).tocsc()`` + +The ``*_arr`` variants pass Pauli strings as ``(N, n_qubits)`` ``uint8`` +arrays of Pauli codes (``0=I, 1=X, 2=Z, 3=Y``) and skip string +construction entirely — at ~10^5 basis rows per evolution step, per-row +``str.join`` dominates wall time. + +Each jump term can be either: + +- a single Hermitian Pauli (`("ZZII", γ)`), routed to a fast diagonal path, + or +- a complex Pauli sum (`([("XIII", 0.5+0j), ("YIII", 0+0.5j)], γ)`) to + describe e.g. amplitude-damping (`σ⁻`) and excitation (`σ⁺`) operators. + +For the general case the dissipator +``γ ( L† p L − ½ {L†L, p} )`` is evaluated directly; the L†L Pauli +expansion is precomputed once at construction. +""" + +from __future__ import annotations + +from collections.abc import Iterable, Sequence +from typing import Union + +import numpy as np + +from ._core import LindbladSpec as _LindbladSpec + +_PAULI_CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} +# Lookup table mapping code -> ASCII byte for vectorised string output. +_CODE_TO_ASCII = np.array([ord("I"), ord("X"), ord("Z"), ord("Y")], dtype=np.uint8) + +# A jump operator is either a Hermitian Pauli (single string) or a complex +# linear combination of Pauli strings. +PauliLincomb = Iterable[tuple[str, complex]] +JumpSpec = Union[tuple[str, float], tuple[PauliLincomb, float]] + + +def _string_to_codes(s: str, n_qubits: int) -> np.ndarray: + """Encode a Pauli string ``"IXYZ..."`` as a length-``n_qubits`` uint8 array. + + Underscores in the input are ignored, matching the Rust parser + (``parse_pauli_string`` in `ppvm-lindblad`) so users can write + ``"X_Y_Z"`` for readability. + """ + s_clean = s.replace("_", "") + if len(s_clean) != n_qubits: + raise ValueError( + f"Pauli string {s!r} has length {len(s_clean)} (after stripping '_') " + f"!= n_qubits {n_qubits}" + ) + try: + return np.array([_PAULI_CODE[c] for c in s_clean], dtype=np.uint8) + except KeyError as exc: + bad = exc.args[0] + raise ValueError( + f"Pauli string {s!r} contains invalid character {bad!r}; " + f"expected one of 'I', 'X', 'Y', 'Z' (and '_' is allowed as a separator)" + ) from None + + +def _codes_to_string(codes: np.ndarray) -> str: + """Decode one length-``n_qubits`` row of Pauli codes back to a string.""" + return _CODE_TO_ASCII[codes].tobytes().decode("ascii") + + +def _basis_to_codes(basis: Sequence[str], n_qubits: int) -> np.ndarray: + """Stack a sequence of Pauli strings into an ``(N, n_qubits)`` uint8 array.""" + arr = np.zeros((len(basis), n_qubits), dtype=np.uint8) + for i, s in enumerate(basis): + arr[i] = _string_to_codes(s, n_qubits) + return arr + + +def _codes_to_basis(arr: np.ndarray) -> list[str]: + """Inverse of `_basis_to_codes`. One call into C per row.""" + bytes_per_row = _CODE_TO_ASCII[arr].tobytes() + n = arr.shape[1] + return [bytes_per_row[i * n : (i + 1) * n].decode("ascii") for i in range(arr.shape[0])] + + +def sigma_plus(site: int, n_qubits: int) -> list[tuple[str, complex]]: + """``σ⁺_q = (X_q + i Y_q) / 2``. Use as a Lindblad jump for excitation.""" + if not 0 <= site < n_qubits: + raise ValueError(f"site {site} out of range for n_qubits={n_qubits}") + x_str = "I" * site + "X" + "I" * (n_qubits - site - 1) + y_str = "I" * site + "Y" + "I" * (n_qubits - site - 1) + return [(x_str, 0.5 + 0.0j), (y_str, 0.0 + 0.5j)] + + +def sigma_minus(site: int, n_qubits: int) -> list[tuple[str, complex]]: + """``σ⁻_q = (X_q − i Y_q) / 2``. Use as a Lindblad jump for amplitude damping.""" + if not 0 <= site < n_qubits: + raise ValueError(f"site {site} out of range for n_qubits={n_qubits}") + x_str = "I" * site + "X" + "I" * (n_qubits - site - 1) + y_str = "I" * site + "Y" + "I" * (n_qubits - site - 1) + return [(x_str, 0.5 + 0.0j), (y_str, 0.0 - 0.5j)] + + +def _normalize_jump(jump_op: str | PauliLincomb) -> list[tuple[str, float, float]]: + """Convert a user-supplied jump operator to ``[(pauli_str, re, im), ...]``. + + Accepts either a single Pauli string (treated as a Hermitian-Pauli jump + with coefficient 1) or an iterable of ``(pauli_str, complex_coeff)`` + pairs. + """ + if isinstance(jump_op, str): + return [(jump_op, 1.0, 0.0)] + out: list[tuple[str, float, float]] = [] + for term in jump_op: + s, c = term + cc = complex(c) + out.append((str(s), float(cc.real), float(cc.imag))) + if not out: + raise ValueError("jump operator lincomb must contain at least one Pauli term") + return out + + +class Lindbladian: + """Pre-compiled adjoint Pauli-Lindbladian acting on Pauli strings. + + Parameters + ---------- + n_qubits: + Number of qubits. + h_terms: + Iterable of ``(pauli_string, coefficient)`` pairs for the + Hermitian Hamiltonian ``H = Σ c_i P_i``. Each ``pauli_string`` is + a length-``n_qubits`` ``str`` over ``"IXYZ"``. + jump_terms: + Iterable of ``(jump_op, rate)`` pairs. ``jump_op`` is either a + Pauli string ``"XYZI..."`` (treated as a Hermitian-Pauli jump + with coefficient 1, hitting the fast path) or an iterable of + ``(pauli_string, complex_coeff)`` pairs for a general complex + Pauli linear combination such as `sigma_plus` or + `sigma_minus`. ``rate`` is the non-negative GKSL rate + ``γ_k``. + + Examples + -------- + Dephasing (Hermitian Pauli): + + >>> Lindbladian(2, [("XX", 1.0)], [("ZI", 0.3), ("IZ", 0.3)]) + + Amplitude damping on site 0 (non-Hermitian): + + >>> jumps = [(sigma_minus(0, 2), 0.5)] + >>> Lindbladian(2, [("XX", 1.0)], jumps) + """ + + def __init__( + self, + n_qubits: int, + h_terms: Iterable[tuple[str, float]], + jump_terms: Iterable[tuple[str | PauliLincomb, float]] = (), + ): + self.n_qubits = int(n_qubits) + h_strs: list[str] = [] + h_coeffs: list[float] = [] + for s, c in h_terms: + h_strs.append(s) + h_coeffs.append(float(c)) + j_lincombs: list[list[tuple[str, float, float]]] = [] + j_rates: list[float] = [] + for jump_op, rate in jump_terms: + j_lincombs.append(_normalize_jump(jump_op)) + j_rates.append(float(rate)) + self._spec = _LindbladSpec(self.n_qubits, h_strs, h_coeffs, j_lincombs, j_rates) + + @property + def num_h_terms(self) -> int: + return self._spec.num_h_terms + + @property + def num_jump_terms(self) -> int: + return self._spec.num_jump_terms + + # ── Pure-ndarray hot path ── + + def action_arr(self, p: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Apply ``L*`` to a single Pauli string given as uint8 codes. + + Returns ``(out_basis, out_coeffs)``: a ``(M, n_qubits)`` uint8 + array and a length-``M`` float64 array. + """ + return self._spec.action(np.ascontiguousarray(p, dtype=np.uint8)) + + def leakage_arr( + self, + basis_arr: np.ndarray, + coeffs: np.ndarray, + protected_arr: np.ndarray | None = None, + ) -> tuple[np.ndarray, np.ndarray]: + """Off-basis component of ``L*( Σ_j coeffs[j] basis[j] )``. + + ``basis_arr``: ``(N, n_qubits)`` uint8. ``coeffs``: length-N float64. + ``protected_arr``: optional ``(K, n_qubits)`` uint8 of Pauli strings + that must NEVER appear in the leakage output. + + Returns ``(out_basis, out_coeffs)`` packed the same way as + `action_arr`. + """ + n = self.n_qubits + if protected_arr is None: + protected_arr = np.zeros((0, n), dtype=np.uint8) + return self._spec.leakage( + np.ascontiguousarray(basis_arr, dtype=np.uint8), + np.ascontiguousarray(coeffs, dtype=np.float64), + np.ascontiguousarray(protected_arr, dtype=np.uint8), + ) + + def pc_step_arr( + self, + basis_arr: np.ndarray, + coeffs: np.ndarray, + dt: float, + max_basis: int, + drop_tol: float = 1e-12, + protected_arr: np.ndarray | None = None, + num_threads: int | None = None, + admit_basis: int | None = None, + tau_add: float | None = None, + ) -> tuple[np.ndarray, np.ndarray]: + """One predictor-corrector adaptive step. + + All work — leakage expansion, matrix-exponential step, second-hop + re-expansion, corrector — runs in Rust; SciPy is not required. + The matrix-exponential action is computed matrix-free via the + external ``quspin-expm`` crate. + + Truncation. ``max_basis`` is a hard rank cap on the live basis: + enrichment adds at most ``max_basis - len(basis)`` of the largest + leakage strings, and the post-step basis is trimmed to the + top-``max_basis`` entries by ``|coeff|`` (``protected`` words always + kept). Pass a large value (e.g. ``10_000_000``) for the near-exact, + uncapped case. ``drop_tol`` additionally prunes basis entries whose + absolute coefficient is below the threshold after the corrector + (unless the word is ``protected``). + + ``num_threads``, when set, pins this call to a freshly-built rayon + pool of that size — useful for benchmarking parallel scaling. + + ``tau_add``, when set, filters leakage admission by an absolute + coefficient-rate threshold (a dt- and drop_tol-independent + parameterization). Off by default — with cap-based truncation it + is at most a modest wall optimization. + + ``admit_basis``, when set (must be >= ``max_basis``), bounds the + enriched working set during the step instead of ``max_basis``: the + step may hold up to ``admit_basis`` strings transiently, and the + final truncation keeps the top-``max_basis`` by evolved ``|coeff|`` + over the whole union (retained + newly admitted) — rank + displacement, so no ``drop_tol`` is needed to sustain membership + turnover. With the default ``None``, admission is bounded by + ``max_basis`` and turnover requires ``drop_tol > 0``. + + Returns ``(new_basis_arr, new_coeffs)``; the basis may have grown + (or shrunk, if ``max_basis`` / ``drop_tol`` pruned entries). + """ + n = self.n_qubits + if protected_arr is None: + protected_arr = np.zeros((0, n), dtype=np.uint8) + return self._spec.pc_step( + np.ascontiguousarray(basis_arr, dtype=np.uint8), + np.ascontiguousarray(coeffs, dtype=np.float64), + float(dt), + int(max_basis), + float(drop_tol), + np.ascontiguousarray(protected_arr, dtype=np.uint8), + None if num_threads is None else int(num_threads), + None if admit_basis is None else int(admit_basis), + None if tau_add is None else float(tau_add), + ) + + def pc_step( + self, + basis: Sequence[str], + coeffs: np.ndarray, + dt: float, + max_basis: int, + drop_tol: float = 1e-12, + protected: Sequence[str] | None = None, + num_threads: int | None = None, + admit_basis: int | None = None, + tau_add: float | None = None, + ) -> tuple[list[str], np.ndarray]: + """String-keyed variant of `pc_step_arr`.""" + n = self.n_qubits + basis_arr = _basis_to_codes(basis, n) + protected_arr = ( + _basis_to_codes(list(protected), n) if protected else np.zeros((0, n), dtype=np.uint8) + ) + new_basis_arr, new_coeffs = self.pc_step_arr( + basis_arr, + coeffs, + dt, + max_basis, + drop_tol, + protected_arr, + num_threads, + admit_basis, + tau_add, + ) + return _codes_to_basis(new_basis_arr), new_coeffs + + def generator_arr(self, basis_arr: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Generator matrix as COO triples ``(rows, cols, vals)``. + + Basis given as uint8 codes. To get a SciPy sparse matrix: + + >>> import scipy.sparse as sp + >>> rows, cols, vals = L_op.generator_arr(basis_arr) + >>> M = sp.coo_matrix( + ... (vals, (rows, cols)), shape=(len(basis_arr), len(basis_arr)) + ... ).tocsc() + """ + return self._spec.generator(np.ascontiguousarray(basis_arr, dtype=np.uint8)) + + # ── String-keyed convenience API (slower; for tests / display) ── + + def action(self, p: str) -> dict[str, float]: + """Apply ``L*`` to a single Pauli string ``p`` (string-keyed dict).""" + codes = _string_to_codes(p, self.n_qubits) + out_basis, out_coeffs = self._spec.action(codes) + keys = _codes_to_basis(out_basis) + return {k: float(v) for k, v in zip(keys, out_coeffs) if v != 0.0} + + def leakage( + self, + basis: Sequence[str], + coeffs: np.ndarray, + protected: Sequence[str] | None = None, + ) -> dict[str, float]: + """Off-basis leakage as a ``dict[str, float]`` (slower API).""" + n = self.n_qubits + basis_arr = _basis_to_codes(basis, n) + protected_arr = ( + _basis_to_codes(list(protected), n) if protected else np.zeros((0, n), dtype=np.uint8) + ) + out_basis, out_coeffs = self._spec.leakage( + basis_arr, + np.ascontiguousarray(coeffs, dtype=np.float64), + protected_arr, + ) + keys = _codes_to_basis(out_basis) + return {k: float(v) for k, v in zip(keys, out_coeffs) if v != 0.0} + + def generator(self, basis: Sequence[str]) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Generator matrix as COO triples ``(rows, cols, vals)``, + basis given as strings. See `generator_arr` for the conversion + to a SciPy sparse matrix.""" + n = self.n_qubits + basis_arr = _basis_to_codes(basis, n) + return self.generator_arr(basis_arr) diff --git a/ppvm-python/test/lindblad/__init__.py b/ppvm-python/test/lindblad/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/ppvm-python/test/lindblad/_helpers.py b/ppvm-python/test/lindblad/_helpers.py new file mode 100644 index 000000000..25484772e --- /dev/null +++ b/ppvm-python/test/lindblad/_helpers.py @@ -0,0 +1,356 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Shared helpers for the Lindbladian tests. + +Two reference kernels live here: + +- :func:`_reference_action` builds `L*(p)` for **Hermitian-Pauli** jumps via + the single-qubit Pauli multiplication table (cheap; only depends on `p`'s + weight). +- :func:`_dense_action` builds `L*(p)` for **arbitrary jump operators** by + constructing the full 2^L × 2^L dense Liouvillian. Only viable for L ≤ 3, + but it makes no assumption about the shape of the jumps. + +The bilinear NN-XY + Z-dephasing reference :func:`_bilinear_nn_xy_z_dephasing_obc` +gives an exact closed-form answer that the predictor-corrector tests converge +toward as dt → 0. +""" + +from __future__ import annotations + +import numpy as np + +# --- Pauli multiplication table for the Hermitian-Pauli reference ----------- +I, X, Y, Z = range(4) # noqa: E741 (I is standard Pauli notation here) +MUL = { + (I, I): (1, I), + (I, X): (1, X), + (I, Y): (1, Y), + (I, Z): (1, Z), + (X, I): (1, X), + (X, X): (1, I), + (X, Y): (1j, Z), + (X, Z): (-1j, Y), + (Y, I): (1, Y), + (Y, X): (-1j, Z), + (Y, Y): (1, I), + (Y, Z): (1j, X), + (Z, I): (1, Z), + (Z, X): (1j, Y), + (Z, Y): (-1j, X), + (Z, Z): (1, I), +} +CODE = {I: "I", X: "X", Y: "Y", Z: "Z"} +ICODE = {"I": I, "X": X, "Y": Y, "Z": Z} + + +def _str_to_tuple(s): + return tuple(ICODE[ch] for ch in s) + + +def _tuple_to_str(t): + return "".join(CODE[ch] for ch in t) + + +def _mul_pauli(p, q): + """Return ``(phase, p·q)`` for two Pauli strings given as code tuples.""" + phase = 1 + r = list(p) + for i, (pi, qi) in enumerate(zip(p, q)): + ph, rr = MUL[(pi, qi)] + phase *= ph + r[i] = rr + return phase, tuple(r) + + +def _reference_action(p_str, h_terms, jump_terms): + """``L*(p) = i[H, p] + sum_k gamma_k (L_k p L_k - p)``, term by term.""" + p = _str_to_tuple(p_str) + out: dict = {} + for h_str, coeff_h in h_terms: + h = _str_to_tuple(h_str) + ph_pp, pp = _mul_pauli(h, p) + ph_pq, _ = _mul_pauli(p, h) + # i[H, p] = i (Hp - pH) = i (ph_pp - ph_pq) r ; real for Hermitian H, p. + coeff = (1j * coeff_h * (ph_pp - ph_pq)).real + if coeff: + out[pp] = out.get(pp, 0.0) + coeff + for j_str, gamma in jump_terms: + j = _str_to_tuple(j_str) + ph_pp, _ = _mul_pauli(j, p) + # Hermitian Pauli L, p: Lp has imaginary phase iff they anti-commute, + # and then L p L = -p, contributing -2 gamma p. + if abs(ph_pp.imag) > 0.5: + out[p] = out.get(p, 0.0) + (-2.0 * gamma) + return {kk: v for kk, v in out.items() if v} + + +def _to_str_dict(d): + return {_tuple_to_str(kk): v for kk, v in d.items()} + + +# --- model builders --------------------------------------------------------- + + +def xy_dephasing(L, alpha, gamma): + """Long-range XY model with PBC + Z dephasing.""" + pairs = [ + (a, b, 1.0 / min(b - a, L - b + a) ** alpha) for a in range(L) for b in range(a + 1, L) + ] + kac = sum(j for _, _, j in pairs) / L + pairs = [(a, b, j / kac) for a, b, j in pairs] + h_terms = [] + for a, b, j in pairs: + for q in "XY": + term = ["I"] * L + term[a] = term[b] = q + h_terms.append(("".join(term), j)) + jump_terms = [("I" * i + "Z" + "I" * (L - i - 1), gamma) for i in range(L)] + return h_terms, jump_terms + + +def tfim_xdeph(L, J, h, gamma): + """TFIM (ZZ + transverse X) with X dephasing.""" + h_terms = [] + for i in range(L - 1): + term = ["I"] * L + term[i] = term[i + 1] = "Z" + h_terms.append(("".join(term), J)) + for i in range(L): + term = ["I"] * L + term[i] = "X" + h_terms.append(("".join(term), h)) + jump_terms = [("I" * i + "X" + "I" * (L - i - 1), gamma) for i in range(L)] + return h_terms, jump_terms + + +def nn_xy_z_dephasing_obc(L, J, gamma): + """Nearest-neighbour XY (OBC) + per-site Z dephasing.""" + h_terms = [] + for i in range(L - 1): + a, b = i, i + 1 + xs = ["I"] * L + xs[a] = xs[b] = "X" + ys = ["I"] * L + ys[a] = ys[b] = "Y" + h_terms += [("".join(xs), J), ("".join(ys), J)] + jump_terms = [("I" * j + "Z" + "I" * (L - j - 1), gamma) for j in range(L)] + return h_terms, jump_terms + + +def random_pauli_str(rng, L): + chars = ["I"] * L + positions = rng.choice(L, size=int(rng.integers(1, L + 1)), replace=False) + for q in positions: + chars[q] = rng.choice(["X", "Y", "Z"]) + return "".join(chars) + + +def assert_action_matches(L_op, h_terms, jump_terms, strings): + """Compare `L_op.action(p)` against :func:`_reference_action` for each `p`.""" + for p in strings: + got = L_op.action(p) + want = _to_str_dict(_reference_action(p, h_terms, jump_terms)) + for kk in set(got) | set(want): + assert abs(got.get(kk, 0.0) - want.get(kk, 0.0)) < 1e-12, ( + f"action mismatch at p={p!r} k={kk!r}: " + f"shim={got.get(kk, 0.0)} ref={want.get(kk, 0.0)}" + ) + + +# --- dense Liouvillian reference (for non-Hermitian jumps) ------------------ + +_DENSE_PAULI = { + "I": np.eye(2, dtype=complex), + "X": np.array([[0, 1], [1, 0]], dtype=complex), + "Y": np.array([[0, -1j], [1j, 0]], dtype=complex), + "Z": np.array([[1, 0], [0, -1]], dtype=complex), +} +SIGMA_MINUS_MAT = np.array([[0, 0], [1, 0]], dtype=complex) +SIGMA_PLUS_MAT = np.array([[0, 1], [0, 0]], dtype=complex) + + +def pauli_mat(s): + """Dense matrix for Pauli string ``s`` (leftmost char = leftmost factor).""" + M = np.array([[1.0]], dtype=complex) + for c in s: + M = np.kron(M, _DENSE_PAULI[c]) + return M + + +def all_strings(L): + if L == 0: + return [""] + sub = all_strings(L - 1) + return [c + s for c in "IXYZ" for s in sub] + + +def dense_action(H, jumps, p_str, L): + """`L*(p)` computed densely, returned as a real Pauli-coefficient dict.""" + p_mat = pauli_mat(p_str) + # macOS Accelerate emits spurious divide warnings on exact zeros. + with np.errstate(divide="ignore", invalid="ignore", over="ignore"): + out_mat = 1j * (H @ p_mat - p_mat @ H) + for Lop, gamma in jumps: + Ld = Lop.conj().T + out_mat += gamma * (Ld @ p_mat @ Lop - 0.5 * (Ld @ Lop @ p_mat + p_mat @ Ld @ Lop)) + d = 2**L + out = {} + for q_str in all_strings(L): + coef = np.trace(pauli_mat(q_str) @ out_mat) / d + assert abs(coef.imag) < 1e-9, f"non-real coefficient for {q_str}: {coef}" + if abs(coef.real) > 1e-12: + out[q_str] = coef.real + return out + + +def embed_op(op_1q, site, L): + """Embed a single-qubit dense operator at ``site`` of an L-qubit register.""" + eye = _DENSE_PAULI["I"] + M = np.array([[1.0]], dtype=complex) + for j in range(L): + M = np.kron(M, op_1q if j == site else eye) + return M + + +# --- closed bilinear evolution for NN-XY + Z-dephasing (OBC) ---------------- + + +def bilinear_nn_xy_z_dephasing_obc(L, J, gamma, times, site0): + """Closed bilinear evolution of `C_j(t)` for the NN XY + Z-dephasing model + with open boundary conditions. + + `F_{mn}(t) = 2^{-L} Tr[B_{mn}(t) Z_i]` evolves as + `∂_t F_{mn} = i·2J·(F_{m+1,n}+F_{m-1,n}-F_{m,n+1}-F_{m,n-1}) - 4γ(1-δ_{mn})F_{mn}` + on the L×L lattice; OBC means edge terms (m=0 or m=L-1, etc.) drop. + Z_j = I - 2 n_j gives C_j = -2 F_{jj}. + """ + dim = L * L + + def idx(m, n): + return m * L + n + + gen = np.zeros((dim, dim), dtype=complex) + for m in range(L): + for n in range(L): + row = idx(m, n) + if m + 1 < L: + gen[row, idx(m + 1, n)] += 1j * 2 * J + if m - 1 >= 0: + gen[row, idx(m - 1, n)] += 1j * 2 * J + if n + 1 < L: + gen[row, idx(m, n + 1)] += -1j * 2 * J + if n - 1 >= 0: + gen[row, idx(m, n - 1)] += -1j * 2 * J + if m != n: + gen[row, row] += -4 * gamma + + evals, evecs = np.linalg.eig(gen) + evecs_inv = np.linalg.inv(evecs) + f0 = np.zeros(dim, dtype=complex) + f0[idx(site0, site0)] = -0.5 + coeffs = evecs_inv @ f0 + + corr = np.empty((len(times), L)) + diag = [idx(j, j) for j in range(L)] + for nt, t in enumerate(times): + ft = evecs @ (np.exp(evals * t) * coeffs) + corr[nt] = np.real(-2 * ft[diag]) + return corr + + +# --- numpy-only matrix exponential reference -------------------------------- + + +def coo_to_dense(triples, n_basis): + """Build a dense ``(n_basis, n_basis)`` array from COO triples.""" + rows, cols, vals = triples + M = np.zeros((n_basis, n_basis), dtype=float) + M[rows, cols] = vals + return M + + +def expm_mv_dense(M, v): + """``exp(M) @ v`` via numpy eigendecomposition. Independent of the Rust + Al-Mohy & Higham implementation; small bases only. + + The Lindbladian is generally diagonalizable, so + ``M = V diag(λ) V^{-1}`` and ``exp(M) v = V (exp(λ) ⊙ (V^{-1} v))``. + """ + evals, evecs = np.linalg.eig(M) + rhs = np.linalg.solve(evecs, v.astype(complex)) + return np.real(evecs @ (np.exp(evals) * rhs)) + + +# --- adaptive PC evolution reference (used by test_adaptive_pc.py) ---------- + + +def _generator_dense(L_op, basis): + """`L_op.generator(basis)` returns COO triples; convert to dense.""" + return coo_to_dense(L_op.generator(basis), len(basis)) + + +def adaptive_z_correlator(L_op, L, site0, dt, n_steps, tau_add): + """Adaptive Heisenberg-picture evolution of Z_{site0} on a growing basis. + + First-hop only: each step adds the strings from `L_op.leakage(...)`, + then matrix-exponentiates the (small, dense) restricted generator via + numpy eigendecomposition. Local truncation is O(dt²) per step. + """ + z_strings = ["I" * j + "Z" + "I" * (L - j - 1) for j in range(L)] + basis = [z_strings[site0]] + coeffs = np.array([1.0]) + protected = [z_strings[site0]] + + corr = np.zeros((n_steps + 1, L)) + corr[0, site0] = 1.0 + + for step in range(n_steps): + leak = L_op.leakage(basis, coeffs, protected=protected) + new = [k for k, v in leak.items() if abs(v) > tau_add] + if new: + basis = basis + new + coeffs = np.concatenate([coeffs, np.zeros(len(new))]) + M = _generator_dense(L_op, basis) + coeffs = expm_mv_dense(dt * M, coeffs) + index = {s: i for i, s in enumerate(basis)} + for j in range(L): + if z_strings[j] in index: + corr[step + 1, j] = coeffs[index[z_strings[j]]] + return corr + + +def adaptive_z_correlator_pc(L_op, L, site0, dt, n_steps, tau_add): + """Same as :func:`adaptive_z_correlator` but with predictor-corrector + basis expansion: predict, then add the second-hop leakage strings before + redoing the step. Lifts the per-step truncation from O(dt²) to O(dt³). + """ + z_strings = ["I" * j + "Z" + "I" * (L - j - 1) for j in range(L)] + basis = [z_strings[site0]] + coeffs = np.array([1.0]) + protected = [z_strings[site0]] + + corr = np.zeros((n_steps + 1, L)) + corr[0, site0] = 1.0 + + for step in range(n_steps): + leak = L_op.leakage(basis, coeffs, protected=protected) + new = [k for k, v in leak.items() if abs(v) > tau_add] + if new: + basis = basis + new + coeffs = np.concatenate([coeffs, np.zeros(len(new))]) + coeffs_pre = coeffs.copy() + M = _generator_dense(L_op, basis) + coeffs_predict = expm_mv_dense(dt * M, coeffs) + leak2 = L_op.leakage(basis, coeffs_predict, protected=protected) + new2 = [k for k, v in leak2.items() if abs(v) > tau_add] + if new2: + basis = basis + new2 + coeffs_pre = np.concatenate([coeffs_pre, np.zeros(len(new2))]) + M = _generator_dense(L_op, basis) + coeffs = expm_mv_dense(dt * M, coeffs_pre) + index = {s: i for i, s in enumerate(basis)} + for j in range(L): + if z_strings[j] in index: + corr[step + 1, j] = coeffs[index[z_strings[j]]] + return corr diff --git a/ppvm-python/test/lindblad/test_action_generator.py b/ppvm-python/test/lindblad/test_action_generator.py new file mode 100644 index 000000000..7b571c9c9 --- /dev/null +++ b/ppvm-python/test/lindblad/test_action_generator.py @@ -0,0 +1,108 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Cross-checks for :meth:`Lindbladian.action` / :meth:`Lindbladian.generator` +/ :meth:`Lindbladian.leakage` against the Hermitian-Pauli reference built from +the single-qubit Pauli multiplication table. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from ppvm import Lindbladian + +from ._helpers import ( + _reference_action, + _tuple_to_str, + assert_action_matches, + coo_to_dense, + expm_mv_dense, + random_pauli_str, + tfim_xdeph, + xy_dephasing, +) + + +def test_action_xy_dephasing(): + L = 8 + h_terms, jump_terms = xy_dephasing(L, alpha=1.0, gamma=0.3) + L_op = Lindbladian(L, h_terms, jump_terms) + rng = np.random.default_rng(42) + strings = ["I" * L] + strings += ["I" * i + "Z" + "I" * (L - i - 1) for i in range(L)] + strings += [random_pauli_str(rng, L) for _ in range(20)] + assert_action_matches(L_op, h_terms, jump_terms, strings) + + +def test_action_tfim_xdephasing(): + L = 6 + h_terms, jump_terms = tfim_xdeph(L, J=0.7, h=0.4, gamma=0.2) + L_op = Lindbladian(L, h_terms, jump_terms) + rng = np.random.default_rng(7) + strings = [random_pauli_str(rng, L) for _ in range(30)] + assert_action_matches(L_op, h_terms, jump_terms, strings) + + +def test_generator_leakage_and_expm(): + L = 5 + dt = 0.1 + h_terms, jump_terms = xy_dephasing(L, alpha=1.0, gamma=0.4) + L_op = Lindbladian(L, h_terms, jump_terms) + basis = ["I" * i + "Z" + "I" * (L - i - 1) for i in range(L)] + basis += ["YIYII", "IYIYI", "ZZIII"] + coeffs = np.array([0.5, -0.3, 0.2, 0.4, -0.1, 0.6, 0.2, 0.1]) + + M_shim = coo_to_dense(L_op.generator(basis), len(basis)) + M_ref = np.zeros((len(basis), len(basis))) + index = {p: i for i, p in enumerate(basis)} + for col, p in enumerate(basis): + for r_tuple, v in _reference_action(p, h_terms, jump_terms).items(): + r = _tuple_to_str(r_tuple) + if r in index: + M_ref[index[r], col] += v + assert np.max(np.abs(M_shim - M_ref)) < 1e-12 + + shim_leak = L_op.leakage(basis, coeffs) + ref_leak: dict = {} + for p, cf in zip(basis, coeffs): + for r_tuple, v in _reference_action(p, h_terms, jump_terms).items(): + r = _tuple_to_str(r_tuple) + if r not in index: + ref_leak[r] = ref_leak.get(r, 0.0) + v * cf + ref_leak = {kk: v for kk, v in ref_leak.items() if v} + for kk in set(shim_leak) | set(ref_leak): + assert abs(shim_leak.get(kk, 0.0) - ref_leak.get(kk, 0.0)) < 1e-12 + + c_shim = expm_mv_dense(dt * M_shim, coeffs) + c_ref = expm_mv_dense(dt * M_ref, coeffs) + assert np.allclose(c_shim, c_ref, atol=1e-13) + + +def test_generator_rejects_duplicate_basis(): + """Duplicate basis rows would silently overwrite each other in the + row-index map and produce an incorrect sparse generator. The user-facing + entry point must reject them with a clear ValueError instead. + """ + L = 4 + h_terms, jump_terms = xy_dephasing(L, alpha=1.0, gamma=0.3) + L_op = Lindbladian(L, h_terms, jump_terms) + basis = ["ZIII", "IZII", "ZIII"] # duplicate at rows 0 and 2 + with pytest.raises(ValueError, match=r"duplicate Pauli word at row 0 and row 2"): + L_op.generator(basis) + # pc_step also builds the row index and must reject duplicates. + with pytest.raises(ValueError, match=r"duplicate Pauli word"): + L_op.pc_step(basis, np.ones(len(basis)), 0.01, 10_000_000) + + +def test_protected_strings_suppressed(): + L = 4 + h_terms, jump_terms = xy_dephasing(L, alpha=1.0, gamma=0.0) + L_op = Lindbladian(L, h_terms, jump_terms) + basis = ["ZIII"] + coeffs = np.array([1.0]) + leak = L_op.leakage(basis, coeffs) + assert leak, "expected some leakage" + protected_key = next(iter(leak)) + leak2 = L_op.leakage(basis, coeffs, protected=[protected_key]) + assert protected_key not in leak2 diff --git a/ppvm-python/test/lindblad/test_adaptive_pc.py b/ppvm-python/test/lindblad/test_adaptive_pc.py new file mode 100644 index 000000000..ab467c90c --- /dev/null +++ b/ppvm-python/test/lindblad/test_adaptive_pc.py @@ -0,0 +1,108 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""End-to-end convergence of the adaptive predictor-corrector evolution +(numpy eigendecomposition reference) against the closed bilinear NN-XY + +Z-dephasing solution. + +For NN interactions the JW bilinears stay closed under the adjoint +Lindbladian, so the spin correlator obeys a tractable L²×L² ODE +(:func:`bilinear_nn_xy_z_dephasing_obc`). OBC keeps this exact; PBC would +introduce a parity-twist that only matches up to 1/L corrections. +""" + +from __future__ import annotations + +from itertools import pairwise + +import numpy as np + +from ppvm import Lindbladian + +from ._helpers import ( + adaptive_z_correlator, + adaptive_z_correlator_pc, + bilinear_nn_xy_z_dephasing_obc, + nn_xy_z_dephasing_obc, +) + + +def test_adaptive_converges_to_nn_xy_z_dephasing_bilinear(): + """Halving dt drives the adaptive stepper toward the closed bilinear solution. + + Single-hop has local truncation O(dt²) per step → global error O(T·dt). + """ + L = 4 + J = 1.0 + gamma = 1.0 + site0 = L // 2 + T = 0.05 + tau_add = 1e-12 # tight enough that integrator (T·dt) dominates + + h_terms, jump_terms = nn_xy_z_dephasing_obc(L, J, gamma) + L_op = Lindbladian(L, h_terms, jump_terms) + + errors = [] + final_corr = None + for dt in (0.01, 0.005, 0.0025): + n_steps = round(T / dt) + times = np.arange(n_steps + 1) * dt + exact = bilinear_nn_xy_z_dephasing_obc(L, J, gamma, times, site0) + shim = adaptive_z_correlator(L_op, L, site0, dt, n_steps, tau_add) + # Endpoint comparison only — independent of which step counts we use. + errors.append(float(np.max(np.abs(shim[-1] - exact[-1])))) + if dt == 0.0025: + final_corr = (shim[-1], exact[-1]) + + assert final_corr is not None + + # Halving dt should roughly halve the error; allow 2× slack. + assert errors[1] < 0.8 * errors[0], f"dt-halving did not help: {errors}" + assert errors[2] < 0.8 * errors[1], f"dt-halving did not help: {errors}" + + # Integrator floor at dt = 0.0025 is ~T·dt = 1.25e-4; expect <1e-3. + assert errors[-1] < 1e-3, ( + f"shim vs bilinear at smallest dt: max abs error = {errors[-1]:.3g}; " + f"shim={final_corr[0]}, exact={final_corr[1]}" + ) + + +def test_predictor_corrector_lifts_dt_scaling_to_cubic(): + """The predictor-corrector basis expansion lifts the single-hop scheme's + local O(dt²) truncation to O(dt³). PC error is also strictly smaller at + every dt we test. + """ + L = 4 + J = 1.0 + gamma = 1.0 + site0 = L // 2 + T = 0.05 + tau_add = 1e-12 + + h_terms, jump_terms = nn_xy_z_dephasing_obc(L, J, gamma) + L_op = Lindbladian(L, h_terms, jump_terms) + + err_single = [] + err_pc = [] + for dt in (0.01, 0.005, 0.0025): + n_steps = round(T / dt) + times = np.arange(n_steps + 1) * dt + exact = bilinear_nn_xy_z_dephasing_obc(L, J, gamma, times, site0) + single = adaptive_z_correlator(L_op, L, site0, dt, n_steps, tau_add) + pc = adaptive_z_correlator_pc(L_op, L, site0, dt, n_steps, tau_add) + err_single.append(float(np.max(np.abs(single[-1] - exact[-1])))) + err_pc.append(float(np.max(np.abs(pc[-1] - exact[-1])))) + + # PC strictly more accurate than single-hop at every dt (by ~100× in this + # regime). Threshold loose enough to absorb expm_multiply tolerance noise. + for s, p, dt in zip(err_single, err_pc, (0.01, 0.005, 0.0025)): + assert p < s / 50, ( + f"PC ({p:.3e}) not meaningfully better than single-hop ({s:.3e}) at dt={dt}" + ) + + # dt-scaling one order steeper: halving dt should drop the error ~8× + # (dt³ vs single-hop's ~4×). Require >5× per halving with safety margin. + for prev, curr in pairwise(err_pc): + assert curr < prev / 5, f"PC dt-halving ratio < 5: errors {err_pc}" + + # Smallest-dt PC error should sit at FP noise of the bilinear reference. + assert err_pc[-1] < 1e-7, f"PC at smallest dt: error = {err_pc[-1]:.3e}" diff --git a/ppvm-python/test/lindblad/test_non_hermitian_jumps.py b/ppvm-python/test/lindblad/test_non_hermitian_jumps.py new file mode 100644 index 000000000..e0170a567 --- /dev/null +++ b/ppvm-python/test/lindblad/test_non_hermitian_jumps.py @@ -0,0 +1,106 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Non-Hermitian dissipators: cross-check :meth:`Lindbladian.action` / +:meth:`Lindbladian.generator` / :meth:`Lindbladian.leakage` against the dense +2^L × 2^L Liouvillian reference. Only viable for L ≤ 3. +""" + +from __future__ import annotations + +import numpy as np + +from ppvm import Lindbladian, sigma_minus, sigma_plus + +from ._helpers import ( + SIGMA_MINUS_MAT, + SIGMA_PLUS_MAT, + all_strings, + coo_to_dense, + dense_action, + embed_op, + pauli_mat, + random_pauli_str, +) + + +def test_amplitude_damping_action(): + L = 3 + gamma = 0.5 + h_terms = [("XXI", 1.0), ("IXX", 0.7), ("ZIZ", 0.3)] + jump_terms = [(sigma_minus(i, L), gamma) for i in range(L)] + L_op = Lindbladian(L, h_terms, jump_terms) + + H = sum(c * pauli_mat(s) for s, c in h_terms) + jumps_dense = [(embed_op(SIGMA_MINUS_MAT, i, L), gamma) for i in range(L)] + + rng = np.random.default_rng(11) + strings = ["III", "ZII", "IZI", "IIZ", "XYZ", "YYZ"] + strings += [random_pauli_str(rng, L) for _ in range(10)] + + for p in strings: + got = L_op.action(p) + want = dense_action(H, jumps_dense, p, L) + for k in set(got) | set(want): + diff = abs(got.get(k, 0.0) - want.get(k, 0.0)) + assert diff < 1e-10, ( + f"sigma_minus action mismatch at p={p!r} k={k!r}: " + f"shim={got.get(k, 0.0)} ref={want.get(k, 0.0)} diff={diff}" + ) + + +def test_thermal_excitation_damping_action(): + """σ⁺ + σ⁻ jumps together (thermal bath at finite temperature).""" + L = 2 + h_terms = [("XX", 1.0), ("ZI", 0.2), ("IZ", 0.1)] + jump_terms = [(sigma_minus(i, L), 0.4) for i in range(L)] + [ + (sigma_plus(i, L), 0.1) for i in range(L) + ] + L_op = Lindbladian(L, h_terms, jump_terms) + + H = sum(c * pauli_mat(s) for s, c in h_terms) + jumps_dense = [(embed_op(SIGMA_MINUS_MAT, i, L), 0.4) for i in range(L)] + [ + (embed_op(SIGMA_PLUS_MAT, i, L), 0.1) for i in range(L) + ] + + for p in all_strings(L): + got = L_op.action(p) + want = dense_action(H, jumps_dense, p, L) + for k in set(got) | set(want): + diff = abs(got.get(k, 0.0) - want.get(k, 0.0)) + assert diff < 1e-10, ( + f"sigma_plus/sigma_minus action mismatch at p={p!r} k={k!r}: " + f"shim={got.get(k, 0.0)} ref={want.get(k, 0.0)}" + ) + + +def test_amplitude_damping_generator_and_leakage(): + L = 3 + gamma = 0.3 + h_terms = [("XXI", 0.5), ("IXX", 0.5), ("ZII", 0.2), ("IZI", 0.2), ("IIZ", 0.2)] + jump_terms = [(sigma_minus(0, L), gamma), (sigma_minus(2, L), gamma)] + L_op = Lindbladian(L, h_terms, jump_terms) + + basis = ["III", "ZII", "IZI", "IIZ", "ZZI", "IZZ"] + coeffs = np.array([0.1, 0.5, -0.3, 0.4, 0.2, -0.1]) + + H = sum(c * pauli_mat(s) for s, c in h_terms) + jumps_dense = [(embed_op(SIGMA_MINUS_MAT, i, L), gamma) for i in (0, 2)] + + M_shim = coo_to_dense(L_op.generator(basis), len(basis)) + M_ref = np.zeros((len(basis), len(basis))) + idx = {p: i for i, p in enumerate(basis)} + leak_ref: dict = {} + for col, p in enumerate(basis): + action_p = dense_action(H, jumps_dense, p, L) + for q, v in action_p.items(): + if q in idx: + M_ref[idx[q], col] += v + else: + leak_ref[q] = leak_ref.get(q, 0.0) + v * coeffs[col] + assert np.max(np.abs(M_shim - M_ref)) < 1e-10 + + leak_shim = L_op.leakage(basis, coeffs) + leak_ref = {k: v for k, v in leak_ref.items() if abs(v) > 1e-14} + for k in set(leak_shim) | set(leak_ref): + diff = abs(leak_shim.get(k, 0.0) - leak_ref.get(k, 0.0)) + assert diff < 1e-10, f"leakage mismatch at k={k!r}: diff={diff}" diff --git a/ppvm-python/test/lindblad/test_pc_step_rust.py b/ppvm-python/test/lindblad/test_pc_step_rust.py new file mode 100644 index 000000000..f2269e474 --- /dev/null +++ b/ppvm-python/test/lindblad/test_pc_step_rust.py @@ -0,0 +1,123 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Pure-Rust :meth:`Lindbladian.pc_step` (Al-Mohy & Higham expm + parallel +SpMV): agrees with the numpy-eigendecomp PC reference at FP precision and +shows the same cubic dt-scaling against the bilinear reference. Plus a +sanity check that a length-1 real lincomb routes to the Hermitian fast path. +""" + +from __future__ import annotations + +from itertools import pairwise + +import numpy as np + +from ppvm import Lindbladian + +from ._helpers import ( + adaptive_z_correlator_pc, + bilinear_nn_xy_z_dephasing_obc, + nn_xy_z_dephasing_obc, + random_pauli_str, + xy_dephasing, +) + + +def _adaptive_z_correlator_pc_rust(L_op, L, site0, dt, n_steps, max_basis): + """Same as :func:`_helpers.adaptive_z_correlator_pc` but the per-step PC + work (leakage expansion, predictor expm, second-hop expansion, corrector + expm) all runs in Rust through :meth:`Lindbladian.pc_step`.""" + z_strings = ["I" * j + "Z" + "I" * (L - j - 1) for j in range(L)] + basis = [z_strings[site0]] + coeffs = np.array([1.0]) + protected = [z_strings[site0]] + + corr = np.zeros((n_steps + 1, L)) + corr[0, site0] = 1.0 + + for step in range(n_steps): + basis, coeffs = L_op.pc_step( + basis, coeffs, dt, max_basis, drop_tol=0.0, protected=protected + ) + index = {s: i for i, s in enumerate(basis)} + for j in range(L): + if z_strings[j] in index: + corr[step + 1, j] = coeffs[index[z_strings[j]]] + return corr + + +def test_pc_step_rust_matches_python_pc(): + """The pure-Rust PC step agrees with the numpy-eigendecomp PC reference + at FP precision. + + Pins the Rust matrix exponential (Al-Mohy & Higham) against an + independent reference (numpy ``eig``) under the exact same + basis-expansion schedule.""" + L = 4 + J = 1.0 + gamma = 1.0 + site0 = L // 2 + dt = 0.01 + n_steps = 5 + tau_add = 1e-12 + # Large max_basis so the rust rank cap never binds: full enrichment, + # matching the python reference's effectively-all-leakage tau_add. + max_basis = 10_000_000 + + h_terms, jump_terms = nn_xy_z_dephasing_obc(L, J, gamma) + L_op = Lindbladian(L, h_terms, jump_terms) + + rust = _adaptive_z_correlator_pc_rust(L_op, L, site0, dt, n_steps, max_basis) + py_ref = adaptive_z_correlator_pc(L_op, L, site0, dt, n_steps, tau_add) + + diff = float(np.max(np.abs(rust - py_ref))) + assert diff < 1e-10, f"Rust PC differs from numpy-eigendecomp PC by {diff:.3e}" + + +def test_pc_step_rust_dt_scaling_is_cubic(): + """End-to-end: the Rust-only PC step matches the bilinear reference with + cubic dt-scaling, confirming the Rust matrix exponential is not the + accuracy bottleneck.""" + L = 4 + J = 1.0 + gamma = 1.0 + site0 = L // 2 + T = 0.05 + max_basis = 10_000_000 # large: rank cap never binds (full enrichment) + + h_terms, jump_terms = nn_xy_z_dephasing_obc(L, J, gamma) + L_op = Lindbladian(L, h_terms, jump_terms) + + err = [] + for dt in (0.01, 0.005, 0.0025): + n_steps = round(T / dt) + times = np.arange(n_steps + 1) * dt + exact = bilinear_nn_xy_z_dephasing_obc(L, J, gamma, times, site0) + rust = _adaptive_z_correlator_pc_rust(L_op, L, site0, dt, n_steps, max_basis) + err.append(float(np.max(np.abs(rust[-1] - exact[-1])))) + + # Halving dt should drop the error by ≥5× (cubic gives 8×). + for prev, curr in pairwise(err): + assert curr < prev / 5, f"Rust PC dt-halving ratio < 5: errors {err}" + assert err[-1] < 1e-7, f"Rust PC tight-dt error too large: {err[-1]:.3e}" + + +def test_lincomb_single_term_matches_hermitian_fast_path(): + """A length-1 real lincomb should route to the Hermitian fast path and + produce numerically identical results to passing the string directly.""" + L = 4 + h_terms, jump_simple = xy_dephasing(L, alpha=1.0, gamma=0.3) + L_simple = Lindbladian(L, h_terms, jump_simple) + # Same operator, expressed as a length-1 complex lincomb. + jump_lincomb = [([(s, 1.0 + 0.0j)], g) for s, g in jump_simple] + L_lincomb = Lindbladian(L, h_terms, jump_lincomb) + + rng = np.random.default_rng(99) + for _ in range(5): + p = random_pauli_str(rng, L) + got_a = L_simple.action(p) + got_b = L_lincomb.action(p) + for k in set(got_a) | set(got_b): + assert got_a.get(k, 0.0) == got_b.get(k, 0.0), ( + f"lincomb fast path mismatch at p={p!r}: {got_a} vs {got_b}" + ) diff --git a/ppvm-python/uv.lock b/ppvm-python/uv.lock index 88e4fa704..f8fae91cc 100644 --- a/ppvm-python/uv.lock +++ b/ppvm-python/uv.lock @@ -906,9 +906,14 @@ source = { editable = "." } dependencies = [ { name = "bloqade-circuit" }, { name = "kirin-toolchain" }, + { name = "numpy" }, ] [package.dev-dependencies] +demo = [ + { name = "scipy", version = "1.15.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "scipy", version = "1.17.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, +] dev = [ { name = "numpy" }, { name = "pytest" }, @@ -919,9 +924,11 @@ dev = [ requires-dist = [ { name = "bloqade-circuit", specifier = ">=0.14.1" }, { name = "kirin-toolchain", specifier = "~=0.22.2" }, + { name = "numpy", specifier = ">=1.26" }, ] [package.metadata.requires-dev] +demo = [{ name = "scipy", specifier = ">=1.13" }] dev = [ { name = "numpy", specifier = ">=2.2.6" }, { name = "pytest", specifier = ">=9.0.2" },