diff --git a/Cargo.lock b/Cargo.lock index c0f74a9c0..c8ba3bebc 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -193,7 +193,7 @@ dependencies = [ "petgraph 0.6.5", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -219,7 +219,7 @@ checksum = "9035ad2d096bed7955a320ee7e2230574d28fd3c3a0f186cbea1ff3c7eed5dbb" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -366,7 +366,7 @@ checksum = "f9abbd1bc6865053c427f7198e6af43bfdedc55ab791faed4fbd361d789575ff" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -544,7 +544,7 @@ dependencies = [ "heck 0.5.0", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -905,7 +905,7 @@ dependencies = [ "proc-macro2", "quote", "scratch", - "syn", + "syn 2.0.109", ] [[package]] @@ -919,7 +919,7 @@ dependencies = [ "indexmap 2.12.0", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -937,7 +937,7 @@ dependencies = [ "indexmap 2.12.0", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -961,7 +961,7 @@ dependencies = [ "proc-macro2", "quote", "strsim", - "syn", + "syn 2.0.109", ] [[package]] @@ -972,7 +972,7 @@ checksum = "d38308df82d1080de0afee5d069fa14b0326a88c14f15c5ccda35b4a6c414c81" dependencies = [ "darling_core", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -997,7 +997,7 @@ checksum = "6178a82cf56c836a3ba61a7935cdb1c49bfaa6fa4327cd5bf554a503087de26b" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -1018,7 +1018,7 @@ checksum = "d65d7ce8132b7c0e54497a4d9a55a1c2a0912a0d786cf894472ba818fba45762" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -1029,7 +1029,7 @@ checksum = "ef941ded77d15ca19b40374869ac6000af1c9f2a4c0f3d4c70926287e6364a8f" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -1042,7 +1042,7 @@ dependencies = [ "proc-macro2", "quote", "rustc_version", - "syn", + "syn 2.0.109", ] [[package]] @@ -1062,7 +1062,7 @@ checksum = "bda628edc44c4bb645fbe0f758797143e4e07926f7ebf4e9bdfbd3d2ce621df3" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", "unicode-xid", ] @@ -1111,7 +1111,7 @@ checksum = "97369cbbc041bc366949bc74d34658d6cda5621039731c6310521892a3a20ae0" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -1171,7 +1171,7 @@ dependencies = [ "once_cell", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -1350,7 +1350,7 @@ checksum = "162ee34ebcb7c64a8abebc059ce0fee27c2262618d7b60ed8faf72fef13c3650" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -1445,6 +1445,102 @@ dependencies = [ "stable_deref_trait", ] +[[package]] +name = "glam" +version = "0.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "333928d5eb103c5d4050533cec0384302db6be8ef7d3cebd30ec6a35350353da" + +[[package]] +name = "glam" +version = "0.15.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3abb554f8ee44336b72d522e0a7fe86a29e09f839a36022fa869a7dfe941a54b" + +[[package]] +name = "glam" +version = "0.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4126c0479ccf7e8664c36a2d719f5f2c140fbb4f9090008098d2c291fa5b3f16" + +[[package]] +name = "glam" +version = "0.17.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e01732b97afd8508eee3333a541b9f7610f454bb818669e66e90f5f57c93a776" + +[[package]] +name = "glam" +version = "0.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "525a3e490ba77b8e326fb67d4b44b4bd2f920f44d4cc73ccec50adc68e3bee34" + +[[package]] +name = "glam" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b8509e6791516e81c1a630d0bd7fbac36d2fa8712a9da8662e716b52d5051ca" + +[[package]] +name = "glam" +version = "0.20.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f43e957e744be03f5801a55472f593d43fabdebf25a4585db250f04d86b1675f" + +[[package]] +name = "glam" +version = "0.21.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "518faa5064866338b013ff9b2350dc318e14cc4fcd6cb8206d7e7c9886c98815" + +[[package]] +name = "glam" +version = "0.22.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "12f597d56c1bd55a811a1be189459e8fad2bbc272616375602443bdfb37fa774" + +[[package]] +name = "glam" +version = "0.23.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e4afd9ad95555081e109fe1d21f2a30c691b5f0919c67dfa690a2e1eb6bd51c" + +[[package]] +name = "glam" +version = "0.24.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945" + +[[package]] +name = "glam" +version = "0.25.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3" + +[[package]] +name = "glam" +version = "0.27.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e05e7e6723e3455f4818c7b26e855439f7546cf617ef669d1adedb8669e5cb9" + +[[package]] +name = "glam" +version = "0.28.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "779ae4bf7e8421cf91c0b3b64e7e8b40b862fba4d393f59150042de7c4965a94" + +[[package]] +name = "glam" +version = "0.29.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee" + +[[package]] +name = "glam" +version = "0.30.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bd47b05dddf0005d850e5644cae7f2b14ac3df487979dbfff3b56f20b1a6ae46" + [[package]] name = "glob" version = "0.3.3" @@ -1507,6 +1603,12 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "hermit-abi" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c" + [[package]] name = "hex" version = "0.4.3" @@ -1919,7 +2021,7 @@ checksum = "f365c8de536236cfdebd0ba2130de22acefed18b1fb99c32783b3840aec5fb46" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -2018,7 +2120,7 @@ checksum = "980af8b43c3ad5d8d349ace167ec8170839f753a42d233ba19e08afe1850fa69" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -2041,6 +2143,16 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "lambert_w" +version = "1.2.31" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c567f2087fc83535a312e683b6ed8811395690ef896df7b82966b21b7526580" +dependencies = [ + "num-complex", + "num-traits", +] + [[package]] name = "lazy_static" version = "1.5.0" @@ -2053,6 +2165,18 @@ version = "0.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "09edd9e8b54e49e587e4f6295a7d29c3ea94d469cb40ab8ca70b288248a81db2" +[[package]] +name = "levenberg-marquardt" +version = "0.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8be7a65739a815308eef33a6d8c78e435a7317305d5b0af0c8c465a2d7ac6fc6" +dependencies = [ + "cfg-if", + "nalgebra", + "num-traits", + "rustc_version", +] + [[package]] name = "libbz2-rs-sys" version = "0.2.2" @@ -2183,7 +2307,10 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08" dependencies = [ "autocfg", + "num_cpus", + "once_cell", "rawpointer", + "thread-tree", ] [[package]] @@ -2231,6 +2358,49 @@ dependencies = [ "windows-sys 0.61.2", ] +[[package]] +name = "nalgebra" +version = "0.34.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c4d5b3eff5cd580f93da45e64715e8c20a3996342f1e466599cf7a267a0c2f5f" +dependencies = [ + "approx 0.5.1", + "glam 0.14.0", + "glam 0.15.2", + "glam 0.16.0", + "glam 0.17.3", + "glam 0.18.0", + "glam 0.19.0", + "glam 0.20.5", + "glam 0.21.3", + "glam 0.22.0", + "glam 0.23.0", + "glam 0.24.2", + "glam 0.25.0", + "glam 0.27.0", + "glam 0.28.0", + "glam 0.29.3", + "glam 0.30.9", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.109", +] + [[package]] name = "ndarray" version = "0.16.1" @@ -2309,6 +2479,33 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", + "libm", +] + +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + +[[package]] +name = "numpy" +version = "0.27.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fa24ffc88cf9d43f7269d6b6a0d0a00010924a8cc90604a21ef9c433b66998d" +dependencies = [ + "libc", + "ndarray", + "num-complex", + "num-integer", + "num-traits", + "pyo3", + "pyo3-build-config", + "rustc-hash", ] [[package]] @@ -2347,6 +2544,12 @@ version = "0.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "04744f49eae99ab78e0d5c0b603ab218f515ea8cfe5a456d7629ad883a3b6e7d" +[[package]] +name = "order-stat" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "efa535d5117d3661134dbf1719b6f0ffe06f2375843b13935db186cd094105eb" + [[package]] name = "ordered-float" version = "5.1.0" @@ -2379,7 +2582,7 @@ dependencies = [ "proc-macro2", "proc-macro2-diagnostics", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -2427,6 +2630,7 @@ dependencies = [ "pecos-engines", "pecos-hugr-qis", "pecos-llvm", + "pecos-num", "pecos-phir", "pecos-phir-json", "pecos-programs", @@ -2590,6 +2794,18 @@ dependencies = [ "xz2", ] +[[package]] +name = "pecos-num" +version = "0.1.1" +dependencies = [ + "levenberg-marquardt", + "log", + "nalgebra", + "ndarray", + "peroxide", + "roots", +] + [[package]] name = "pecos-phir" version = "0.1.1" @@ -2766,6 +2982,7 @@ dependencies = [ "inkwell", "libc", "log", + "numpy", "parking_lot", "pecos", "pyo3", @@ -2791,6 +3008,39 @@ version = "2.3.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9b4f627cb1b25917193a259e49bdad08f671f8d9708acfd5fe0a8c1455d87220" +[[package]] +name = "peroxide" +version = "0.40.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "394b064b6069af31d72b1c211181d0962c677446e52c7edd50f3d672c83b7394" +dependencies = [ + "anyhow", + "matrixmultiply", + "order-stat", + "paste", + "peroxide-ad", + "peroxide-num", + "puruspe", + "rand 0.9.2", + "rand_distr", +] + +[[package]] +name = "peroxide-ad" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f6fba8ff3f40b67996f7c745f699babaa3e57ef5c8178ec999daf7eedc51dc8c" +dependencies = [ + "quote", + "syn 1.0.109", +] + +[[package]] +name = "peroxide-num" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6b17ddf7141892147b48b5d0e2a3ab8ec7fcbaa06f186d01118f7c933a77863" + [[package]] name = "pest" version = "2.8.3" @@ -2821,7 +3071,7 @@ dependencies = [ "pest_meta", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3044,7 +3294,7 @@ checksum = "af066a9c399a26e020ada66a034357a868728e72cd426f3adcd35f80d88d88c8" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", "version_check", "yansi", ] @@ -3069,7 +3319,17 @@ checksum = "572f69980fc11dd3c07ab054974330844cac436bacb79a69dfda9c2e5c72cba4" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", +] + +[[package]] +name = "puruspe" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69a530540df16fa0cfb1f5a7c26c4d24544a2a9fb94c60831afa776360782d98" +dependencies = [ + "lambert_w", + "num-complex", ] [[package]] @@ -3118,7 +3378,7 @@ dependencies = [ "proc-macro2", "pyo3-macros-backend", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3131,7 +3391,7 @@ dependencies = [ "proc-macro2", "pyo3-build-config", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3267,6 +3527,16 @@ dependencies = [ "getrandom 0.3.4", ] +[[package]] +name = "rand_distr" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a8615d50dcf34fa31f7ab52692afec947c4dd0ab803cc87cb3b0b4570ff7463" +dependencies = [ + "num-traits", + "rand 0.9.2", +] + [[package]] name = "rawpointer" version = "0.2.1" @@ -3330,7 +3600,7 @@ checksum = "b7186006dcb21920990093f30e3dea63b7d6e977bf1256be20c3563a5db070da" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3465,6 +3735,12 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "roots" +version = "0.0.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "082f11ffa03bbef6c2c6ea6bea1acafaade2fd9050ae0234ab44a2153742b058" + [[package]] name = "rstest" version = "0.26.1" @@ -3490,7 +3766,7 @@ dependencies = [ "regex", "relative-path", "rustc_version", - "syn", + "syn 2.0.109", "unicode-ident", ] @@ -3569,6 +3845,15 @@ version = "1.0.20" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "28d3b2b1366ec20994f1fd18c3c594f05c5dd4bc44d8bb0c1c632c8d6829481f" +[[package]] +name = "safe_arch" +version = "0.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" +dependencies = [ + "bytemuck", +] + [[package]] name = "same-file" version = "1.0.6" @@ -3612,7 +3897,7 @@ dependencies = [ "proc-macro2", "quote", "serde_derive_internals", - "syn", + "syn 2.0.109", ] [[package]] @@ -3695,7 +3980,7 @@ checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3706,7 +3991,7 @@ checksum = "18d26a20a969b9e3fdf2fc2d9f21eda6c40e2de84c9408bb5d3b05d499aae711" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3762,7 +4047,7 @@ dependencies = [ "darling", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3799,6 +4084,19 @@ version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" +[[package]] +name = "simba" +version = "0.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c99284beb21666094ba2b75bbceda012e610f5479dfcc2d6e2426f53197ffd95" +dependencies = [ + "approx 0.5.1", + "num-complex", + "num-traits", + "paste", + "wide", +] + [[package]] name = "simd-adler32" version = "0.3.7" @@ -3892,7 +4190,7 @@ dependencies = [ "heck 0.5.0", "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -3901,6 +4199,17 @@ version = "2.6.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "13c2bddecc57b384dee18652358fb23172facb8a2c51ccc10d74c157bdea3292" +[[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]] name = "syn" version = "2.0.109" @@ -3929,7 +4238,7 @@ checksum = "728a70f3dbaf5bab7f0c4b1ac8d7ae5ea60a4b5549c8a5914361c99147a709d2" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -4009,7 +4318,7 @@ checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -4020,7 +4329,16 @@ checksum = "3ff15c8ecd7de3849db632e14d18d2571fa09dfc5ed93479bc4485c7a517c913" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", +] + +[[package]] +name = "thread-tree" +version = "0.3.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ffbd370cb847953a25954d9f63e14824a36113f8c72eecf6eccef5dc4b45d630" +dependencies = [ + "crossbeam-channel", ] [[package]] @@ -4279,7 +4597,7 @@ checksum = "81383ab64e72a7a8b8e13130c49e3dab29def6d0c7d76a03087b3cf71c5c6903" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -4336,7 +4654,7 @@ checksum = "27a7a9b72ba121f6f1f6c3632b85604cac41aedb5ddc70accbebb6cac83de846" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -4523,7 +4841,7 @@ dependencies = [ "bumpalo", "proc-macro2", "quote", - "syn", + "syn 2.0.109", "wasm-bindgen-shared", ] @@ -4757,7 +5075,7 @@ checksum = "8d57f08c4d8acde5550bcd4b45baa16daba411eb6f715d21dbfc26b535c9a17f" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -4811,6 +5129,16 @@ dependencies = [ "rustls-pki-types", ] +[[package]] +name = "wide" +version = "0.7.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" +dependencies = [ + "bytemuck", + "safe_arch", +] + [[package]] name = "winapi-util" version = "0.1.11" @@ -4841,7 +5169,7 @@ checksum = "053e2e040ab57b9dc951b72c264860db7eb3b0200ba345b4e4c3b14f67855ddf" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -4852,7 +5180,7 @@ checksum = "3f316c4a2570ba26bbec722032c4099d8c8bc095efccdc15688708623367e358" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -5118,7 +5446,7 @@ checksum = "b659052874eb698efe5b9e8cf382204678a0086ebf46982b79d6ca3182927e5d" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", "synstructure", ] @@ -5139,7 +5467,7 @@ checksum = "88d2b8d9c68ad2b9e4340d7832716a4d21a22a1154777ad56ea55c51a9cf3831" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] @@ -5159,7 +5487,7 @@ checksum = "d71e5d6e06ab090c67b5e44993ec16b72dcbaabc526db883a360057678b48502" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", "synstructure", ] @@ -5199,7 +5527,7 @@ checksum = "eadce39539ca5cb3985590102671f2567e659fca9666581ad3411d59207951f3" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.109", ] [[package]] diff --git a/Cargo.toml b/Cargo.toml index 42cda71b6..08d0842de 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -85,6 +85,12 @@ tket = "0.16" tket-qsystem = { version = "0.22", default-features = false } ndarray = "0.16" anyhow = "1" + +# Numerical computing dependencies (for pecos-num) +peroxide = "0.40" +roots = "0.0.8" +levenberg-marquardt = "0.15" +nalgebra = "0.34" cxx = "1.0.187" cxx-build = "1.0.187" reqwest = { version = "0.12", default-features = false, features = ["blocking", "rustls-tls"] } @@ -126,6 +132,7 @@ pecos-rslib = { version = "0.1.1", path = "python/pecos-rslib/rust" } pecos-wasm = { version = "0.1.1", path = "crates/pecos-wasm" } pecos-build-utils = { version = "0.1.1", path = "crates/pecos-build-utils" } pecos-llvm-utils = { version = "0.1.1", path = "crates/pecos-llvm-utils" } +pecos-num = { version = "0.1.1", path = "crates/pecos-num" } # Decoder crates pecos-decoder-core = { version = "0.1.1", path = "crates/pecos-decoder-core" } diff --git a/Makefile b/Makefile index 34b21bb0e..af6f5bedf 100644 --- a/Makefile +++ b/Makefile @@ -327,7 +327,7 @@ decoder-cache-clean: ## Clean decoder download cache .PHONY: pytest pytest: ## Run tests on the Python package (not including optional dependencies). ASSUMES: previous build command @$(ADD_LLVM_TO_PATH) uv run pytest ./python/quantum-pecos/tests/ --doctest-modules -m "not optional_dependency" - @$(ADD_LLVM_TO_PATH) uv run pytest ./python/pecos-rslib/tests/ + @$(ADD_LLVM_TO_PATH) uv run --with scipy --with numpy pytest ./python/pecos-rslib/tests/ @$(ADD_LLVM_TO_PATH) uv run pytest ./python/slr-tests/ -m "not optional_dependency" .PHONY: pytest-dep @@ -337,7 +337,7 @@ pytest-dep: ## Run tests on the Python package only for optional dependencies. A .PHONY: pytest-all pytest-all: ## Run all tests on the Python package ASSUMES: previous build command @$(ADD_LLVM_TO_PATH) uv run pytest ./python/quantum-pecos/tests/ -m "" - @$(ADD_LLVM_TO_PATH) uv run pytest ./python/pecos-rslib/tests/ + @$(ADD_LLVM_TO_PATH) uv run --with scipy --with numpy pytest ./python/pecos-rslib/tests/ # .PHONY: pytest-doc # pydoctest: ## Run doctests with pytest. ASSUMES: A build command was ran previously. ASSUMES: previous build command diff --git a/crates/pecos-num/Cargo.toml b/crates/pecos-num/Cargo.toml new file mode 100644 index 000000000..efaaedb8f --- /dev/null +++ b/crates/pecos-num/Cargo.toml @@ -0,0 +1,32 @@ +[package] +name = "pecos-num" +version.workspace = true +edition.workspace = true +authors.workspace = true +homepage.workspace = true +repository.workspace = true +license.workspace = true +keywords.workspace = true +categories.workspace = true +readme.workspace = true +description = "Numerical computing support for PECOS quantum error correction simulations" + +[dependencies] +# Comprehensive numerical computing library (for Newton, polynomial fitting) +peroxide.workspace = true + +# Root finding (for Brent's method only - Peroxide doesn't have it) +roots.workspace = true + +# Curve fitting (scipy-compatible API, better match than Peroxide's AD-based optimizer) +levenberg-marquardt.workspace = true +nalgebra.workspace = true # Required by levenberg-marquardt + +# Array interface (for API compatibility and return types) +ndarray.workspace = true + +# Logging +log.workspace = true + +[lints] +workspace = true diff --git a/crates/pecos-num/README.md b/crates/pecos-num/README.md new file mode 100644 index 000000000..7f88b6763 --- /dev/null +++ b/crates/pecos-num/README.md @@ -0,0 +1,36 @@ +# pecos-num + +`pecos-num` provides numerical computing support for PECOS quantum error correction simulations. + +This crate brings together numerical computing dependencies and implements functionality needed for QEC analysis in PECOS, including threshold fitting, data analysis, and optimization. It provides APIs with similar functionality to scipy and numpy for numerical operations. + +## Features + +- Root finding algorithms (Brent's method, Newton-Raphson) +- Non-linear curve fitting (Levenberg-Marquardt) +- Polynomial fitting and evaluation +- Built on robust Rust numerical libraries (Peroxide, levenberg-marquardt, nalgebra) + +## Usage + +This is an **internal crate** used by: +- `pecos` - The main PECOS metacrate (via prelude) +- `pecos-rslib` - Python bindings exposing numerical functions + +For direct usage in Rust: + +```rust +use pecos_num::prelude::*; + +// Root finding with Brent's method +let root = brentq(|x| x * x - 2.0, 0.0, 2.0, None).unwrap(); + +// Curve fitting +let result = curve_fit( + |x, params| params[0] * x + params[1], + xdata.view(), + ydata.view(), + p0.view(), + None +).unwrap(); +``` diff --git a/crates/pecos-num/src/curve_fit.rs b/crates/pecos-num/src/curve_fit.rs new file mode 100644 index 000000000..35ce78b46 --- /dev/null +++ b/crates/pecos-num/src/curve_fit.rs @@ -0,0 +1,414 @@ +// Copyright 2024 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +//! Non-linear curve fitting using Levenberg-Marquardt algorithm. +//! +//! This module provides a Rust implementation of `scipy.optimize.curve_fit` +//! using the well-tested `levenberg-marquardt` crate. +//! +//! Note: We use `levenberg-marquardt` instead of Peroxide's optimizer because +//! Peroxide requires AD (automatic differentiation) types, while `scipy.optimize.curve_fit` +//! uses simple float functions. The levenberg-marquardt crate provides a better API match. + +use levenberg_marquardt::{LeastSquaresProblem, LevenbergMarquardt}; +use nalgebra::{DMatrix, DVector, Dyn, Owned}; +use ndarray::{Array1, Array2, ArrayView1}; + +/// Error type for curve fitting operations. +#[derive(Debug, Clone)] +pub enum CurveFitError { + /// Optimization failed to converge + ConvergenceError { message: String }, + /// Invalid input data + InvalidInput { message: String }, + /// Numerical issues during fitting + NumericalIssue { message: String }, +} + +impl std::fmt::Display for CurveFitError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::ConvergenceError { message } => write!(f, "Convergence error: {message}"), + Self::InvalidInput { message } => write!(f, "Invalid input: {message}"), + Self::NumericalIssue { message } => write!(f, "Numerical issue: {message}"), + } + } +} + +impl std::error::Error for CurveFitError {} + +/// Problem struct for Levenberg-Marquardt optimization. +struct CurveFitProblem +where + F: Fn(f64, &[f64]) -> f64, +{ + func: F, + xdata: Vec, + ydata: Vec, + params: DVector, +} + +impl LeastSquaresProblem for CurveFitProblem +where + F: Fn(f64, &[f64]) -> f64, +{ + type ParameterStorage = Owned; + type ResidualStorage = Owned; + type JacobianStorage = Owned; + + fn set_params(&mut self, p: &DVector) { + self.params.copy_from(p); + } + + fn params(&self) -> DVector { + self.params.clone() + } + + fn residuals(&self) -> Option> { + let n = self.xdata.len(); + let mut residuals = DVector::zeros(n); + let param_slice = self.params.as_slice(); + + for (i, (&x, &y)) in self.xdata.iter().zip(self.ydata.iter()).enumerate() { + residuals[i] = (self.func)(x, param_slice) - y; + } + + Some(residuals) + } + + fn jacobian(&self) -> Option> { + let n = self.xdata.len(); + let n_params = self.params.len(); + let mut jacobian = DMatrix::zeros(n, n_params); + + let eps = 1e-8; + let residuals = self.residuals()?; + let param_slice = self.params.as_slice(); + + for j in 0..n_params { + let step = eps * (1.0 + param_slice[j].abs()).max(eps); + + // Create perturbed parameters + let mut params_plus = self.params.clone(); + params_plus[j] += step; + let params_plus_slice = params_plus.as_slice(); + + // Compute residuals with perturbed parameters + for (i, &x) in self.xdata.iter().enumerate() { + let residual_plus = (self.func)(x, params_plus_slice) - self.ydata[i]; + jacobian[(i, j)] = (residual_plus - residuals[i]) / step; + } + } + + Some(jacobian) + } +} + +/// Options for curve fitting. +#[derive(Debug, Clone)] +pub struct CurveFitOptions { + /// Maximum number of iterations + pub maxfev: usize, + /// Tolerance for parameter changes + pub xtol: f64, + /// Tolerance for cost changes + pub ftol: f64, + /// Initial damping parameter (ignored, using crate defaults) + pub lambda: f64, +} + +impl Default for CurveFitOptions { + fn default() -> Self { + Self { + maxfev: 1000, + xtol: 1e-8, + ftol: 1e-8, + lambda: 0.01, + } + } +} + +/// Result from curve fitting. +#[derive(Debug, Clone)] +pub struct CurveFitResult { + /// Optimal parameters + pub params: Array1, + /// Covariance matrix (if available) + pub pcov: Option>, + /// Number of function evaluations + pub nfev: usize, + /// Final cost value + pub cost: f64, +} + +// No problem struct needed - we'll use a closure directly + +/// Fit a non-linear function to data using Levenberg-Marquardt. +/// +/// This is a Rust implementation of `scipy.optimize.curve_fit` using the +/// `levenberg-marquardt` crate for robust, well-tested optimization. +/// +/// # Arguments +/// +/// * `func` - Model function: f(x, params) -> y +/// * `xdata` - Independent variable data +/// * `ydata` - Dependent variable data +/// * `p0` - Initial guess for parameters +/// * `options` - Optional fitting options +/// +/// # Returns +/// +/// Returns the optimal parameters and covariance matrix. +/// +/// # Errors +/// +/// Returns an error if: +/// - Data arrays have different lengths +/// - Optimization fails to converge +/// - Numerical issues during fitting +/// +/// # Examples +/// +/// ``` +/// use pecos_num::curve_fit::{curve_fit, CurveFitOptions}; +/// use ndarray::array; +/// +/// // Fit linear: y = a * x + b +/// fn linear(x: f64, params: &[f64]) -> f64 { +/// params[0] * x + params[1] +/// } +/// +/// let xdata = array![0.0, 1.0, 2.0, 3.0, 4.0]; +/// let ydata = array![1.0, 3.0, 5.0, 7.0, 9.0]; +/// let p0 = array![1.0, 0.0]; +/// +/// let result = curve_fit(linear, xdata.view(), ydata.view(), p0.view(), None).unwrap(); +/// // result.params ≈ [2.0, 1.0] (for y = 2*x + 1) +/// ``` +pub fn curve_fit( + func: F, + xdata: ArrayView1, + ydata: ArrayView1, + p0: ArrayView1, + options: Option, +) -> Result +where + F: Fn(f64, &[f64]) -> f64, +{ + let n = xdata.len(); + + if n != ydata.len() { + return Err(CurveFitError::InvalidInput { + message: format!( + "xdata and ydata must have same length: x={n}, y={}", + ydata.len() + ), + }); + } + + if n < p0.len() { + return Err(CurveFitError::InvalidInput { + message: format!( + "Need at least {} data points for {} parameters, got {n}", + p0.len(), + p0.len() + ), + }); + } + + let opts = options.unwrap_or_default(); + + // Create problem for Levenberg-Marquardt + let problem = CurveFitProblem { + func, + xdata: xdata.to_vec(), + ydata: ydata.to_vec(), + params: DVector::from_vec(p0.to_vec()), + }; + + // Run Levenberg-Marquardt optimization + let (result, report) = LevenbergMarquardt::new() + .with_stepbound(100.0) + .with_patience(opts.maxfev) + .minimize(problem); + + // Check convergence + if !report.termination.was_successful() { + return Err(CurveFitError::ConvergenceError { + message: format!("Optimization did not converge: {:?}", report.termination), + }); + } + + // Get final parameters and residuals + let final_params = result.params(); + let final_residuals = result + .residuals() + .ok_or_else(|| CurveFitError::NumericalIssue { + message: "Failed to compute final residuals".to_string(), + })?; + let cost = final_residuals.dot(&final_residuals); + + // Get Jacobian at solution + let jacobian = result + .jacobian() + .ok_or_else(|| CurveFitError::NumericalIssue { + message: "Failed to compute Jacobian".to_string(), + })?; + + // Compute covariance matrix: (J^T * J)^-1 * variance + let jt_j = jacobian.transpose() * &jacobian; + let pcov = match jt_j.svd(true, true).pseudo_inverse(1e-15) { + Ok(inv) => { + let n_params = final_params.len(); + // Cast to f64 is safe for reasonable dataset sizes (< 2^53 points) + #[allow(clippy::cast_precision_loss)] + let dof = (n as f64 - n_params as f64).max(1.0); + let variance = cost / dof; + let cov_mat = inv * variance; + + // Convert to ndarray + let mut pcov_array = Array2::zeros((n_params, n_params)); + for i in 0..n_params { + for j in 0..n_params { + pcov_array[[i, j]] = cov_mat[(i, j)]; + } + } + Some(pcov_array) + } + Err(_) => None, + }; + + log::debug!( + "curve_fit: converged after {} evaluations with cost={:.6e}", + report.number_of_evaluations, + cost + ); + + Ok(CurveFitResult { + params: Array1::from_vec(final_params.as_slice().to_vec()), + pcov, + nfev: report.number_of_evaluations, + cost, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + #[test] + fn test_curve_fit_linear() { + // Fit y = a*x + b + fn linear(x: f64, params: &[f64]) -> f64 { + params[0] * x + params[1] + } + + let xdata = array![0.0, 1.0, 2.0, 3.0, 4.0]; + let ydata = array![1.0, 3.0, 5.0, 7.0, 9.0]; // y = 2*x + 1 + + let p0 = array![1.0, 0.0]; + + let result = curve_fit(linear, xdata.view(), ydata.view(), p0.view(), None).unwrap(); + + assert!( + (result.params[0] - 2.0).abs() < 1e-6, + "slope should be 2.0, got {}", + result.params[0] + ); + assert!( + (result.params[1] - 1.0).abs() < 1e-6, + "intercept should be 1.0, got {}", + result.params[1] + ); + } + + #[test] + fn test_curve_fit_exponential() { + // Fit y = a * exp(b * x) + fn exponential(x: f64, params: &[f64]) -> f64 { + params[0] * (params[1] * x).exp() + } + + let xdata = array![0.0, 1.0, 2.0, 3.0, 4.0]; + // y = e^x: use std::f64::consts::E for accurate test data + let ydata = array![ + 1.0_f64.exp(), + std::f64::consts::E, + (2.0_f64).exp(), + (3.0_f64).exp(), + (4.0_f64).exp() + ]; + + let p0 = array![1.0, 1.0]; + + let result = curve_fit(exponential, xdata.view(), ydata.view(), p0.view(), None).unwrap(); + + assert!( + (result.params[0] - 1.0).abs() < 0.05, + "coefficient should be ~1.0, got {}", + result.params[0] + ); + assert!( + (result.params[1] - 1.0).abs() < 0.05, + "exponent should be ~1.0, got {}", + result.params[1] + ); + } + + #[test] + fn test_curve_fit_quadratic() { + // Fit y = a*x^2 + b*x + c + fn quadratic(x: f64, params: &[f64]) -> f64 { + params[0] * x * x + params[1] * x + params[2] + } + + let xdata = array![0.0, 1.0, 2.0, 3.0, 4.0]; + let ydata = array![3.0, 6.0, 11.0, 18.0, 27.0]; // y = x^2 + 2*x + 3 + + let p0 = array![1.0, 1.0, 1.0]; + + let result = curve_fit(quadratic, xdata.view(), ydata.view(), p0.view(), None).unwrap(); + + assert!( + (result.params[0] - 1.0).abs() < 1e-6, + "x^2 coef should be 1.0, got {}", + result.params[0] + ); + assert!( + (result.params[1] - 2.0).abs() < 1e-6, + "x coef should be 2.0, got {}", + result.params[1] + ); + assert!( + (result.params[2] - 3.0).abs() < 1e-6, + "constant should be 3.0, got {}", + result.params[2] + ); + } + + #[test] + fn test_curve_fit_insufficient_data() { + fn linear(x: f64, params: &[f64]) -> f64 { + params[0] * x + params[1] + } + + let xdata = array![0.0]; + let ydata = array![1.0]; + let p0 = array![1.0, 0.0]; + + let result = curve_fit(linear, xdata.view(), ydata.view(), p0.view(), None); + assert!(matches!(result, Err(CurveFitError::InvalidInput { .. }))); + } +} diff --git a/crates/pecos-num/src/lib.rs b/crates/pecos-num/src/lib.rs new file mode 100644 index 000000000..bdf374355 --- /dev/null +++ b/crates/pecos-num/src/lib.rs @@ -0,0 +1,37 @@ +// Copyright 2024 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +//! # pecos-num: Numerical Computing for PECOS +//! +//! This crate provides numerical computing functionality for PECOS, serving as a +//! Rust-based replacement for scipy.optimize functions. It offers: +//! +//! - Root finding algorithms (Brent's method, Newton-Raphson) +//! - Curve fitting (Levenberg-Marquardt, polynomial fitting) +//! - Performance improvements over scipy +//! - Better cross-platform support +//! +//! ## Usage +//! +//! This crate is typically accessed through the `pecos::prelude`. Python bindings +//! are provided separately in `pecos-rslib`. + +pub mod curve_fit; +pub mod optimize; +pub mod polynomial; +pub mod prelude; + +pub use curve_fit::{CurveFitError, CurveFitOptions, CurveFitResult, curve_fit}; +pub use optimize::{BrentqOptions, NewtonOptions, OptimizeError, brentq, newton}; +pub use polynomial::{Poly1d, PolynomialError, polyfit}; diff --git a/crates/pecos-num/src/optimize.rs b/crates/pecos-num/src/optimize.rs new file mode 100644 index 000000000..4e27ba21f --- /dev/null +++ b/crates/pecos-num/src/optimize.rs @@ -0,0 +1,374 @@ +// Copyright 2024 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +//! Root finding and optimization algorithms. +//! +//! This module provides implementations of common numerical optimization +//! algorithms, compatible with scipy.optimize API. +//! +//! Uses Peroxide for Newton's method implementation, with scipy-compatible +//! functional wrappers. + +use peroxide::fuga::{NewtonMethod, RootFinder, RootFindingProblem, anyhow}; +use std::fmt; + +/// Error type for optimization functions. +#[derive(Debug, Clone)] +pub enum OptimizeError { + /// Function values at interval endpoints have the same sign + SameSigns { fa: f64, fb: f64 }, + /// Maximum iterations exceeded without convergence + MaxIterations { iterations: usize }, + /// Derivative is zero or near-zero + ZeroDerivative { x: f64, derivative: f64 }, + /// Numerical issues (NaN, Inf encountered) + NumericalIssue { message: String }, + /// Convergence criterion not met + ConvergenceFailed { message: String }, +} + +impl fmt::Display for OptimizeError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::SameSigns { fa, fb } => { + write!( + f, + "f(a) and f(b) must have opposite signs. Got f(a)={fa}, f(b)={fb}" + ) + } + Self::MaxIterations { iterations } => { + write!(f, "Maximum iterations ({iterations}) exceeded") + } + Self::ZeroDerivative { x, derivative } => { + write!(f, "Derivative is zero at x={x} (derivative={derivative})") + } + Self::NumericalIssue { message } => { + write!(f, "Numerical issue: {message}") + } + Self::ConvergenceFailed { message } => { + write!(f, "Convergence failed: {message}") + } + } + } +} + +impl std::error::Error for OptimizeError {} + +/// Options for Brent's method root finding. +#[derive(Debug, Clone)] +pub struct BrentqOptions { + /// Absolute tolerance for root finding + pub xtol: f64, + /// Relative tolerance for root finding + pub rtol: f64, + /// Maximum number of iterations + pub maxiter: usize, +} + +impl Default for BrentqOptions { + fn default() -> Self { + Self { + xtol: 2e-12, + rtol: 8.881_784_197_001_252e-16, // scipy default + maxiter: 100, + } + } +} + +/// Options for Newton-Raphson method. +#[derive(Debug, Clone)] +pub struct NewtonOptions { + /// Absolute tolerance for convergence + pub tol: f64, + /// Maximum number of iterations + pub maxiter: usize, + /// Step size for numerical derivative (if fprime not provided) + pub eps: f64, +} + +impl Default for NewtonOptions { + fn default() -> Self { + Self { + tol: 1.48e-8, // scipy default + maxiter: 50, + eps: 1e-8, + } + } +} + +/// Find root of a function using Brent's method. +/// +/// This is a Rust implementation of scipy.optimize.brentq. +/// +/// Brent's method combines root bracketing, bisection, and inverse quadratic +/// interpolation. It is generally considered one of the best methods for +/// finding roots of a continuous function. +/// +/// # Arguments +/// +/// * `f` - Function for which to find root +/// * `a` - Lower bound of interval +/// * `b` - Upper bound of interval +/// * `options` - Optional configuration parameters +/// +/// # Returns +/// +/// Returns the root of the function within the interval [a, b]. +/// +/// # Errors +/// +/// Returns an error if: +/// - f(a) and f(b) have the same sign +/// - Maximum iterations exceeded +/// - Numerical issues encountered +/// +/// # Examples +/// +/// ``` +/// use pecos_num::optimize::{brentq, BrentqOptions}; +/// +/// // Find root of f(x) = x^2 - 2 (should be sqrt(2)) +/// let root = brentq(|x| x * x - 2.0, 0.0, 2.0, None).unwrap(); +/// assert!((root - 2f64.sqrt()).abs() < 1e-10); +/// ``` +pub fn brentq(f: F, a: f64, b: f64, options: Option) -> Result +where + F: Fn(f64) -> f64, +{ + let opts = options.unwrap_or_default(); + + // Use roots crate for Brent's method + let mut convergency = roots::SimpleConvergency { + eps: opts.xtol, + max_iter: opts.maxiter, + }; + + let result = roots::find_root_brent(a, b, &f, &mut convergency); + + match result { + Ok(root) => { + log::debug!("brentq converged to root={root}"); + Ok(root) + } + Err(e) => { + log::warn!("brentq failed: {e:?}"); + // Check if it's a sign issue + let fa = f(a); + let fb = f(b); + if fa * fb > 0.0 { + Err(OptimizeError::SameSigns { fa, fb }) + } else { + Err(OptimizeError::MaxIterations { + iterations: opts.maxiter, + }) + } + } + } +} + +/// Internal wrapper for Newton's method using Peroxide. +struct NewtonProblem +where + F: Fn(f64) -> f64, + G: Fn(f64) -> f64, +{ + f: F, + fprime: Option, + eps: f64, + x0: f64, +} + +impl RootFindingProblem<1, 1, f64> for NewtonProblem +where + F: Fn(f64) -> f64, + G: Fn(f64) -> f64, +{ + fn function(&self, x: [f64; 1]) -> Result<[f64; 1], anyhow::Error> { + Ok([(self.f)(x[0])]) + } + + fn derivative(&self, x: [f64; 1]) -> Result<[[f64; 1]; 1], anyhow::Error> { + let fprime_x = if let Some(ref fprime_fn) = self.fprime { + (fprime_fn)(x[0]) + } else { + // Numerical derivative using finite differences + let h = self.eps; + let fx = (self.f)(x[0]); + let fx_plus_h = (self.f)(x[0] + h); + (fx_plus_h - fx) / h + }; + + Ok([[fprime_x]]) + } + + fn initial_guess(&self) -> f64 { + self.x0 + } +} + +/// Find root using Newton-Raphson method. +/// +/// This is a scipy.optimize.newton-compatible wrapper around Peroxide's Newton implementation. +/// +/// Newton's method uses the function value and its derivative to iteratively +/// converge to a root. It typically converges quickly when close to the root, +/// but may fail if the initial guess is poor or the derivative is zero. +/// +/// # Arguments +/// +/// * `f` - Function for which to find root +/// * `x0` - Initial guess +/// * `fprime` - Optional derivative function. If None, uses numerical derivative. +/// * `options` - Optional configuration parameters +/// +/// # Returns +/// +/// Returns the root of the function. +/// +/// # Errors +/// +/// Returns an error if: +/// - Maximum iterations exceeded +/// - Derivative is zero or near-zero +/// - Numerical issues encountered +/// +/// # Examples +/// +/// ``` +/// use pecos_num::optimize::{newton, NewtonOptions}; +/// +/// // Find root of f(x) = x^2 - 2 (should be sqrt(2)) +/// // With derivative f'(x) = 2x +/// let root = newton( +/// |x| x * x - 2.0, +/// 1.0, // initial guess +/// Some(|x| 2.0 * x), // derivative +/// None +/// ).unwrap(); +/// assert!((root - 2f64.sqrt()).abs() < 1e-10); +/// ``` +pub fn newton( + f: F, + x0: f64, + fprime: Option, + options: Option, +) -> Result +where + F: Fn(f64) -> f64, + G: Fn(f64) -> f64, +{ + let opts = options.unwrap_or_default(); + + log::debug!("newton starting from x0={x0}"); + + // Create Peroxide problem + let problem = NewtonProblem { + f, + fprime, + eps: opts.eps, + x0, + }; + + // Create Peroxide Newton method + let method = NewtonMethod { + max_iter: opts.maxiter, + tol: opts.tol, + }; + + // Solve using Peroxide + let result = method.find(&problem); + + match result { + Ok(root) => { + log::debug!("newton converged to root={}", root[0]); + Ok(root[0]) + } + Err(e) => { + log::warn!("newton failed: {e:?}"); + Err(OptimizeError::ConvergenceFailed { + message: format!("{e:?}"), + }) + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_brentq_sqrt2() { + // Find sqrt(2) by solving x^2 - 2 = 0 + let root = brentq(|x| x * x - 2.0, 0.0, 2.0, None).unwrap(); + assert!((root - 2f64.sqrt()).abs() < 1e-10); + } + + #[test] + fn test_brentq_cubic() { + // Find root of x^3 - x - 2 = 0 (root is approximately 1.52138) + let root = brentq(|x| x.powi(3) - x - 2.0, 1.0, 2.0, None).unwrap(); + let expected = 1.521_379_706_804_567_6; + assert!((root - expected).abs() < 1e-10); + } + + #[test] + fn test_brentq_same_signs() { + // Should fail when f(a) and f(b) have same sign + let result = brentq(|x| x * x + 1.0, -1.0, 1.0, None); + assert!(matches!(result, Err(OptimizeError::SameSigns { .. }))); + } + + #[test] + fn test_newton_sqrt2() { + // Find sqrt(2) using Newton's method with derivative + let root = newton(|x| x * x - 2.0, 1.0, Some(|x: f64| 2.0 * x), None).unwrap(); + assert!((root - 2f64.sqrt()).abs() < 1e-10); + } + + #[test] + fn test_newton_numerical_derivative() { + // Find sqrt(2) using Newton's method with numerical derivative + let root = newton(|x| x * x - 2.0, 1.0, None:: f64>, None).unwrap(); + assert!((root - 2f64.sqrt()).abs() < 1e-8); + } + + #[test] + fn test_newton_cubic() { + // Find root of x^3 - x - 2 = 0 + let root = newton( + |x| x.powi(3) - x - 2.0, + 1.5, + Some(|x: f64| 3.0 * x.powi(2) - 1.0), + None, + ) + .unwrap(); + let expected = 1.521_379_706_804_567_6; + assert!((root - expected).abs() < 1e-10); + } + + #[test] + fn test_newton_polynomial_root() { + // Find root of (x-3)(x-5) = x^2 - 8x + 15 = 0 + // Should find root near initial guess (close to 5) + let root = newton( + |x| x * x - 8.0 * x + 15.0, + 4.5, // Start at 4.5, not 4.0 (which has zero derivative) + Some(|x: f64| 2.0 * x - 8.0), + None, + ) + .unwrap(); + // Should converge to 5 since we start at 4.5 + assert!((root - 5.0).abs() < 1e-10); + } +} diff --git a/crates/pecos-num/src/polynomial.rs b/crates/pecos-num/src/polynomial.rs new file mode 100644 index 000000000..74148f239 --- /dev/null +++ b/crates/pecos-num/src/polynomial.rs @@ -0,0 +1,299 @@ +// Copyright 2024 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +//! Polynomial fitting and evaluation. +//! +//! This module provides implementations of polynomial operations, +//! compatible with numpy.polyfit and numpy.poly1d API. +//! +//! Uses Peroxide for linear algebra (SVD solving). + +use ndarray::{Array1, ArrayView1}; +use peroxide::fuga::{Col, LU, LinearAlgebra, MatrixTrait, Row, matrix}; + +/// Error type for polynomial operations. +#[derive(Debug, Clone)] +pub enum PolynomialError { + /// Insufficient data points for the requested degree + InsufficientData { num_points: usize, degree: usize }, + /// Numerical issues during fitting + NumericalIssue { message: String }, + /// Linear algebra error + LinAlgError { message: String }, +} + +impl std::fmt::Display for PolynomialError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::InsufficientData { num_points, degree } => { + write!( + f, + "Insufficient data: need at least {} points for degree {}, got {}", + degree + 1, + degree, + num_points + ) + } + Self::NumericalIssue { message } => write!(f, "Numerical issue: {message}"), + Self::LinAlgError { message } => write!(f, "Linear algebra error: {message}"), + } + } +} + +impl std::error::Error for PolynomialError {} + +/// Fit a polynomial of given degree to data points. +/// +/// This is a Rust implementation of numpy.polyfit. +/// +/// # Arguments +/// +/// * `x` - x-coordinates of data points +/// * `y` - y-coordinates of data points +/// * `deg` - Degree of the polynomial fit +/// +/// # Returns +/// +/// Returns the polynomial coefficients in decreasing order of degree. +/// For example, for degree 2: [c0, c1, c2] where y = c0*x^2 + c1*x + c2 +/// +/// # Errors +/// +/// Returns an error if: +/// - Not enough data points for the requested degree +/// - Numerical issues during fitting +/// +/// # Examples +/// +/// ``` +/// use pecos_num::polynomial::polyfit; +/// use ndarray::array; +/// +/// // Fit y = 2x + 1 +/// let x = array![0.0, 1.0, 2.0, 3.0]; +/// let y = array![1.0, 3.0, 5.0, 7.0]; +/// let coeffs = polyfit(x.view(), y.view(), 1).unwrap(); +/// assert!((coeffs[0] - 2.0).abs() < 1e-10); // slope +/// assert!((coeffs[1] - 1.0).abs() < 1e-10); // intercept +/// ``` +pub fn polyfit( + x: ArrayView1, + y: ArrayView1, + deg: usize, +) -> Result, PolynomialError> { + let n = x.len(); + + if n != y.len() { + return Err(PolynomialError::NumericalIssue { + message: format!("x and y must have same length: x={n}, y={}", y.len()), + }); + } + + if n < deg + 1 { + return Err(PolynomialError::InsufficientData { + num_points: n, + degree: deg, + }); + } + + // Build Vandermonde matrix using Peroxide + // For degree 2: [[x0^2, x0, 1], [x1^2, x1, 1], ...] + // Flatten to 1D vec for Peroxide's matrix constructor + let mut vandermonde_data = Vec::with_capacity(n * (deg + 1)); + for &xi in x { + for j in 0..=deg { + // Cast is safe: polynomial degrees are always << i32::MAX + #[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)] + let power = (deg - j) as i32; + vandermonde_data.push(xi.powi(power)); + } + } + let vandermonde = matrix(vandermonde_data, n, deg + 1, Row); + + // Convert y to vector and then to column matrix + let y_vec: Vec = y.iter().copied().collect(); + let y_mat = matrix(y_vec.clone(), n, 1, Col); + + // Solve least squares: coeffs = (A^T A)^{-1} A^T y + // where A is the Vandermonde matrix + let at = vandermonde.t(); // A^T + let gram_matrix = &at * &vandermonde; // A^T A (Gram matrix) + let at_y = &at * &y_mat; // A^T y + + // Solve the normal equations + let at_y_vec: Vec = at_y.data.clone(); + let coeffs_vec = gram_matrix.solve(&at_y_vec, LU); + + // Convert back to ndarray + let coeffs = Array1::from_vec(coeffs_vec); + + log::debug!("polyfit: fitted polynomial of degree {deg} with coeffs: {coeffs:?}"); + + Ok(coeffs) +} + +/// Polynomial class for evaluation. +/// +/// This is a Rust implementation of numpy.poly1d functionality. +#[derive(Debug, Clone)] +pub struct Poly1d { + /// Polynomial coefficients in decreasing order of degree + /// For [c0, c1, c2]: y = c0*x^2 + c1*x + c2 + coeffs: Array1, +} + +impl Poly1d { + /// Create a new polynomial from coefficients. + /// + /// # Arguments + /// + /// * `coeffs` - Coefficients in decreasing order of degree + /// + /// # Examples + /// + /// ``` + /// use pecos_num::polynomial::Poly1d; + /// use ndarray::array; + /// + /// // Create polynomial: 2x^2 + 3x + 1 + /// let p = Poly1d::new(array![2.0, 3.0, 1.0]); + /// assert_eq!(p.eval(0.0), 1.0); // p(0) = 1 + /// assert_eq!(p.eval(1.0), 6.0); // p(1) = 2 + 3 + 1 = 6 + /// ``` + #[must_use] + pub fn new(coeffs: Array1) -> Self { + Self { coeffs } + } + + /// Evaluate the polynomial at a given value. + /// + /// Uses Horner's method for efficient evaluation. + /// + /// # Arguments + /// + /// * `x` - Value at which to evaluate the polynomial + /// + /// # Returns + /// + /// The value of the polynomial at x + /// + /// # Panics + /// + /// Panics if the coefficient array is not in standard layout (contiguous in memory). + #[must_use] + pub fn eval(&self, x: f64) -> f64 { + if self.coeffs.is_empty() { + return 0.0; + } + + // Horner's method: a0 + x(a1 + x(a2 + x(...))) + let mut result = self.coeffs[0]; + for &coeff in &self.coeffs.as_slice().unwrap()[1..] { + result = result * x + coeff; + } + result + } + + /// Get the degree of the polynomial. + #[must_use] + pub fn degree(&self) -> usize { + if self.coeffs.is_empty() { + 0 + } else { + self.coeffs.len() - 1 + } + } + + /// Get the coefficients. + #[must_use] + pub fn coefficients(&self) -> &Array1 { + &self.coeffs + } +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + #[test] + fn test_polyfit_linear() { + // Fit y = 2x + 1 + let x = array![0.0, 1.0, 2.0, 3.0, 4.0]; + let y = array![1.0, 3.0, 5.0, 7.0, 9.0]; + + let coeffs = polyfit(x.view(), y.view(), 1).unwrap(); + + assert_eq!(coeffs.len(), 2); + assert!((coeffs[0] - 2.0).abs() < 1e-10); // slope + assert!((coeffs[1] - 1.0).abs() < 1e-10); // intercept + } + + #[test] + fn test_polyfit_quadratic() { + // Fit y = x^2 + 2x + 3 + let x = array![0.0, 1.0, 2.0, 3.0, 4.0]; + let y = array![3.0, 6.0, 11.0, 18.0, 27.0]; + + let coeffs = polyfit(x.view(), y.view(), 2).unwrap(); + + assert_eq!(coeffs.len(), 3); + assert!((coeffs[0] - 1.0).abs() < 1e-10); // x^2 + assert!((coeffs[1] - 2.0).abs() < 1e-10); // x + assert!((coeffs[2] - 3.0).abs() < 1e-10); // constant + } + + #[test] + fn test_poly1d_eval() { + // Test polynomial: 2x^2 + 3x + 1 + let p = Poly1d::new(array![2.0, 3.0, 1.0]); + + // Allow exact float comparison for simple polynomial evaluations with integer coefficients + #[allow(clippy::float_cmp)] + { + assert_eq!(p.eval(0.0), 1.0); // p(0) = 1 + assert_eq!(p.eval(1.0), 6.0); // p(1) = 2 + 3 + 1 = 6 + assert_eq!(p.eval(2.0), 15.0); // p(2) = 8 + 6 + 1 = 15 + assert_eq!(p.eval(-1.0), 0.0); // p(-1) = 2 - 3 + 1 = 0 + } + } + + #[test] + fn test_polyfit_and_eval() { + // Fit a polynomial and check evaluation + let x = array![0.0, 1.0, 2.0, 3.0, 4.0]; + let y = array![1.0, 3.0, 5.0, 7.0, 9.0]; + + let coeffs = polyfit(x.view(), y.view(), 1).unwrap(); + let p = Poly1d::new(coeffs); + + // Check that polynomial evaluates correctly at training points + for (xi, yi) in x.iter().zip(y.iter()) { + assert!((p.eval(*xi) - yi).abs() < 1e-10); + } + } + + #[test] + fn test_polyfit_insufficient_data() { + let x = array![0.0, 1.0]; + let y = array![1.0, 2.0]; + + // Try to fit degree 3 polynomial with only 2 points + let result = polyfit(x.view(), y.view(), 3); + assert!(matches!( + result, + Err(PolynomialError::InsufficientData { .. }) + )); + } +} diff --git a/crates/pecos-num/src/prelude.rs b/crates/pecos-num/src/prelude.rs new file mode 100644 index 000000000..1f253e737 --- /dev/null +++ b/crates/pecos-num/src/prelude.rs @@ -0,0 +1,24 @@ +// Copyright 2025 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except +// in compliance with the License.You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software distributed under the License +// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express +// or implied. See the License for the specific language governing permissions and limitations under +// the License. + +//! A prelude for users of the `pecos-num` crate. +//! +//! This prelude re-exports numerical computing functions that replace scipy.optimize. + +// Re-export curve fitting +pub use crate::curve_fit::{CurveFitError, CurveFitOptions, CurveFitResult, curve_fit}; + +// Re-export optimization algorithms +pub use crate::optimize::{BrentqOptions, NewtonOptions, OptimizeError, brentq, newton}; + +// Re-export polynomial fitting +pub use crate::polynomial::{Poly1d, PolynomialError, polyfit}; diff --git a/crates/pecos/Cargo.toml b/crates/pecos/Cargo.toml index d6268b0ce..d1a75367d 100644 --- a/crates/pecos/Cargo.toml +++ b/crates/pecos/Cargo.toml @@ -29,6 +29,7 @@ pecos-llvm = { workspace = true, optional = true } pecos-hugr-qis = { workspace = true, optional = true } pecos-phir = { workspace = true, features = ["hugr"] } pecos-rng.workspace = true +pecos-num.workspace = true log.workspace = true tempfile.workspace = true serde_json.workspace = true diff --git a/crates/pecos/src/prelude.rs b/crates/pecos/src/prelude.rs index aefb9a330..c5d5b3026 100644 --- a/crates/pecos/src/prelude.rs +++ b/crates/pecos/src/prelude.rs @@ -48,6 +48,7 @@ //! - `pecos_qis_selene::prelude` - Selene-based QIS interface (when `selene` feature enabled) //! - `pecos_programs::prelude` - Program type definitions //! - `pecos_rng::prelude` - Random number generation +//! - `pecos_num::prelude` - Numerical computing (scipy.optimize replacement) //! - `pecos_hugr_qis::prelude` - HUGR to QIS compilation //! - `pecos_phir_json::prelude` - PHIR-JSON format support //! @@ -88,6 +89,9 @@ pub use pecos_programs::prelude::*; // Re-export RNG prelude pub use pecos_rng::prelude::*; +// Re-export numerical computing prelude +pub use pecos_num::prelude::*; + // Re-export HUGR compiler prelude #[cfg(feature = "llvm")] pub use pecos_hugr_qis::prelude::*; diff --git a/python/pecos-rslib/pyproject.toml b/python/pecos-rslib/pyproject.toml index 50168efed..c2ad44386 100644 --- a/python/pecos-rslib/pyproject.toml +++ b/python/pecos-rslib/pyproject.toml @@ -43,6 +43,11 @@ manifest-path = "rust/Cargo.toml" dev = [ "patchelf; platform_system != 'Windows'", # For setting rpath in shared libraries during development (Linux/macOS only) ] +test = [ + "pytest>=7.0", + "numpy>=1.20", + "scipy>=1.7", # For comparison tests only +] [tool.uv.sources] pecos-rslib = { workspace = true } diff --git a/python/pecos-rslib/rust/Cargo.toml b/python/pecos-rslib/rust/Cargo.toml index 755252fab..c7018d825 100644 --- a/python/pecos-rslib/rust/Cargo.toml +++ b/python/pecos-rslib/rust/Cargo.toml @@ -25,10 +25,12 @@ wasm = ["pecos/wasm"] [dependencies] # Use the pecos metacrate which includes all simulators and runtimes by default -# Enable llvm and wasm features explicitly for full Python functionality +# Enable llvm and wasm features for full Python functionality +# pecos-num is included by default in pecos pecos = { workspace = true, features = ["llvm", "wasm"] } pyo3 = { workspace=true, features = ["extension-module", "abi3-py310", "generate-import-lib"] } +numpy = "0.27" parking_lot.workspace = true regex.workspace = true serde_json.workspace = true diff --git a/python/pecos-rslib/rust/src/lib.rs b/python/pecos-rslib/rust/src/lib.rs index 47011cbd7..d4c086bdc 100644 --- a/python/pecos-rslib/rust/src/lib.rs +++ b/python/pecos-rslib/rust/src/lib.rs @@ -22,6 +22,7 @@ mod cpp_sparse_sim_bindings; mod engine_bindings; mod engine_builders; mod noise_helpers; +mod num_bindings; mod pauli_prop_bindings; // mod pcg_bindings; mod hugr_compilation_bindings; @@ -150,6 +151,9 @@ fn _pecos_rslib(_py: Python<'_>, m: &Bound<'_, PyModule>) -> PyResult<()> { // Register binding module for LLVM bitcode generation llvm_bindings::register_binding_module(m)?; + // Register numerical computing module (scipy.optimize replacements) + num_bindings::register_num_module(m)?; + // Register program types m.add_class::()?; m.add_class::()?; diff --git a/python/pecos-rslib/rust/src/num_bindings.rs b/python/pecos-rslib/rust/src/num_bindings.rs new file mode 100644 index 000000000..97afa7b22 --- /dev/null +++ b/python/pecos-rslib/rust/src/num_bindings.rs @@ -0,0 +1,643 @@ +// Copyright 2024 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +//! Python bindings for pecos-num numerical computing functions. +//! +//! This module provides drop-in replacements for scipy.optimize functions, +//! implemented in Rust for better performance and easier deployment. + +use numpy::ndarray::Array1; +use numpy::{PyArray1, PyArray2, PyReadonlyArray1}; +use pyo3::prelude::*; +use pyo3::types::PyTuple; + +// Import numerical computing types from pecos prelude +// Functions are accessed via pecos::prelude module +use pecos::prelude::{ + BrentqOptions, CurveFitError, CurveFitOptions, NewtonOptions, Poly1d as RustPoly1d, +}; + +/// Helper function to convert `CurveFitError` to appropriate Python exception. +/// +/// Maps Rust errors to Python exceptions following `scipy.optimize.curve_fit` conventions: +/// - `ConvergenceError` -> `RuntimeError` (scipy raises `RuntimeError` for convergence failures) +/// - `InvalidInput` -> `ValueError` (standard Python convention for invalid inputs) +/// - `NumericalIssue` -> `RuntimeError` (similar to convergence issues) +fn map_curve_fit_error(error: CurveFitError) -> PyErr { + match error { + CurveFitError::InvalidInput { message } => { + PyErr::new::(format!("curve_fit failed: {message}")) + } + CurveFitError::ConvergenceError { message } | CurveFitError::NumericalIssue { message } => { + PyErr::new::(format!( + "curve_fit failed: {message}" + )) + } + } +} + +/// Find root of a function using Brent's method. +/// +/// This is a drop-in replacement for scipy.optimize.brentq. +/// +/// Args: +/// f: Callable[[float], float] - Function for which to find root +/// a: float - Lower bound of interval +/// b: float - Upper bound of interval +/// xtol: float - Absolute tolerance (default: 2e-12) +/// rtol: float - Relative tolerance (default: 8.881784197001252e-16) +/// maxiter: int - Maximum iterations (default: 100) +/// +/// Returns: +/// float: The root of the function +/// +/// Raises: +/// `ValueError`: If f(a) and f(b) have the same sign +/// `RuntimeError`: If maximum iterations exceeded +/// +/// Examples: +/// >>> from `pecos_rslib.num` import brentq +/// >>> # Find sqrt(2) by solving x^2 - 2 = 0 +/// >>> root = brentq(lambda x: x**2 - 2, 0, 2) +/// >>> abs(root - 2**0.5) < 1e-10 +/// True +#[pyfunction] +#[pyo3(signature = (f, a, b, xtol=None, rtol=None, maxiter=None))] +#[allow(clippy::needless_pass_by_value)] // Py is a cheap ref-counted pointer; closure needs ownership +fn brentq( + _py: Python<'_>, + f: Py, + a: f64, + b: f64, + xtol: Option, + rtol: Option, + maxiter: Option, +) -> PyResult { + // Create closure that calls Python function + let func = |x: f64| -> f64 { + Python::attach(|py| { + f.call1(py, (x,)) + .and_then(|result| result.extract::(py)) + .unwrap_or(f64::NAN) + }) + }; + + // Configure options + let opts = BrentqOptions { + xtol: xtol.unwrap_or(2e-12), + rtol: rtol.unwrap_or(8.881_784_197_001_252e-16), + maxiter: maxiter.unwrap_or(100), + }; + + // Call Rust implementation + pecos::prelude::brentq(func, a, b, Some(opts)) + .map_err(|e| PyErr::new::(format!("brentq failed: {e}"))) +} + +/// Find root using Newton-Raphson method. +/// +/// This is a drop-in replacement for scipy.optimize.newton. +/// +/// Args: +/// func: Callable[[float], float] - Function for which to find root +/// x0: float - Initial guess +/// fprime: Optional[Callable[[float], float]] - Derivative function (default: None uses numerical derivative) +/// tol: float - Convergence tolerance (default: 1.48e-8) +/// maxiter: int - Maximum iterations (default: 50) +/// +/// Returns: +/// float: The root of the function +/// +/// Raises: +/// `ValueError`: If derivative is zero +/// `RuntimeError`: If maximum iterations exceeded or convergence fails +/// +/// Examples: +/// >>> from `pecos_rslib.num` import newton +/// >>> # Find sqrt(2) by solving x^2 - 2 = 0 +/// >>> root = newton(lambda x: x**2 - 2, x0=1.0, fprime=lambda x: 2*x) +/// >>> abs(root - 2**0.5) < 1e-10 +/// True +#[pyfunction] +#[pyo3(signature = (func, x0, fprime=None, tol=None, maxiter=None))] +#[allow(clippy::needless_pass_by_value)] // Py is a cheap ref-counted pointer; closures need ownership +fn newton( + _py: Python<'_>, + func: Py, + x0: f64, + fprime: Option>, + tol: Option, + maxiter: Option, +) -> PyResult { + // Create closure for function + let f = |x: f64| -> f64 { + Python::attach(|py| { + func.call1(py, (x,)) + .and_then(|result| result.extract::(py)) + .unwrap_or(f64::NAN) + }) + }; + + // Configure options + let opts = NewtonOptions { + tol: tol.unwrap_or(1.48e-8), + maxiter: maxiter.unwrap_or(50), + eps: 1e-8, + }; + + // Call Rust implementation + let result = if let Some(fprime_fn) = fprime { + // Use provided derivative + let fprime_closure = |x: f64| -> f64 { + Python::attach(|py| { + fprime_fn + .call1(py, (x,)) + .and_then(|result| result.extract::(py)) + .unwrap_or(f64::NAN) + }) + }; + pecos::prelude::newton(f, x0, Some(fprime_closure), Some(opts)) + } else { + // Use numerical derivative + pecos::prelude::newton(f, x0, None:: f64>, Some(opts)) + }; + + result.map_err(|e| { + PyErr::new::(format!("newton failed: {e}")) + }) +} + +/// Fit a polynomial of given degree to data points. +/// +/// This is a drop-in replacement for numpy.polyfit. +/// +/// Args: +/// x: `array_like` - x-coordinates of data points +/// y: `array_like` - y-coordinates of data points +/// deg: int - Degree of the polynomial fit +/// +/// Returns: +/// ndarray: Polynomial coefficients in decreasing order of degree +/// For example, for degree 2: [c0, c1, c2] where y = c0*x^2 + c1*x + c2 +/// +/// Raises: +/// `ValueError`: If not enough data points for the requested degree +/// `RuntimeError`: If numerical issues during fitting +/// +/// Examples: +/// >>> from `pecos_rslib.num` import polyfit +/// >>> import numpy as np +/// >>> # Fit y = 2x + 1 +/// >>> x = np.array([0.0, 1.0, 2.0, 3.0]) +/// >>> y = np.array([1.0, 3.0, 5.0, 7.0]) +/// >>> coeffs = polyfit(x, y, 1) +/// >>> # coeffs ≈ [2.0, 1.0] (slope, intercept) +#[pyfunction] +#[allow(clippy::needless_pass_by_value)] // PyReadonlyArray1 is a lightweight wrapper +fn polyfit( + py: Python<'_>, + x: PyReadonlyArray1, + y: PyReadonlyArray1, + deg: usize, +) -> PyResult>> { + let x_view = x.as_array(); + let y_view = y.as_array(); + + let coeffs = pecos::prelude::polyfit(x_view, y_view, deg).map_err(|e| { + PyErr::new::(format!("polyfit failed: {e}")) + })?; + + Ok(PyArray1::from_array(py, &coeffs).unbind()) +} + +/// Polynomial class for evaluation. +/// +/// This is a drop-in replacement for numpy.poly1d. +/// +/// Examples: +/// >>> from `pecos_rslib.num` import Poly1d +/// >>> import numpy as np +/// >>> # Create polynomial: 2x^2 + 3x + 1 +/// >>> p = Poly1d(np.array([2.0, 3.0, 1.0])) +/// >>> p.eval(0.0) # p(0) = 1 +/// 1.0 +/// >>> p.eval(1.0) # p(1) = 2 + 3 + 1 = 6 +/// 6.0 +#[pyclass] +struct Poly1d { + inner: RustPoly1d, +} + +#[pymethods] +impl Poly1d { + /// Create a new polynomial from coefficients. + /// + /// Args: + /// coeffs: `array_like` - Coefficients in decreasing order of degree + #[new] + #[allow(clippy::needless_pass_by_value)] // PyReadonlyArray1 is a lightweight wrapper + fn new(coeffs: PyReadonlyArray1) -> Self { + let coeffs_array = coeffs.as_array().to_owned(); + Self { + inner: RustPoly1d::new(coeffs_array), + } + } + + /// Evaluate the polynomial at a given value. + /// + /// Args: + /// x: float - Value at which to evaluate the polynomial + /// + /// Returns: + /// float: The value of the polynomial at x + fn eval(&self, x: f64) -> f64 { + self.inner.eval(x) + } + + /// Get the degree of the polynomial. + /// + /// Returns: + /// int: Degree of the polynomial + fn degree(&self) -> usize { + self.inner.degree() + } + + /// Get the polynomial coefficients. + /// + /// Returns: + /// ndarray: Coefficients in decreasing order of degree + fn coefficients(&self, py: Python<'_>) -> Py> { + PyArray1::from_array(py, self.inner.coefficients()).unbind() + } + + /// Call the polynomial (same as eval). + fn __call__(&self, x: f64) -> f64 { + self.inner.eval(x) + } + + /// String representation of the polynomial. + fn __repr__(&self) -> String { + format!("Poly1d(coefficients={:?})", self.inner.coefficients()) + } +} + +/// Fit a non-linear function to data using Levenberg-Marquardt. +/// +/// This is a drop-in replacement for `scipy.optimize.curve_fit`. +/// +/// Args: +/// f: Callable[[float, array], float] - Model function f(x, params) or f((x1, x2, ...), params) +/// xdata: `array_like` or tuple of arrays - Independent variable data (can be single array or tuple of arrays) +/// ydata: `array_like` - Dependent variable data +/// p0: `array_like` - Initial guess for parameters +/// maxfev: int - Maximum function evaluations (default: 1000) +/// xtol: float - Parameter tolerance (default: 1e-8) +/// ftol: float - Cost tolerance (default: 1e-8) +/// +/// Returns: +/// tuple: (popt, pcov) - Optimal parameters and covariance matrix +/// +/// Raises: +/// `ValueError`: If data arrays have different lengths +/// `RuntimeError`: If optimization fails to converge +/// +/// Examples: +/// >>> from `pecos_rslib.num` import `curve_fit` +/// >>> import numpy as np +/// >>> # Example 1: Single independent variable +/// >>> def func(x, a, b): +/// ... return a * x + b +/// >>> xdata = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) +/// >>> ydata = np.array([1.0, 3.0, 5.0, 7.0, 9.0]) +/// >>> p0 = np.array([1.0, 0.0]) +/// >>> popt, pcov = `curve_fit(func`, xdata, ydata, p0) +/// >>> # popt ≈ [2.0, 1.0] +/// >>> +/// >>> # Example 2: Multiple independent variables (tuple of arrays) +/// >>> def func2(x, a, b): +/// ... p, d = x # Unpack tuple +/// ... return a * p ** (b / d) +/// >>> pdata = np.array([0.1, 0.2, 0.3]) +/// >>> ddata = np.array([3.0, 3.0, 3.0]) +/// >>> ydata2 = np.array([0.5, 0.7, 0.9]) +/// >>> popt2, pcov2 = `curve_fit(func2`, (pdata, ddata), ydata2, np.array([1.0, 1.0])) +#[pyfunction] +#[pyo3(signature = (f, xdata, ydata, p0, maxfev=None, xtol=None, ftol=None))] +#[allow(clippy::type_complexity)] // Complex return type required for scipy compatibility +#[allow(clippy::too_many_arguments)] // scipy.optimize.curve_fit has many parameters +#[allow(clippy::needless_pass_by_value)] // PyReadonlyArray1 is a lightweight wrapper +fn curve_fit<'py>( + py: Python<'py>, + f: Py, + xdata: &Bound<'py, PyAny>, + ydata: PyReadonlyArray1, + p0: &Bound<'py, PyAny>, + maxfev: Option, + xtol: Option, + ftol: Option, +) -> PyResult<(Py>, Py>)> { + // Convert p0 to array (accept array, tuple, or list) + let p0_array = if let Ok(array) = p0.extract::>() { + array + } else if let Ok(tuple) = p0.cast() { + // Convert tuple to array + let values: Vec = tuple.extract()?; + let np = py.import("numpy")?; + let array = np.call_method1("array", (values,))?; + array.extract::>()? + } else if let Ok(list) = p0.extract::>() { + // Convert list to array + let np = py.import("numpy")?; + let array = np.call_method1("array", (list,))?; + array.extract::>()? + } else { + return Err(PyErr::new::( + "p0 must be an array, tuple, or list", + )); + }; + + // Check if xdata is a tuple or a single array + if let Ok(tuple) = xdata.cast() { + // Handle tuple case (multiple independent variables) + curve_fit_tuple(py, f, tuple, ydata, p0_array, maxfev, xtol, ftol) + } else if let Ok(array) = xdata.extract::>() { + // Handle single array case + curve_fit_array(py, f, array, ydata, p0_array, maxfev, xtol, ftol) + } else { + Err(PyErr::new::( + "xdata must be an array or tuple of arrays", + )) + } +} + +/// Helper function for `curve_fit` with single array xdata. +#[allow(clippy::type_complexity)] // Complex return type required for scipy compatibility +#[allow(clippy::too_many_arguments)] // Matches scipy.optimize.curve_fit parameters +#[allow(clippy::needless_pass_by_value)] // PyReadonlyArray1 is a lightweight wrapper +fn curve_fit_array( + py: Python<'_>, + f: Py, + xdata: PyReadonlyArray1, + ydata: PyReadonlyArray1, + p0: PyReadonlyArray1, + maxfev: Option, + xtol: Option, + ftol: Option, +) -> PyResult<(Py>, Py>)> { + let xdata_view = xdata.as_array(); + let ydata_view = ydata.as_array(); + let p0_view = p0.as_array(); + + // Create closure that calls Python function + // The Python function signature is f(x, *params) + let func = move |x: f64, params: &[f64]| -> f64 { + Python::attach(|py| { + // Build arguments tuple: (x, *params) + let mut args_vec = Vec::with_capacity(1 + params.len()); + args_vec.push(x); + args_vec.extend_from_slice(params); + + let Ok(tuple) = pyo3::types::PyTuple::new(py, &args_vec) else { + return f64::NAN; + }; + + match f.call1(py, tuple) { + Ok(result) => result.extract::(py).unwrap_or(f64::NAN), + Err(_) => f64::NAN, + } + }) + }; + + // Configure options + let opts = CurveFitOptions { + maxfev: maxfev.unwrap_or(1000), + xtol: xtol.unwrap_or(1e-8), + ftol: ftol.unwrap_or(1e-8), + lambda: 0.01, + }; + + // Call Rust implementation + let result = pecos::prelude::curve_fit(func, xdata_view, ydata_view, p0_view, Some(opts)) + .map_err(map_curve_fit_error)?; + + // Convert results to Python arrays + let popt = PyArray1::from_array(py, &result.params).unbind(); + + // If covariance is available, return it; otherwise create identity matrix + let pcov = if let Some(cov) = result.pcov { + PyArray2::from_array(py, &cov).unbind() + } else { + // Return identity matrix if covariance not available + let n = result.params.len(); + let mut cov_array = vec![vec![0.0; n]; n]; + for (i, row) in cov_array.iter_mut().enumerate().take(n) { + row[i] = 1.0; + } + PyArray2::from_vec2(py, &cov_array).unwrap().unbind() + }; + + Ok((popt, pcov)) +} + +/// Helper function for `curve_fit` with tuple of arrays as xdata. +/// +/// This handles the scipy behavior where xdata can be a tuple of arrays, +/// and the function f receives tuples of x values. +#[allow(clippy::type_complexity)] // Complex return type required for scipy compatibility +#[allow(clippy::too_many_arguments)] // Matches scipy.optimize.curve_fit parameters +#[allow(clippy::too_many_lines)] // Complex scipy compatibility logic required +#[allow(clippy::needless_pass_by_value)] // PyReadonlyArray1 is a lightweight wrapper +fn curve_fit_tuple<'py>( + py: Python<'py>, + f: Py, + xdata_tuple: &Bound<'py, PyTuple>, + ydata: PyReadonlyArray1, + p0: PyReadonlyArray1, + maxfev: Option, + xtol: Option, + ftol: Option, +) -> PyResult<(Py>, Py>)> { + // Extract arrays from tuple + let mut xdata_arrays: Vec> = Vec::new(); + for item in xdata_tuple.iter() { + // Try to extract as f64 array first + if let Ok(array) = item.extract::>() { + xdata_arrays.push(array.as_array().to_owned()); + } else if let Ok(int_array) = item.extract::>() { + // Handle integer arrays by converting to f64 + #[allow(clippy::cast_precision_loss)] + // Accepting precision loss for large integers in scientific data + let float_array: Array1 = int_array.as_array().mapv(|x| x as f64); + xdata_arrays.push(float_array); + } else if let Ok(int_array) = item.extract::>() { + // Handle i32 arrays + let float_array: Array1 = int_array.as_array().mapv(f64::from); + xdata_arrays.push(float_array); + } else { + return Err(PyErr::new::( + "Each element in xdata tuple must be a numeric array (int or float)", + )); + } + } + + if xdata_arrays.is_empty() { + return Err(PyErr::new::( + "xdata tuple must contain at least one array", + )); + } + + // Verify all arrays have the same length + let n = xdata_arrays[0].len(); + for (i, arr) in xdata_arrays.iter().enumerate().skip(1) { + if arr.len() != n { + return Err(PyErr::new::(format!( + "All xdata arrays must have the same length. Array 0 has length {}, array {} has length {}", + n, + i, + arr.len() + ))); + } + } + + let ydata_view = ydata.as_array(); + if ydata_view.len() != n { + return Err(PyErr::new::(format!( + "xdata and ydata must have the same length: xdata has {}, ydata has {}", + n, + ydata_view.len() + ))); + } + + // Create a "virtual" xdata that's just indices, and modify the function wrapper + // to look up the actual values from the tuple of arrays + #[allow(clippy::cast_precision_loss)] // Array indices are always small enough for f64 + let xdata_indices: Array1 = Array1::from_iter((0..n).map(|i| i as f64)); + + // Clone the arrays for use in closure + let xdata_arrays_clone = xdata_arrays.clone(); + + // Create closure that calls Python function with tuple of x values + // The Python function signature is f((x1, x2, ...), *params) + let func = move |idx: f64, params: &[f64]| -> f64 { + #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] + let i = idx as usize; // idx is always a valid non-negative array index + + Python::attach(|py| { + // Build tuple of x values at index i + let x_values: Vec = xdata_arrays_clone.iter().map(|arr| arr[i]).collect(); + + // Create Python tuple for x values + let Ok(x_tuple) = PyTuple::new(py, &x_values) else { + return f64::NAN; + }; + + // Build complete arguments: First create a Vec of all arguments + // Then convert to PyTuple + // Arguments are: (x_tuple, *params) + + // Create Python list to build arguments + let Ok(list_module) = py.import("builtins") else { + return f64::NAN; + }; + + let py_list = match list_module.getattr("list") { + Ok(list_func) => match list_func.call0() { + Ok(l) => l, + Err(_) => return f64::NAN, + }, + Err(_) => return f64::NAN, + }; + + // Append x_tuple as first element + if py_list.call_method1("append", (x_tuple,)).is_err() { + return f64::NAN; + } + + // Append each param + for ¶m in params { + if py_list.call_method1("append", (param,)).is_err() { + return f64::NAN; + } + } + + // Convert list to tuple + let Ok(tuple_func) = list_module.getattr("tuple") else { + return f64::NAN; + }; + + let Ok(args_tuple) = tuple_func.call1((py_list,)) else { + return f64::NAN; + }; + + // Downcast to PyTuple + let Ok(args_as_tuple) = args_tuple.cast() else { + return f64::NAN; + }; + + // Call function with arguments + match f.call1(py, args_as_tuple) { + Ok(result) => result.extract::(py).unwrap_or(f64::NAN), + Err(e) => { + let () = e.print(py); + f64::NAN + } + } + }) + }; + + // Configure options + let opts = CurveFitOptions { + maxfev: maxfev.unwrap_or(1000), + xtol: xtol.unwrap_or(1e-8), + ftol: ftol.unwrap_or(1e-8), + lambda: 0.01, + }; + + let p0_view = p0.as_array(); + + // Call Rust implementation with index-based xdata + let result = + pecos::prelude::curve_fit(func, xdata_indices.view(), ydata_view, p0_view, Some(opts)) + .map_err(map_curve_fit_error)?; + + // Convert results to Python arrays + let popt = PyArray1::from_array(py, &result.params).unbind(); + + // If covariance is available, return it; otherwise create identity matrix + let pcov = if let Some(cov) = result.pcov { + PyArray2::from_array(py, &cov).unbind() + } else { + // Return identity matrix if covariance not available + let n = result.params.len(); + let mut cov_array = vec![vec![0.0; n]; n]; + for (i, row) in cov_array.iter_mut().enumerate().take(n) { + row[i] = 1.0; + } + PyArray2::from_vec2(py, &cov_array).unwrap().unbind() + }; + + Ok((popt, pcov)) +} + +/// Register the num submodule with Python bindings. +pub fn register_num_module(m: &Bound<'_, PyModule>) -> PyResult<()> { + let num_module = PyModule::new(m.py(), "num")?; + num_module.add_function(wrap_pyfunction!(brentq, &num_module)?)?; + num_module.add_function(wrap_pyfunction!(newton, &num_module)?)?; + num_module.add_function(wrap_pyfunction!(polyfit, &num_module)?)?; + num_module.add_function(wrap_pyfunction!(curve_fit, &num_module)?)?; + num_module.add_class::()?; + m.add_submodule(&num_module)?; + Ok(()) +} diff --git a/python/pecos-rslib/src/pecos_rslib/__init__.py b/python/pecos-rslib/src/pecos_rslib/__init__.py index e8ae6b59c..66233618c 100644 --- a/python/pecos-rslib/src/pecos_rslib/__init__.py +++ b/python/pecos-rslib/src/pecos_rslib/__init__.py @@ -17,6 +17,7 @@ import ctypes import logging +import sys from importlib.metadata import PackageNotFoundError, version from pathlib import Path from typing import Any, NoReturn @@ -34,6 +35,7 @@ StateVecEngineRs, binding, # llvmlite-compatible binding module for bitcode ir, # llvmlite-compatible LLVM IR module + num, # Numerical computing functions (scipy.optimize replacements) ) from pecos_rslib.cppsparse_sim import CppSparseSimRs from pecos_rslib.rscoin_toss import CoinToss @@ -41,6 +43,9 @@ from pecos_rslib.rssparse_sim import SparseSimRs from pecos_rslib.rsstate_vec import StateVecRs +# Register num module in sys.modules to enable "from pecos_rslib.num import ..." syntax +sys.modules["pecos_rslib.num"] = num + # HUGR compilation functions - explicit, no automatic fallback try: from pecos_rslib._pecos_rslib import ( @@ -438,6 +443,8 @@ def get_compilation_backends() -> dict[str, Any]: # llvmlite-compatible modules "ir", "binding", + # Numerical computing (scipy.optimize replacements) + "num", # QuEST simulators "QuestStateVec", "QuestDensityMatrix", diff --git a/python/pecos-rslib/tests/test_scipy_comparison.py b/python/pecos-rslib/tests/test_scipy_comparison.py new file mode 100644 index 000000000..5f8a0a9d7 --- /dev/null +++ b/python/pecos-rslib/tests/test_scipy_comparison.py @@ -0,0 +1,597 @@ +""" +Comprehensive comparison tests between pecos_rslib.num and scipy.optimize. + +These tests verify that our Rust implementations produce results that match +scipy within reasonable numerical tolerances. +""" + +import numpy as np +import pytest + +# Import both our implementation and scipy +from pecos_rslib.num import brentq as pecos_brentq +from pecos_rslib.num import newton as pecos_newton +from pecos_rslib.num import curve_fit as pecos_curve_fit +from pecos_rslib.num import polyfit as pecos_polyfit +from pecos_rslib.num import Poly1d as PecosPoly1d + +from scipy.optimize import brentq as scipy_brentq +from scipy.optimize import newton as scipy_newton +from scipy.optimize import curve_fit as scipy_curve_fit + + +class TestBrentqComparison: + """Compare brentq implementations.""" + + def test_sqrt2(self): + """Find sqrt(2) by solving x^2 - 2 = 0.""" + + def f(x): + return x * x - 2.0 + + pecos_root = pecos_brentq(f, 0.0, 2.0) + scipy_root = scipy_brentq(f, 0.0, 2.0) + + assert abs(pecos_root - scipy_root) < 1e-10 + assert abs(pecos_root - np.sqrt(2)) < 1e-10 + + def test_cubic(self): + """Find root of x^3 - x - 2 = 0.""" + + def f(x): + return x**3 - x - 2.0 + + pecos_root = pecos_brentq(f, 1.0, 2.0) + scipy_root = scipy_brentq(f, 1.0, 2.0) + + assert abs(pecos_root - scipy_root) < 1e-10 + # Verify both found the correct root + assert abs(f(pecos_root)) < 1e-10 + assert abs(f(scipy_root)) < 1e-10 + + def test_transcendental(self): + """Find root of cos(x) = x.""" + + def f(x): + return np.cos(x) - x + + pecos_root = pecos_brentq(f, 0.0, 1.0) + scipy_root = scipy_brentq(f, 0.0, 1.0) + + assert abs(pecos_root - scipy_root) < 1e-10 + + def test_exponential(self): + """Find root of e^x = 3.""" + + def f(x): + return np.exp(x) - 3.0 + + pecos_root = pecos_brentq(f, 0.0, 2.0) + scipy_root = scipy_brentq(f, 0.0, 2.0) + + assert abs(pecos_root - scipy_root) < 1e-10 + assert abs(pecos_root - np.log(3)) < 1e-10 + + def test_polynomial_near_zero(self): + """Test with polynomial that has root near zero.""" + + def f(x): + return x**3 - 0.001 + + pecos_root = pecos_brentq(f, 0.0, 1.0) + scipy_root = scipy_brentq(f, 0.0, 1.0) + + assert abs(pecos_root - scipy_root) < 1e-10 + + def test_steep_function(self): + """Test with steep function.""" + + def f(x): + return np.tanh(10 * x) + + pecos_root = pecos_brentq(f, -1.0, 1.0) + scipy_root = scipy_brentq(f, -1.0, 1.0) + + assert abs(pecos_root - scipy_root) < 1e-10 + + def test_same_sign_error(self): + """Verify both implementations raise error for same-sign endpoints.""" + + def f(x): + return x * x + 1.0 # No real roots + + with pytest.raises(ValueError, match="opposite signs"): + pecos_brentq(f, -1.0, 1.0) + + with pytest.raises(ValueError, match="sign"): + scipy_brentq(f, -1.0, 1.0) + + +class TestNewtonComparison: + """Compare newton implementations.""" + + def test_sqrt2_with_derivative(self): + """Find sqrt(2) with analytical derivative.""" + + def f(x): + return x * x - 2.0 + + def fprime(x): + return 2.0 * x + + pecos_root = pecos_newton(f, 1.0, fprime=fprime) + scipy_root = scipy_newton(f, 1.0, fprime=fprime) + + assert abs(pecos_root - scipy_root) < 1e-8 + assert abs(pecos_root - np.sqrt(2)) < 1e-8 + + def test_sqrt2_numerical_derivative(self): + """Find sqrt(2) with numerical derivative.""" + + def f(x): + return x * x - 2.0 + + pecos_root = pecos_newton(f, 1.0) + scipy_root = scipy_newton(f, 1.0) + + # Numerical derivatives may differ slightly, so use larger tolerance + assert abs(pecos_root - scipy_root) < 1e-6 + assert abs(pecos_root - np.sqrt(2)) < 1e-6 + + def test_cubic_polynomial(self): + """Find root of x^3 - x - 2 = 0.""" + + def f(x): + return x**3 - x - 2.0 + + def fprime(x): + return 3.0 * x**2 - 1.0 + + pecos_root = pecos_newton(f, 1.5, fprime=fprime) + scipy_root = scipy_newton(f, 1.5, fprime=fprime) + + assert abs(pecos_root - scipy_root) < 1e-8 + + def test_exponential_function(self): + """Find root of e^x - 3 = 0.""" + + def f(x): + return np.exp(x) - 3.0 + + def fprime(x): + return np.exp(x) + + pecos_root = pecos_newton(f, 1.0, fprime=fprime) + scipy_root = scipy_newton(f, 1.0, fprime=fprime) + + assert abs(pecos_root - scipy_root) < 1e-8 + assert abs(pecos_root - np.log(3)) < 1e-8 + + def test_transcendental(self): + """Find root of cos(x) = x.""" + + def f(x): + return np.cos(x) - x + + def fprime(x): + return -np.sin(x) - 1.0 + + pecos_root = pecos_newton(f, 0.5, fprime=fprime) + scipy_root = scipy_newton(f, 0.5, fprime=fprime) + + assert abs(pecos_root - scipy_root) < 1e-8 + + def test_difficult_initial_guess(self): + """Test convergence from a non-ideal starting point.""" + + def f(x): + return x**3 - 2 * x - 5 + + def fprime(x): + return 3 * x**2 - 2 + + # Start far from the root + pecos_root = pecos_newton(f, 3.0, fprime=fprime) + scipy_root = scipy_newton(f, 3.0, fprime=fprime) + + assert abs(pecos_root - scipy_root) < 1e-8 + + +class TestCurveFitComparison: + """Compare curve_fit implementations.""" + + def test_linear_fit(self): + """Fit y = a*x + b.""" + + def linear(x, a, b): + return a * x + b + + xdata = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + ydata = np.array([1.0, 3.0, 5.0, 7.0, 9.0]) # y = 2*x + 1 + p0 = np.array([1.0, 0.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(linear, xdata, ydata, p0) + scipy_popt, scipy_pcov = scipy_curve_fit(linear, xdata, ydata, p0) + + # Parameters should match closely + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=1e-6, atol=1e-8) + # Covariances should match (may have small numerical differences) + np.testing.assert_allclose(pecos_pcov, scipy_pcov, rtol=1e-4, atol=1e-8) + + def test_exponential_fit(self): + """Fit y = a * exp(b * x).""" + + def exponential(x, a, b): + return a * np.exp(b * x) + + xdata = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + ydata = np.array([1.0, 2.718, 7.389, 20.086, 54.598]) + p0 = np.array([1.0, 1.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(exponential, xdata, ydata, p0) + scipy_popt, scipy_pcov = scipy_curve_fit(exponential, xdata, ydata, p0) + + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=1e-3, atol=1e-4) + np.testing.assert_allclose(pecos_pcov, scipy_pcov, rtol=0.1, atol=1e-6) + + def test_quadratic_fit(self): + """Fit y = a*x^2 + b*x + c.""" + + def quadratic(x, a, b, c): + return a * x**2 + b * x + c + + xdata = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + ydata = np.array([3.0, 6.0, 11.0, 18.0, 27.0]) # y = x^2 + 2*x + 3 + p0 = np.array([1.0, 1.0, 1.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(quadratic, xdata, ydata, p0) + scipy_popt, scipy_pcov = scipy_curve_fit(quadratic, xdata, ydata, p0) + + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(pecos_pcov, scipy_pcov, rtol=1e-4, atol=1e-8) + + def test_gaussian_fit(self): + """Fit Gaussian function.""" + + def gaussian(x, amp, mu, sigma): + return amp * np.exp(-((x - mu) ** 2) / (2 * sigma**2)) + + xdata = np.linspace(-5, 5, 50) + ydata = gaussian(xdata, 2.0, 1.0, 1.5) + 0.01 * np.random.randn(50) + p0 = np.array([1.0, 0.0, 1.0]) + + # Set random seed for reproducibility + np.random.seed(42) + ydata = gaussian(xdata, 2.0, 1.0, 1.5) + 0.01 * np.random.randn(50) + + pecos_popt, pecos_pcov = pecos_curve_fit( + gaussian, xdata, ydata, p0, maxfev=5000 + ) + scipy_popt, scipy_pcov = scipy_curve_fit( + gaussian, xdata, ydata, p0, maxfev=5000 + ) + + # Parameters should be close (some variation due to optimization differences) + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=0.1, atol=0.1) + + def test_tuple_xdata_quantum_pecos_pattern(self): + """Test curve_fit with tuple xdata (quantum-pecos pattern).""" + + def func(x, a, b, c): + p, d = x + return a * p ** (b / d + c) + + p = np.array([0.001, 0.002, 0.003, 0.004, 0.005]) + d = np.array([3, 3, 3, 3, 3]) # Integer array! + plog = np.array([0.01, 0.015, 0.02, 0.025, 0.03]) + p0 = np.array([1.0, 1.0, 1.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(func, (p, d), plog, p0, maxfev=5000) + scipy_popt, scipy_pcov = scipy_curve_fit(func, (p, d), plog, p0, maxfev=5000) + + # This is a difficult optimization problem - different optimizers may converge + # to different local minima. Instead of comparing parameters, verify both + # implementations produce good fits to the data. + pecos_residuals = [] + scipy_residuals = [] + for i in range(len(p)): + pecos_pred = func((p[i], d[i]), *pecos_popt) + scipy_pred = func((p[i], d[i]), *scipy_popt) + pecos_residuals.append((pecos_pred - plog[i]) ** 2) + scipy_residuals.append((scipy_pred - plog[i]) ** 2) + + pecos_rmse = np.sqrt(np.mean(pecos_residuals)) + scipy_rmse = np.sqrt(np.mean(scipy_residuals)) + + # Both should produce reasonable fits + assert pecos_rmse < 0.01, f"PECOS fit too poor: RMSE={pecos_rmse}" + assert scipy_rmse < 0.01, f"Scipy fit too poor: RMSE={scipy_rmse}" + # And similar fit quality + assert abs(pecos_rmse - scipy_rmse) < 0.005, "Fit quality differs too much" + + def test_sine_fit(self): + """Fit sine wave.""" + + def sine(x, amp, freq, phase): + return amp * np.sin(2 * np.pi * freq * x + phase) + + xdata = np.linspace(0, 2, 100) + np.random.seed(42) + ydata = sine(xdata, 1.5, 2.0, 0.5) + 0.05 * np.random.randn(100) + p0 = np.array([1.0, 2.0, 0.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(sine, xdata, ydata, p0, maxfev=5000) + scipy_popt, scipy_pcov = scipy_curve_fit(sine, xdata, ydata, p0, maxfev=5000) + + # Parameters should be similar + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=0.1, atol=0.1) + + def test_power_law_fit(self): + """Fit power law y = a * x^b.""" + + def power_law(x, a, b): + return a * x**b + + xdata = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + ydata = np.array([2.0, 5.66, 10.39, 16.0, 22.36]) # y ≈ 2*x^1.5 + p0 = np.array([1.0, 1.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(power_law, xdata, ydata, p0) + scipy_popt, scipy_pcov = scipy_curve_fit(power_law, xdata, ydata, p0) + + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=1e-3, atol=1e-4) + + def test_noisy_data(self): + """Test with noisy data.""" + + def linear(x, a, b): + return a * x + b + + np.random.seed(123) + xdata = np.linspace(0, 10, 50) + ydata = 2.5 * xdata + 1.3 + np.random.normal(0, 0.5, 50) + p0 = np.array([1.0, 0.0]) + + pecos_popt, pecos_pcov = pecos_curve_fit(linear, xdata, ydata, p0) + scipy_popt, scipy_pcov = scipy_curve_fit(linear, xdata, ydata, p0) + + # Should converge to similar values + np.testing.assert_allclose(pecos_popt, scipy_popt, rtol=1e-4, atol=1e-6) + np.testing.assert_allclose(pecos_pcov, scipy_pcov, rtol=0.01, atol=1e-8) + + def test_p0_accepts_sequence_types(self): + """Test that p0 accepts tuple, list, and array (scipy compatibility).""" + + def quadratic(x, a, b, c): + return a * x**2 + b * x + c + + xdata = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + ydata = np.array([1.0, 2.0, 5.0, 10.0, 17.0]) # y = x^2 + 1 + + # Test with tuple (quantum-pecos usage pattern) + p0_tuple = (1.0, 0.0, 0.0) + popt_tuple, _ = pecos_curve_fit(quadratic, xdata, ydata, p0_tuple) + + # Test with list + p0_list = [1.0, 0.0, 0.0] + popt_list, _ = pecos_curve_fit(quadratic, xdata, ydata, p0_list) + + # Test with array + p0_array = np.array([1.0, 0.0, 0.0]) + popt_array, _ = pecos_curve_fit(quadratic, xdata, ydata, p0_array) + + # All should produce the same result + np.testing.assert_allclose(popt_tuple, popt_array, rtol=1e-10, atol=1e-12) + np.testing.assert_allclose(popt_list, popt_array, rtol=1e-10, atol=1e-12) + + # Should match expected values + np.testing.assert_allclose(popt_array, [1.0, 0.0, 1.0], rtol=1e-6, atol=1e-8) + + +class TestPolyfitComparison: + """Compare polyfit implementations.""" + + def test_linear_fit(self): + """Fit degree 1 polynomial (line).""" + x = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + y = np.array([1.0, 3.0, 5.0, 7.0, 9.0]) # y = 2*x + 1 + + pecos_coeffs = pecos_polyfit(x, y, 1) + scipy_coeffs = np.polyfit(x, y, 1) + + np.testing.assert_allclose(pecos_coeffs, scipy_coeffs, rtol=1e-10, atol=1e-12) + + def test_quadratic_fit(self): + """Fit degree 2 polynomial.""" + x = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + y = np.array([3.0, 6.0, 11.0, 18.0, 27.0]) # y = x^2 + 2*x + 3 + + pecos_coeffs = pecos_polyfit(x, y, 2) + scipy_coeffs = np.polyfit(x, y, 2) + + np.testing.assert_allclose(pecos_coeffs, scipy_coeffs, rtol=1e-10, atol=1e-12) + + def test_cubic_fit(self): + """Fit degree 3 polynomial.""" + x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) + y = np.array([1.0, 3.0, 17.0, 55.0, 129.0, 251.0]) # y = x^3 + 2*x^2 + 3*x + 1 + + pecos_coeffs = pecos_polyfit(x, y, 3) + scipy_coeffs = np.polyfit(x, y, 3) + + np.testing.assert_allclose(pecos_coeffs, scipy_coeffs, rtol=1e-9, atol=1e-10) + + def test_high_degree(self): + """Test higher degree polynomial.""" + np.random.seed(42) + x = np.linspace(0, 1, 20) + # Generate y from a degree 5 polynomial + true_coeffs = np.array([1.0, -2.0, 3.0, -1.0, 2.0, 1.0]) + y = np.polyval(true_coeffs, x) + + pecos_coeffs = pecos_polyfit(x, y, 5) + scipy_coeffs = np.polyfit(x, y, 5) + + np.testing.assert_allclose(pecos_coeffs, scipy_coeffs, rtol=1e-8, atol=1e-10) + # Verify we recovered the original coefficients + np.testing.assert_allclose(pecos_coeffs, true_coeffs, rtol=1e-8, atol=1e-10) + + def test_noisy_data(self): + """Test polyfit with noisy data.""" + np.random.seed(456) + x = np.linspace(0, 5, 30) + y = 2 * x**2 - 3 * x + 1 + np.random.normal(0, 0.5, 30) + + pecos_coeffs = pecos_polyfit(x, y, 2) + scipy_coeffs = np.polyfit(x, y, 2) + + # Should get similar coefficients + np.testing.assert_allclose(pecos_coeffs, scipy_coeffs, rtol=1e-8, atol=1e-10) + + def test_overdetermined_system(self): + """Test with many more data points than parameters.""" + np.random.seed(789) + x = np.linspace(-2, 2, 100) + y = 1.5 * x + 2.0 + np.random.normal(0, 0.1, 100) + + pecos_coeffs = pecos_polyfit(x, y, 1) + scipy_coeffs = np.polyfit(x, y, 1) + + np.testing.assert_allclose(pecos_coeffs, scipy_coeffs, rtol=1e-8, atol=1e-10) + + +class TestPoly1dComparison: + """Compare Poly1d implementations.""" + + def test_evaluation(self): + """Test polynomial evaluation.""" + coeffs = np.array([2.0, 3.0, 1.0]) # 2*x^2 + 3*x + 1 + + pecos_poly = PecosPoly1d(coeffs) + scipy_poly = np.poly1d(coeffs) + + test_points = [-2.0, -1.0, 0.0, 1.0, 2.0, 3.5] + for x in test_points: + pecos_val = pecos_poly.eval(x) + scipy_val = scipy_poly(x) + assert abs(pecos_val - scipy_val) < 1e-12 + + def test_degree(self): + """Test degree calculation.""" + coeffs = np.array([1.0, 2.0, 3.0, 4.0]) # degree 3 + + pecos_poly = PecosPoly1d(coeffs) + scipy_poly = np.poly1d(coeffs) + + assert pecos_poly.degree() == len(coeffs) - 1 + assert pecos_poly.degree() == scipy_poly.order + + def test_fit_and_evaluate(self): + """Test fitting then evaluating.""" + x = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + y = np.array([1.0, 3.0, 5.0, 7.0, 9.0]) + + pecos_coeffs = pecos_polyfit(x, y, 1) + scipy_coeffs = np.polyfit(x, y, 1) + + pecos_poly = PecosPoly1d(pecos_coeffs) + scipy_poly = np.poly1d(scipy_coeffs) + + # Evaluate at original points + for xi, yi in zip(x, y, strict=False): + pecos_val = pecos_poly.eval(xi) + scipy_val = scipy_poly(xi) + assert abs(pecos_val - scipy_val) < 1e-10 + assert abs(pecos_val - yi) < 1e-10 + + def test_complex_polynomial(self): + """Test with complex polynomial.""" + coeffs = np.array([1.0, -2.5, 3.7, -1.2, 0.5]) + + pecos_poly = PecosPoly1d(coeffs) + scipy_poly = np.poly1d(coeffs) + + test_points = np.linspace(-3, 3, 20) + for x in test_points: + pecos_val = pecos_poly.eval(x) + scipy_val = scipy_poly(x) + assert abs(pecos_val - scipy_val) < 1e-10 + + +class TestEdgeCases: + """Test edge cases and error handling.""" + + def test_brentq_narrow_interval(self): + """Test brentq with very narrow interval.""" + + def f(x): + return x - 0.5 + + pecos_root = pecos_brentq(f, 0.4999, 0.5001) + scipy_root = scipy_brentq(f, 0.4999, 0.5001) + + assert abs(pecos_root - scipy_root) < 1e-10 + + def test_newton_near_zero_derivative(self): + """Test newton with function that has small derivative. + + Note: This is a pathological case where x^3 has a triple root at 0. + Newton's method may struggle to converge precisely due to the flat derivative. + """ + + def f(x): + return x**3 + + def fprime(x): + return 3 * x**2 + + # Both should converge to something close to 0 + pecos_root = pecos_newton(f, 0.1, fprime=fprime) + scipy_root = scipy_newton(f, 0.1, fprime=fprime) + + # Verify both find a root (may not be exactly 0 due to numerical issues) + assert ( + abs(f(pecos_root)) < 1e-6 + ), f"PECOS didn't find root: f({pecos_root})={f(pecos_root)}" + assert ( + abs(f(scipy_root)) < 1e-6 + ), f"Scipy didn't find root: f({scipy_root})={f(scipy_root)}" + + # Both should be close to 0 (allow larger tolerance due to pathological case) + assert abs(pecos_root) < 0.01, f"PECOS root too far from 0: {pecos_root}" + assert abs(scipy_root) < 0.01, f"Scipy root too far from 0: {scipy_root}" + + def test_curve_fit_exact_fit(self): + """Test curve_fit with data that fits model exactly.""" + + def linear(x, a, b): + return a * x + b + + xdata = np.array([0.0, 1.0, 2.0]) + ydata = np.array([1.0, 3.0, 5.0]) # Exactly y = 2*x + 1 + p0 = np.array([1.0, 0.0]) + + pecos_popt, _ = pecos_curve_fit(linear, xdata, ydata, p0) + scipy_popt, _ = scipy_curve_fit(linear, xdata, ydata, p0) + + # Should get exact solution + np.testing.assert_allclose(pecos_popt, [2.0, 1.0], rtol=1e-10, atol=1e-12) + np.testing.assert_allclose(scipy_popt, [2.0, 1.0], rtol=1e-10, atol=1e-12) + + def test_polyfit_exact_degree(self): + """Test polyfit when data is exact polynomial.""" + # Generate data from exact polynomial + x = np.array([0.0, 1.0, 2.0, 3.0]) + coeffs_true = np.array([2.0, -1.0, 3.0]) # 2*x^2 - x + 3 + y = np.polyval(coeffs_true, x) + + pecos_coeffs = pecos_polyfit(x, y, 2) + scipy_coeffs = np.polyfit(x, y, 2) + + # Should recover exact coefficients + np.testing.assert_allclose(pecos_coeffs, coeffs_true, rtol=1e-12, atol=1e-14) + np.testing.assert_allclose(scipy_coeffs, coeffs_true, rtol=1e-12, atol=1e-14) + + +if __name__ == "__main__": + pytest.main([__file__, "-v"]) diff --git a/python/quantum-pecos/pyproject.toml b/python/quantum-pecos/pyproject.toml index 63d45cf6c..92366820d 100644 --- a/python/quantum-pecos/pyproject.toml +++ b/python/quantum-pecos/pyproject.toml @@ -31,7 +31,6 @@ dependencies = [ "pecos-rslib==0.7.0.dev4", "phir>=0.3.3", "numpy>=1.15.0", - "scipy>=1.1.0", "networkx>=2.1.0", ] classifiers = [ diff --git a/python/quantum-pecos/src/pecos/misc/threshold_curve.py b/python/quantum-pecos/src/pecos/misc/threshold_curve.py index 517c0838d..5c73fdae1 100644 --- a/python/quantum-pecos/src/pecos/misc/threshold_curve.py +++ b/python/quantum-pecos/src/pecos/misc/threshold_curve.py @@ -22,7 +22,7 @@ from typing import TYPE_CHECKING import numpy as np -from scipy.optimize import curve_fit +from pecos_rslib.num import curve_fit if TYPE_CHECKING: from collections.abc import Callable diff --git a/python/quantum-pecos/src/pecos/tools/pseudo_threshold_tools.py b/python/quantum-pecos/src/pecos/tools/pseudo_threshold_tools.py index e0a0ec92e..3ed28b5ce 100644 --- a/python/quantum-pecos/src/pecos/tools/pseudo_threshold_tools.py +++ b/python/quantum-pecos/src/pecos/tools/pseudo_threshold_tools.py @@ -23,7 +23,7 @@ from typing import TYPE_CHECKING import numpy as np -from scipy.optimize import brentq, curve_fit, newton +from pecos_rslib.num import brentq, curve_fit, newton from pecos.decoders import MWPM2D from pecos.engines import circuit_runners diff --git a/uv.lock b/uv.lock index 061bed69a..5dfb7cf14 100644 --- a/uv.lock +++ b/uv.lock @@ -2621,11 +2621,23 @@ source = { editable = "python/pecos-rslib" } dev = [ { name = "patchelf", marker = "sys_platform != 'win32'" }, ] +test = [ + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "numpy", version = "2.3.4", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "pytest" }, + { name = "scipy", version = "1.15.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "scipy", version = "1.16.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, +] [package.metadata] [package.metadata.requires-dev] dev = [{ name = "patchelf", marker = "sys_platform != 'win32'" }] +test = [ + { name = "numpy", specifier = ">=1.20" }, + { name = "pytest", specifier = ">=7.0" }, + { name = "scipy", specifier = ">=1.7" }, +] [[package]] name = "pecos-workspace" @@ -3521,8 +3533,6 @@ dependencies = [ { name = "numpy", version = "2.3.4", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "pecos-rslib" }, { name = "phir" }, - { name = "scipy", version = "1.15.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, - { name = "scipy", version = "1.16.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, ] [package.optional-dependencies] @@ -3560,7 +3570,6 @@ requires-dist = [ { name = "pytket-cutensornet", marker = "python_full_version >= '3.11' and extra == 'cuda'", specifier = ">=0.12.0" }, { name = "quantum-pecos", extras = ["guppy"], marker = "extra == 'all'" }, { name = "quantum-pecos", extras = ["visualization"], marker = "extra == 'all'" }, - { name = "scipy", specifier = ">=1.1.0" }, { name = "selene-sim", marker = "extra == 'guppy'", specifier = "~=0.2.0" }, ] provides-extras = ["guppy", "visualization", "all", "cuda"]