diff --git a/.cargo/config.toml b/.cargo/config.toml new file mode 100644 index 0000000..b4692be --- /dev/null +++ b/.cargo/config.toml @@ -0,0 +1,13 @@ +# PyO3's `extension-module` feature (the default for this crate, see +# Cargo.toml) tells PyO3 not to link against libpython directly -- those +# symbols are meant to be resolved by the Python interpreter at import +# time, the way `maturin` builds and loads this crate. Plain `cargo +# build`/`cargo test` on macOS don't get maturin's linker setup, so +# Apple's `ld` fails on the unresolved Python C-API symbols unless told to +# defer them to runtime, same as `maturin` already does under the hood. +# See the PyO3 FAQ: https://pyo3.rs/main/faq +[target.x86_64-apple-darwin] +rustflags = ["-C", "link-arg=-undefined", "-C", "link-arg=dynamic_lookup"] + +[target.aarch64-apple-darwin] +rustflags = ["-C", "link-arg=-undefined", "-C", "link-arg=dynamic_lookup"] \ No newline at end of file diff --git a/.gitignore b/.gitignore index 8130c3a..fb0329c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ debug/ target/ +.DS_Store \ No newline at end of file diff --git a/Cargo.lock b/Cargo.lock index 3922a32..32bf657 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -20,6 +20,15 @@ dependencies = [ "cc", ] +[[package]] +name = "android_system_properties" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae221649c9976a6f6c56ae1facf410f3ddb33cc661c4b7b61020a912d4237fbc" +dependencies = [ + "libc", +] + [[package]] name = "anes" version = "0.1.6" @@ -47,12 +56,24 @@ version = "1.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" +[[package]] +name = "bitflags" +version = "2.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06" + [[package]] name = "bumpalo" version = "3.20.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649" +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" + [[package]] name = "cast" version = "0.3.0" @@ -75,6 +96,31 @@ version = "1.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" +[[package]] +name = "chacha20" +version = "0.10.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65c35e4b699c7e15ccbe7ee35c005e4fc0a278d22238a2857e6ce2dadeda1b06" +dependencies = [ + "cfg-if", + "cpufeatures", + "rand_core 0.10.1", +] + +[[package]] +name = "chrono" +version = "0.4.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1aa79e62e7697b8e29b513a68abacf485adcd1fe8284a4316c5ae868e6633327" +dependencies = [ + "iana-time-zone", + "js-sys", + "num-traits", + "serde", + "wasm-bindgen", + "windows-link", +] + [[package]] name = "ciborium" version = "0.2.2" @@ -127,6 +173,30 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9" +[[package]] +name = "core-foundation-sys" +version = "0.8.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" + +[[package]] +name = "cpufeatures" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ca28b0ae3115b884660db4118d803791fd6756b6e88f39c0f3f7859060d7566" +dependencies = [ + "libc", +] + +[[package]] +name = "crc32fast" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8498c871161e1742aaa9d52551b2d6ebdd4c3d45a3be423e3728f33b955be550" +dependencies = [ + "cfg-if", +] + [[package]] name = "criterion" version = "0.8.2" @@ -229,6 +299,17 @@ dependencies = [ "slab", ] +[[package]] +name = "getrandom" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + [[package]] name = "getrandom" version = "0.3.4" @@ -236,22 +317,67 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" dependencies = [ "cfg-if", + "js-sys", "libc", - "r-efi", + "r-efi 5.3.0", "wasip2", + "wasm-bindgen", +] + +[[package]] +name = "getrandom" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "300e883d756b2e4ec94e02791f39b04b522276138852cfc41d9fb7e904106099" +dependencies = [ + "cfg-if", + "js-sys", + "libc", + "r-efi 6.0.0", + "rand_core 0.10.1", + "wasm-bindgen", ] +[[package]] +name = "glam" +version = "0.30.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19fc433e8437a212d1b6f1e68c7824af3aed907da60afa994e7f542d18d12aa9" + +[[package]] +name = "glam" +version = "0.31.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "556f6b2ea90b8d15a74e0e7bb41671c9bdf38cd9f78c284d750b9ce58a2b5be7" + +[[package]] +name = "glam" +version = "0.32.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f70749695b063ecbf6b62949ccccde2e733ec3ecbbd71d467dca4e5c6c97cca0" + +[[package]] +name = "glam" +version = "0.33.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e762b9b188634fc508d4ec54b62ea3d352c43fd431d18d05d46f559f217f4725" + [[package]] name = "gstools-core" -version = "1.3.0" +version = "1.4.0" dependencies = [ "approx", "criterion", + "kdtree", + "nalgebra 0.35.0", "ndarray", "ndarray-rand", "numpy", "pyo3", "rayon", + "scirs2-special", + "serde", + "serde_json", ] [[package]] @@ -271,6 +397,30 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "iana-time-zone" +version = "0.1.65" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e31bc9ad994ba00e440a8aa5c9ef0ec67d5cb5e5cb0cc7f8b744a35b389cc470" +dependencies = [ + "android_system_properties", + "core-foundation-sys", + "iana-time-zone-haiku", + "js-sys", + "log", + "wasm-bindgen", + "windows-core", +] + +[[package]] +name = "iana-time-zone-haiku" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f31827a206f56af32e590ba56d5d2d085f558508192593743f16b2306495269f" +dependencies = [ + "cc", +] + [[package]] name = "itertools" version = "0.13.0" @@ -298,18 +448,59 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "kdtree" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9cd666c2aae5dde60d2f9bfa4e1f8c17698fce13d599c0f9b76725641fec13d9" +dependencies = [ + "num-traits", + "thiserror", +] + +[[package]] +name = "lazy_static" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" + [[package]] name = "libc" version = "0.2.186" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "68ab91017fe16c622486840e4c83c9a37afeff978bd239b5293d61ece587de66" +[[package]] +name = "libloading" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "754ca22de805bb5744484a5b151a9e1a8e837d5dc232c2d7d8c2e3492edc8b60" +dependencies = [ + "cfg-if", + "windows-link", +] + [[package]] name = "libm" version = "0.2.16" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" +[[package]] +name = "lock_api" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "224399e74b87b5f3557511d98dff8b14089b3dadafcab6bb93eab67d3aace965" +dependencies = [ + "scopeguard", +] + +[[package]] +name = "log" +version = "0.4.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f9f8bd3e56ce4dfc153cf470fffbfa98c7620958b312ca5c3a4b8d5181fd13c6" + [[package]] name = "matrixmultiply" version = "0.3.10" @@ -326,6 +517,54 @@ version = "2.8.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6b947ae49db0d222b1dbc6b113ce7248a3fc3a6ca21b696717bfc000ba4484d8" +[[package]] +name = "nalgebra" +version = "0.33.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9d43ddcacf343185dfd6de2ee786d9e8b1c2301622afab66b6c73baf9882abfd" +dependencies = [ + "approx", + "matrixmultiply", + "num-complex", + "num-rational", + "num-traits", + "rand 0.8.8", + "rand_distr 0.4.3", + "simba 0.9.1", + "typenum", +] + +[[package]] +name = "nalgebra" +version = "0.35.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adc43a60c217b0c6ff46e47f26911015ad8d2e5a8be1af668c67e370d99a4346" +dependencies = [ + "approx", + "glam 0.30.10", + "glam 0.31.1", + "glam 0.32.1", + "glam 0.33.7", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba 0.10.2", + "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.117", +] + [[package]] name = "ndarray" version = "0.17.2" @@ -341,6 +580,7 @@ dependencies = [ "portable-atomic-util", "rawpointer", "rayon", + "serde", ] [[package]] @@ -350,8 +590,18 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "180f724d496e84764e8ecf28fbe1da74ef231ec4ba15be65a9100be8445d73e3" dependencies = [ "ndarray", - "rand", - "rand_distr", + "rand 0.9.4", + "rand_distr 0.5.1", +] + +[[package]] +name = "num-bigint" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367" +dependencies = [ + "num-integer", + "num-traits", ] [[package]] @@ -372,6 +622,17 @@ dependencies = [ "num-traits", ] +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + [[package]] name = "num-traits" version = "0.2.19" @@ -420,6 +681,35 @@ dependencies = [ "winapi", ] +[[package]] +name = "parking_lot" +version = "0.12.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93857453250e3077bd71ff98b6a65ea6621a19bb0f559a85248955ac12c45a1a" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2621685985a2ebf1c516881c026032ac7deafcda1a2c9b7850dc81e3dfcb64c1" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-link", +] + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + [[package]] name = "pin-project-lite" version = "0.2.17" @@ -529,7 +819,7 @@ dependencies = [ "proc-macro2", "pyo3-macros-backend", "quote", - "syn", + "syn 2.0.117", ] [[package]] @@ -542,7 +832,7 @@ dependencies = [ "proc-macro2", "pyo3-build-config", "quote", - "syn", + "syn 2.0.117", ] [[package]] @@ -560,14 +850,52 @@ version = "5.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" +[[package]] +name = "r-efi" +version = "6.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" + +[[package]] +name = "rand" +version = "0.8.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e058c7de0b26af77780c769414d6257830bb240f3c38477dbc2c16e5f54d6d4c" +dependencies = [ + "libc", + "rand_chacha 0.3.1", + "rand_core 0.6.4", +] + [[package]] name = "rand" version = "0.9.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "44c5af06bb1b7d3216d91932aed5265164bf384dc89cd6ba05cf59a35f5f76ea" dependencies = [ - "rand_chacha", - "rand_core", + "rand_chacha 0.9.0", + "rand_core 0.9.5", +] + +[[package]] +name = "rand" +version = "0.10.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7f5fa3a058cd35567ef9bfa5e75732bee0f9e4c55fa90477bef2dfcdbc4be80" +dependencies = [ + "chacha20", + "getrandom 0.4.3", + "rand_core 0.10.1", +] + +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core 0.6.4", ] [[package]] @@ -577,7 +905,26 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" dependencies = [ "ppv-lite86", - "rand_core", + "rand_core 0.9.5", +] + +[[package]] +name = "rand_chacha" +version = "0.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3e6af7f3e25ded52c41df4e0b1af2d047e45896c2f3281792ed68a1c243daedb" +dependencies = [ + "ppv-lite86", + "rand_core 0.10.1", +] + +[[package]] +name = "rand_core" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom 0.2.17", ] [[package]] @@ -586,7 +933,23 @@ version = "0.9.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c" dependencies = [ - "getrandom", + "getrandom 0.3.4", +] + +[[package]] +name = "rand_core" +version = "0.10.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "63b8176103e19a2643978565ca18b50549f6101881c443590420e4dc998a3c69" + +[[package]] +name = "rand_distr" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31" +dependencies = [ + "num-traits", + "rand 0.8.8", ] [[package]] @@ -596,7 +959,17 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6a8615d50dcf34fa31f7ab52692afec947c4dd0ab803cc87cb3b0b4570ff7463" dependencies = [ "num-traits", - "rand", + "rand 0.9.4", +] + +[[package]] +name = "rand_distr" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4d431c2703ccf129de4d45253c03f49ebb22b97d6ad79ee3ecfc7e3f4862c1d8" +dependencies = [ + "num-traits", + "rand 0.10.2", ] [[package]] @@ -625,11 +998,20 @@ dependencies = [ "crossbeam-utils", ] +[[package]] +name = "redox_syscall" +version = "0.5.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed2bf2547551a7053d6fdfafda3f938979645c44812fbfcda098faae3f1a362d" +dependencies = [ + "bitflags", +] + [[package]] name = "regex" -version = "1.12.3" +version = "1.13.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e10754a14b9137dd7b1e3e5b0493cc9171fdd105e0ab477f51b72e7f3ac0e276" +checksum = "f020237b6c8eed93db2e2cb53c00c60a8e1bc73da7d073199a1180401450218d" dependencies = [ "aho-corasick", "memchr", @@ -639,9 +1021,9 @@ dependencies = [ [[package]] name = "regex-automata" -version = "0.4.14" +version = "0.4.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6e1dd4122fc1595e8162618945476892eefca7b88c52820e74af6262213cae8f" +checksum = "ad8553b9b26413251cbf30e620595c7a41b3887f03da04579c0e6b0d6a06b4b2" dependencies = [ "aho-corasick", "memchr", @@ -650,9 +1032,9 @@ dependencies = [ [[package]] name = "regex-syntax" -version = "0.8.10" +version = "0.8.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dc897dd8d9e8bd1ed8cdad82b5966c3e0ecae09fb1907d58efaa013543185d0a" +checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" [[package]] name = "rustc-hash" @@ -666,6 +1048,24 @@ version = "1.0.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d" +[[package]] +name = "safe_arch" +version = "0.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "safe_arch" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42c6efa15875e6ecb39ca61fb0b0c1a40b84fac5a5ffe71eef7d1000c8eb3f5f" +dependencies = [ + "bytemuck", +] + [[package]] name = "same-file" version = "1.0.6" @@ -675,11 +1075,62 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "scirs2-core" +version = "0.6.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f7d5e6d0c09bd4a0307e073f6526d8b7b534e0fd9cd9ee644a73e169e08b45bd" +dependencies = [ + "chrono", + "crc32fast", + "getrandom 0.3.4", + "getrandom 0.4.3", + "libloading", + "ndarray", + "ndarray-rand", + "num-complex", + "num-integer", + "num-traits", + "once_cell", + "parking_lot", + "rand 0.10.2", + "rand_chacha 0.10.0", + "rand_distr 0.6.0", + "regex", + "statrs", + "thiserror", + "uuid", + "wide 1.7.0", + "windows-sys", +] + +[[package]] +name = "scirs2-special" +version = "0.6.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94350d85540394d20e3eec505052ed052df28c7c0e0d77bebc02d7a9533d2670" +dependencies = [ + "approx", + "lazy_static", + "log", + "regex", + "scirs2-core", + "serde", + "statrs", + "thiserror", +] + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + [[package]] name = "serde" -version = "1.0.228" +version = "1.0.229" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9a8e94ea7f378bd32cbbd37198a4a91436180c5bb472411e48b5ec2e2124ae9e" +checksum = "4148590afebada386688f18773da617792bf2ef03ffc1e4cbd2b1d45b023e0ba" dependencies = [ "serde_core", "serde_derive", @@ -687,29 +1138,29 @@ dependencies = [ [[package]] name = "serde_core" -version = "1.0.228" +version = "1.0.229" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "41d385c7d4ca58e59fc732af25c3983b67ac852c1a25000afe1175de458b67ad" +checksum = "67dca2c9c51e58a4791a4b1ed58308b39c64224d349a935ab5039aa360942a48" dependencies = [ "serde_derive", ] [[package]] name = "serde_derive" -version = "1.0.228" +version = "1.0.229" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79" +checksum = "e7a5d71263a5a7d47b41f6b3f06ba276f10cc18b0931f1799f710578e2309348" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 3.0.5", ] [[package]] name = "serde_json" -version = "1.0.150" +version = "1.0.151" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e8014e44b4736ed0538adeecded0fce2a272f22dc9578a7eb6b2d9993c74cfb9" +checksum = "c841b55ecdae098c80dcae9cf767f6f8a0c2cdb3416bbef72181df4d0fe73f14" dependencies = [ "itoa", "memchr", @@ -724,12 +1175,55 @@ version = "2.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8fadd59c855ef2080decdef8ff161eb6661b86933c9d82e5ba29dc602a55aba" +[[package]] +name = "simba" +version = "0.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c99284beb21666094ba2b75bbceda012e610f5479dfcc2d6e2426f53197ffd95" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide 0.7.33", +] + +[[package]] +name = "simba" +version = "0.10.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1a7200d82ff1c7b4235efd12f11e2537a402c91cf83a5cb97ce80eef787fde7" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "wide 1.7.0", +] + [[package]] name = "slab" version = "0.4.12" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5" +[[package]] +name = "smallvec" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9be42f50aa861c555654aa3a37f52f4b1074bacf4e48fe0ef7fa584e80f1f0f" + +[[package]] +name = "statrs" +version = "0.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2a3fe7c28c6512e766b0874335db33c94ad7b8f9054228ae1c2abd47ce7d335e" +dependencies = [ + "approx", + "nalgebra 0.33.3", + "num-traits", + "rand 0.8.8", +] + [[package]] name = "syn" version = "2.0.117" @@ -741,12 +1235,43 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "syn" +version = "3.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "12df2e0110f65b775f769bb17ef989067a1d931b2eb822bd4346631eeada89f9" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + [[package]] name = "target-lexicon" version = "0.13.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "adb6935a6f5c20170eeceb1a3835a49e12e19d792f6dd344ccc76a985ca5a6ca" +[[package]] +name = "thiserror" +version = "2.0.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec86235f5fcc2a73650310756d2ac5b138a5780bbbdfae3eeccec992c435ba4f" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "2.0.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bc04cd3e1236dd4a98afca4569f2deb3f120e5422a4023be2cb683f8486292af" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.5", +] + [[package]] name = "tinytemplate" version = "1.2.1" @@ -757,12 +1282,30 @@ dependencies = [ "serde_json", ] +[[package]] +name = "typenum" +version = "1.20.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6f5e870be6c3b371b77fe0ee0bafb859fa4964b4404c27de1d380043c4dda20" + [[package]] name = "unicode-ident" version = "1.0.24" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75" +[[package]] +name = "uuid" +version = "1.26.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2ef6dac1e96601b4fb3acccccff2139741fcb757cb9a36089bf5be91cfb285ce" +dependencies = [ + "getrandom 0.4.3", + "js-sys", + "serde_core", + "wasm-bindgen", +] + [[package]] name = "walkdir" version = "2.5.0" @@ -773,6 +1316,12 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "wasi" +version = "0.11.1+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b" + [[package]] name = "wasip2" version = "1.0.3+wasi-0.2.9" @@ -814,7 +1363,7 @@ dependencies = [ "bumpalo", "proc-macro2", "quote", - "syn", + "syn 2.0.117", "wasm-bindgen-shared", ] @@ -837,6 +1386,26 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "wide" +version = "0.7.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" +dependencies = [ + "bytemuck", + "safe_arch 0.7.4", +] + +[[package]] +name = "wide" +version = "1.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8cf05ca94c9fba0c51316899caa1d1e74a064fc310a5efa7375e614343d415be" +dependencies = [ + "bytemuck", + "safe_arch 1.2.0", +] + [[package]] name = "winapi" version = "0.3.9" @@ -868,12 +1437,65 @@ version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" +[[package]] +name = "windows-core" +version = "0.62.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b8e83a14d34d0623b51dce9581199302a221863196a1dde71a7663a4c2be9deb" +dependencies = [ + "windows-implement", + "windows-interface", + "windows-link", + "windows-result", + "windows-strings", +] + +[[package]] +name = "windows-implement" +version = "0.60.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "053e2e040ab57b9dc951b72c264860db7eb3b0200ba345b4e4c3b14f67855ddf" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.117", +] + +[[package]] +name = "windows-interface" +version = "0.59.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f316c4a2570ba26bbec722032c4099d8c8bc095efccdc15688708623367e358" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.117", +] + [[package]] name = "windows-link" version = "0.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" +[[package]] +name = "windows-result" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7781fa89eaf60850ac3d2da7af8e5242a5ea78d1a11c49bf2910bb5a73853eb5" +dependencies = [ + "windows-link", +] + +[[package]] +name = "windows-strings" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7837d08f69c77cf6b07689544538e017c1bfcf57e34b4c0ff58e6c2cd3b37091" +dependencies = [ + "windows-link", +] + [[package]] name = "windows-sys" version = "0.61.2" @@ -906,7 +1528,7 @@ checksum = "8fd425244944f4ab65ccff928e7323354c5a018c75838362fdce749dfad2ee1e" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.117", ] [[package]] diff --git a/Cargo.toml b/Cargo.toml index 5a40838..1aab800 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "gstools-core" -version = "1.3.0" +version = "1.4.0" authors = ["Lennart Schüler "] edition = "2021" description = "The core functions of GSTools" @@ -33,6 +33,11 @@ pyo3 = { version = "0.28", features = ["abi3-py38", "extension-module"] } numpy = "0.28" ndarray = { version = "0.17", features = ["rayon", "approx"] } rayon = "1.10" +kdtree = "0.8.1" +nalgebra = "0.35.0" +scirs2-special = "0.6.5" +serde = { version = "1.0.229", features = ["derive"] } +serde_json = "1.0.151" [dev-dependencies] criterion = { version = "0.8", features = ["html_reports"] } diff --git a/src/covmodel.rs b/src/covmodel.rs new file mode 100644 index 0000000..d4a6bfc --- /dev/null +++ b/src/covmodel.rs @@ -0,0 +1,534 @@ +//! Runtime covariance models. +//! +//! A [`CovModel`] is *built* from a [`CovModelSpec`] via [`CovModel::build`], +//! which validates the raw parameters and precomputes everything that would +//! otherwise be recomputed on every one of the O(n²) pairwise evaluations +//! needed to fill a kriging LHS/RHS: the rescaled length scale, and for +//! Matérn the polynomial coefficients of the half-integer closed form. +//! +//! Formulas follow gstools' own `CovModel.cor` convention exactly (see +//! `gstools/covmodel/models.py`): `correlation(r) = cor(r / len_rescaled)`, +//! `covariance(r) = var * correlation(r)`, with `len_rescaled = len_scale / +//! rescale` and a model-specific rescale factor. + +use crate::covmodel_spec::CovModelSpec; +use std::fmt; + +#[derive(Debug, Clone, PartialEq)] +pub enum ModelError { + InvalidVar(f64), + InvalidLenScale(f64), + InvalidNugget(f64), + InvalidNu(f64), +} + +impl fmt::Display for ModelError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::InvalidVar(v) => write!(f, "var must be >= 0, got {v}"), + Self::InvalidLenScale(v) => write!(f, "len_scale must be > 0, got {v}"), + Self::InvalidNugget(v) => write!(f, "nugget must be >= 0, got {v}"), + Self::InvalidNu(v) => write!(f, "nu must be > 0, got {v}"), + } + } +} + +impl std::error::Error for ModelError {} + +/// Error building a [`CovModel`] from a JSON-encoded [`CovModelSpec`] via +/// [`CovModel::from_json`]: either the JSON itself is malformed/doesn't +/// match the wire format, or it parses fine but [`CovModel::build`] +/// rejects the parameters it describes. +#[derive(Debug)] +pub enum CovModelJsonError { + Json(serde_json::Error), + Model(ModelError), +} + +impl fmt::Display for CovModelJsonError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::Json(e) => write!(f, "invalid covariance model JSON: {e}"), + Self::Model(e) => write!(f, "{e}"), + } + } +} + +impl std::error::Error for CovModelJsonError {} + +impl From for CovModelJsonError { + fn from(e: serde_json::Error) -> Self { + Self::Json(e) + } +} + +impl From for CovModelJsonError { + fn from(e: ModelError) -> Self { + Self::Model(e) + } +} + +/// Gaussian's standard rescale factor `s = sqrt(pi) / 2`, see +/// `Gaussian.default_rescale` in gstools. +fn gaussian_rescale() -> f64 { + std::f64::consts::PI.sqrt() / 2.0 +} + +#[derive(Debug, Clone, PartialEq)] +pub enum CovModel { + Gaussian { + var: f64, + len_rescaled: f64, + nugget: f64, + }, + Exponential { + var: f64, + len_rescaled: f64, + nugget: f64, + }, + Matern { + var: f64, + len_rescaled: f64, + nugget: f64, + nu: f64, + sqrt_nu: f64, + }, + Spherical { + var: f64, + len_rescaled: f64, + nugget: f64, + }, +} + +/// Matérn normalized correlation function, following gstools' +/// `Matern.cor` (`gstools/covmodel/models.py`) term for term: +/// `cor(h) = 2^(1-nu)/Gamma(nu) * (sqrt(nu) h)^nu * K_nu(sqrt(nu) h)`, +/// evaluated via a single `exp(...)` of the log-transformed prefactor to +/// avoid overflow in `Gamma(nu)` for large `nu`, and via the Gaussian +/// far-field approximation for `nu > 20` where the exact term underflows. +fn matern_cor(nu: f64, sqrt_nu: f64, h: f64) -> f64 { + // for nu > 20 we just use the gaussian model + if nu > 20.0 { + return (-(h / 2.0).powi(2)).exp(); + } + if h <= 0.0 { + return 1.0; + } + // calculate by log-transformation to prevent numerical errors + let x = sqrt_nu * h; + let res = ((1.0 - nu) * std::f64::consts::LN_2 - scirs2_special::loggamma(nu) + nu * x.ln()) + .exp() + * scirs2_special::kv(nu, x); + // if nu >> 1 we get errors for the farfield, there 0 is approached; + // covariance is positive + if res.is_finite() { + res.max(0.0) + } else { + 0.0 + } +} + +impl CovModel { + /// Build a [`CovModel`] straight from the JSON wire format a + /// [`CovModelSpec`] deserializes from -- the single entry point Python + /// callers (via PyO3) go through, so the wire format's parsing and its + /// parameter validation both stay colocated with the model types + /// themselves rather than in the PyO3 binding layer. + pub fn from_json(json: &str) -> Result { + let spec: CovModelSpec = serde_json::from_str(json)?; + Ok(Self::build(spec)?) + } + + pub fn build(spec: CovModelSpec) -> Result { + match spec { + CovModelSpec::Gaussian { + var, + len_scale, + nugget, + } => { + check_common(var, len_scale, nugget)?; + Ok(Self::Gaussian { + var, + len_rescaled: len_scale / gaussian_rescale(), + nugget, + }) + } + CovModelSpec::Exponential { + var, + len_scale, + nugget, + } => { + check_common(var, len_scale, nugget)?; + Ok(Self::Exponential { + var, + len_rescaled: len_scale, // rescale = 1 + nugget, + }) + } + CovModelSpec::Matern { + var, + len_scale, + nugget, + nu, + } => { + check_common(var, len_scale, nugget)?; + if nu.is_nan() || nu <= 0.0 { + return Err(ModelError::InvalidNu(nu)); + } + Ok(Self::Matern { + var, + len_rescaled: len_scale, // rescale = 1 + nugget, + nu, + sqrt_nu: nu.sqrt(), + }) + } + CovModelSpec::Spherical { + var, + len_scale, + nugget, + } => { + check_common(var, len_scale, nugget)?; + Ok(Self::Spherical { + var, + len_rescaled: len_scale, // rescale = 1 + nugget, + }) + } + } + } + + fn var(&self) -> f64 { + match self { + Self::Gaussian { var, .. } + | Self::Exponential { var, .. } + | Self::Matern { var, .. } + | Self::Spherical { var, .. } => *var, + } + } + + fn nugget(&self) -> f64 { + match self { + Self::Gaussian { nugget, .. } + | Self::Exponential { nugget, .. } + | Self::Matern { nugget, .. } + | Self::Spherical { nugget, .. } => *nugget, + } + } + + /// Normalized correlation function, `correlation(0) = 1`. + pub fn correlation(&self, r: f64) -> f64 { + let r = r.abs(); + match self { + Self::Gaussian { len_rescaled, .. } => { + let h = r / len_rescaled; + (-(h * h)).exp() + } + Self::Exponential { len_rescaled, .. } => { + let h = r / len_rescaled; + (-h).exp() + } + Self::Matern { + len_rescaled, + nu, + sqrt_nu, + .. + } => { + let h = r / len_rescaled; + matern_cor(*nu, *sqrt_nu, h) + } + Self::Spherical { len_rescaled, .. } => { + // compact support: the polynomial is exactly 0 at h = 1 + let h = (r / len_rescaled).min(1.0); + 1.0 - 1.5 * h + 0.5 * h * h * h + } + } + } + + /// `covariance(r) = var * correlation(r)`. + pub fn covariance(&self, r: f64) -> f64 { + self.var() * self.correlation(r) + } + + /// Isotropic variogram respecting the nugget at `r = 0` (matches + /// gstools' `vario_nugget`: exactly 0 at the origin, `var + nugget - + /// covariance(r)` elsewhere). + pub fn variogram(&self, r: f64) -> f64 { + if r == 0.0 { + 0.0 + } else { + self.var() + self.nugget() - self.covariance(r) + } + } + + /// Sill: `var + nugget`, i.e. `covariance(0)` including the nugget + /// (matches gstools' `cov_nugget` at `r = 0`). + pub fn sill(&self) -> f64 { + self.var() + self.nugget() + } + + /// `covariance(r)`, except at `r = 0` where the nugget is added back + /// in (matches gstools' `cov_nugget`: `sill()` at the origin, + /// `covariance(r)` everywhere else). Used for the "exact" kriging + /// variant, where a target point coinciding with a conditioning point + /// must reproduce that point's value exactly. + pub fn cov_nugget(&self, r: f64) -> f64 { + if r == 0.0 { + self.sill() + } else { + self.covariance(r) + } + } +} + +fn check_common(var: f64, len_scale: f64, nugget: f64) -> Result<(), ModelError> { + if var < 0.0 { + return Err(ModelError::InvalidVar(var)); + } + if len_scale <= 0.0 { + return Err(ModelError::InvalidLenScale(len_scale)); + } + if nugget < 0.0 { + return Err(ModelError::InvalidNugget(nugget)); + } + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn build(spec: CovModelSpec) -> CovModel { + CovModel::build(spec).expect("valid spec") + } + + #[test] + fn correlation_is_one_at_origin() { + for spec in [ + CovModelSpec::Gaussian { + var: 2.0, + len_scale: 3.0, + nugget: 0.1, + }, + CovModelSpec::Exponential { + var: 2.0, + len_scale: 3.0, + nugget: 0.1, + }, + CovModelSpec::Matern { + var: 2.0, + len_scale: 3.0, + nugget: 0.1, + nu: 2.5, + }, + CovModelSpec::Spherical { + var: 2.0, + len_scale: 3.0, + nugget: 0.1, + }, + ] { + let model = build(spec); + assert!((model.correlation(0.0) - 1.0).abs() < 1e-12); + assert!((model.covariance(0.0) - model.var()).abs() < 1e-12); + } + } + + #[test] + fn matern_half_integer_matches_hand_derived_closed_forms() { + // cor(h) = P_n(x) * exp(-x), x = sqrt(nu) * h, len_scale = 1 so + // len_rescaled = h directly. + type ClosedForm = fn(f64) -> f64; + let cases: [(f64, ClosedForm); 3] = [ + (0.5, |x| (-x).exp()), + (1.5, |x| (1.0 + x) * (-x).exp()), + (2.5, |x| (1.0 + x + x * x / 3.0) * (-x).exp()), + ]; + for (nu, closed_form) in cases { + let model = build(CovModelSpec::Matern { + var: 1.0, + len_scale: 1.0, + nugget: 0.0, + nu, + }); + for h in [0.1, 0.5, 1.0, 2.0, 5.0] { + let x = nu.sqrt() * h; + let expected = closed_form(x); + let got = model.correlation(h); + assert!( + (got - expected).abs() < 1e-10, + "nu={nu} h={h}: got {got}, expected {expected}" + ); + } + } + } + + #[test] + fn matern_matches_gstools_reference_for_arbitrary_nu() { + // Reference values from gstools' own `Matern.cor` (scipy.special + // loggamma/kv), computed independently in Python. `nu=25.0` exceeds + // the nu > 20 cutoff, so it exercises the Gaussian far-field + // fallback branch instead of loggamma/kv. + let cases: [(f64, [(f64, f64); 5]); 4] = [ + ( + 0.3, + [ + (0.1, 0.833946507020263), + (0.5, 0.5818747118092827), + (1.0, 0.40558180689489187), + (2.0, 0.21171026561635017), + (5.0, 0.0350803447153922), + ], + ), + ( + 1.3, + [ + (0.1, 0.9913055665037315), + (0.5, 0.8601699813844405), + (1.0, 0.6366803171828707), + (2.0, 0.2917074599714008), + (5.0, 0.01714790657015862), + ], + ), + ( + 3.7, + [ + (0.1, 0.996583355039885), + (0.5, 0.9196887797390918), + (1.0, 0.727966937380819), + (2.0, 0.3296887012120744), + (5.0, 0.00832512681464761), + ], + ), + ( + 25.0, + [ + (0.1, 0.9975031223974601), + (0.5, 0.9394130628134758), + (1.0, 0.7788007830714049), + (2.0, 0.36787944117144233), + (5.0, 0.0019304541362277093), + ], + ), + ]; + for (nu, points) in cases { + let model = build(CovModelSpec::Matern { + var: 1.0, + len_scale: 1.0, + nugget: 0.0, + nu, + }); + for (h, expected) in points { + let got = model.correlation(h); + assert!( + (got - expected).abs() < 1e-9, + "nu={nu} h={h}: got {got}, expected {expected}" + ); + } + } + } + #[test] + fn spherical_matches_reference() { + // cor(h) = 1 - 1.5 h + 0.5 h^3 for h < 1, else 0; rescale = 1, so + // with len_scale = 2 the support ends at r = 2. + let model = build(CovModelSpec::Spherical { + var: 1.0, + len_scale: 2.0, + nugget: 0.0, + }); + let cases = [ + (0.5, 0.6328125), + (1.0, 0.3125), + (-1.0, 0.3125), + (1.8, 0.0145), + (2.0, 0.0), + (4.0, 0.0), + (100.0, 0.0), + ]; + for (r, expected) in cases { + let got = model.correlation(r); + assert!( + (got - expected).abs() < 1e-12, + "r={r}: got {got}, expected {expected}" + ); + } + } + + #[test] + fn spherical_covariance_and_variogram_reach_sill_at_range() { + let model = build(CovModelSpec::Spherical { + var: 2.0, + len_scale: 3.0, + nugget: 0.5, + }); + assert!((model.covariance(1.5) - 2.0 * 0.3125).abs() < 1e-12); + assert_eq!(model.covariance(3.0), 0.0); + assert_eq!(model.covariance(10.0), 0.0); + assert_eq!(model.variogram(0.0), 0.0); + assert!((model.variogram(10.0) - model.sill()).abs() < 1e-12); + assert!((model.cov_nugget(0.0) - 2.5).abs() < 1e-12); + } + + #[test] + fn spherical_builds_from_json() { + let model = CovModel::from_json( + r#"{"type": "spherical", "var": 1.5, "len_scale": 4.0, "nugget": 0.2}"#, + ) + .expect("valid spherical JSON"); + assert_eq!( + model, + CovModel::Spherical { + var: 1.5, + len_rescaled: 4.0, + nugget: 0.2, + } + ); + } + + #[test] + fn spherical_rejects_invalid_len_scale() { + let err = CovModel::build(CovModelSpec::Spherical { + var: 1.0, + len_scale: 0.0, + nugget: 0.0, + }) + .unwrap_err(); + assert_eq!(err, ModelError::InvalidLenScale(0.0)); + } + + #[test] + fn rejects_non_positive_nu() { + for nu in [0.0, -1.5] { + let err = CovModel::build(CovModelSpec::Matern { + var: 1.0, + len_scale: 1.0, + nugget: 0.0, + nu, + }) + .unwrap_err(); + assert_eq!(err, ModelError::InvalidNu(nu)); + } + } + + #[test] + fn correlation_decays_with_distance() { + let model = build(CovModelSpec::Gaussian { + var: 1.0, + len_scale: 5.0, + nugget: 0.0, + }); + assert!(model.correlation(1.0) > model.correlation(5.0)); + assert!(model.correlation(5.0) > model.correlation(20.0)); + assert!(model.correlation(100.0) < 1e-6); + } + + #[test] + fn variogram_is_zero_at_origin_even_with_nugget() { + let model = build(CovModelSpec::Exponential { + var: 1.0, + len_scale: 2.0, + nugget: 0.3, + }); + assert_eq!(model.variogram(0.0), 0.0); + assert!((model.sill() - 1.3).abs() < 1e-12); + } +} diff --git a/src/covmodel_spec.rs b/src/covmodel_spec.rs new file mode 100644 index 0000000..c1d0d1e --- /dev/null +++ b/src/covmodel_spec.rs @@ -0,0 +1,35 @@ +//! Wire format for covariance models: exactly what comes in over JSON. +//! +//! Mirrors the parameter names gstools itself uses (`var`, `len_scale`, +//! `nugget`, `nu`), so a small `to_dict()` on the Python side would produce +//! compatible JSON directly. No validation happens here — a `CovModelSpec` +//! is just "what the wire said", turned into a usable [`crate::covmodel::CovModel`] +//! via [`crate::covmodel::CovModel::build`]. + +use serde::{Deserialize, Serialize}; + +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case")] +pub enum CovModelSpec { + Gaussian { + var: f64, + len_scale: f64, + nugget: f64, + }, + Exponential { + var: f64, + len_scale: f64, + nugget: f64, + }, + Matern { + var: f64, + len_scale: f64, + nugget: f64, + nu: f64, + }, + Spherical { + var: f64, + len_scale: f64, + nugget: f64, + }, +} diff --git a/src/lib.rs b/src/lib.rs index d2523d3..acff9ee 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -14,6 +14,9 @@ use pyo3::prelude::pymodule; +pub mod covmodel; +pub mod local_krige; +pub mod covmodel_spec; pub mod field; pub mod krige; mod short_vec; @@ -21,6 +24,8 @@ pub mod variogram; #[pymodule] mod gstools_core { + use crate::covmodel::CovModel; + use crate::local_krige::calc_field_krige_local; use crate::field::{summator, summator_fourier, summator_incompr}; use crate::krige::{calculator_field_krige, calculator_field_krige_and_variance}; use crate::variogram::{ @@ -28,6 +33,7 @@ mod gstools_core { variogram_unstructured, }; use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2}; + use pyo3::exceptions::PyValueError; use pyo3::prelude::*; #[pymodule_init] @@ -118,6 +124,47 @@ mod gstools_core { calculator_field_krige(krige_mat, krig_vecs, cond, num_threads).into_pyarray(py) } + #[pyfunction(name = "calc_field_krige_local")] + #[allow(clippy::too_many_arguments)] + fn calc_field_krige_local_py<'py>( + py: Python<'py>, + cond_pos: PyReadonlyArray2, + cond_val: PyReadonlyArray1, + target_pos: PyReadonlyArray2, + cov_model_json: &str, + cond_err: PyReadonlyArray1, + drift_cond: PyReadonlyArray2, + drift_target: PyReadonlyArray2, + unbiased: bool, + exact: bool, + local_radius: f64, + num_threads: Option, + ) -> PyResult<(Bound<'py, PyArray1>, Bound<'py, PyArray1>)> { + let model = CovModel::from_json(cov_model_json) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + let cond_pos = cond_pos.as_array(); + let cond_val = cond_val.as_array(); + let target_pos = target_pos.as_array(); + let cond_err = cond_err.as_array(); + let drift_cond = drift_cond.as_array(); + let drift_target = drift_target.as_array(); + let (field, error) = calc_field_krige_local( + cond_pos, + cond_val, + target_pos, + &model, + cond_err, + drift_cond, + drift_target, + unbiased, + exact, + local_radius, + num_threads, + ) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + Ok((field.into_pyarray(py), error.into_pyarray(py))) + } + #[pyfunction(name = "variogram_structured")] fn variogram_structured_py<'py>( py: Python<'py>, diff --git a/src/local_krige.rs b/src/local_krige.rs new file mode 100644 index 0000000..737fb12 --- /dev/null +++ b/src/local_krige.rs @@ -0,0 +1,702 @@ +//! Local (moving-neighborhood) kriging. +//! +//! Rust port of `calc_field_krige_local` from GSTools' +//! `dev/local_kriging/local_kriging_numpy.py`. That Python function is +//! itself already written as a stand-in for the planned +//! `gstools_core::krige::calc_field_krige_local` -- see its docstring: it +//! only ever receives plain arrays and scalars, no Python callables, no +//! `Krige`/`CovModel` objects with behaviour beyond `covariance`/ +//! `cov_nugget`. Anisotropy/rotation stay on the Python side: positions +//! here are assumed already isometrized (isotropic), so every distance is +//! a plain Euclidean norm. +//! +//! The per-target-point loop in [`calc_field_krige_local`] is deliberately +//! *not* parallelized yet: it is written as a `.map()` over an iterator of +//! independent, side-effect-free target-point computations (each only +//! borrows shared, immutable data), so switching to +//! `rayon`'s `.into_par_iter()` later -- as `krige.rs` does for the global +//! case -- is a one-line change, not a restructuring. +//! +//! Unlike global kriging (see `krige.rs`), local kriging needs an actual +//! dense solve per target point (each neighborhood is a different small +//! system) -- there is no pre-inverted matrix to just contract against. +//! [`solve_dense`] uses `nalgebra`'s partial-pivoted `LU` for the common +//! (well-conditioned) case, falling back to an `SVD`-based least-squares +//! solve -- analogous to NumPy's `lstsq` in the Python prototype -- for +//! singular systems. + +use std::fmt; + +use kdtree::{distance::squared_euclidean, KdTree}; +use nalgebra::{DMatrix, DVector}; +use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis}; +use rayon::prelude::*; + +use crate::covmodel::CovModel; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum LocalKrigeError { + /// No conditioning point fell within `local_radius` of this target + /// point (mirrors the `ValueError` raised in the Python prototype). + NoNeighbors { target_index: usize }, + /// The local kriging matrix was singular even for [`solve_dense`]'s + /// least-squares fallback -- in practice essentially unreachable, see + /// its docs. + SingularSystem { target_index: usize }, +} + +impl fmt::Display for LocalKrigeError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::NoNeighbors { target_index } => write!( + f, + "no conditioning points within local_radius for target point {target_index}; \ + increase local_radius" + ), + Self::SingularSystem { target_index } => write!( + f, + "local kriging matrix for target point {target_index} is singular even for the \ + least-squares fallback" + ), + } + } +} + +impl std::error::Error for LocalKrigeError {} + +/// Assemble the local kriging matrix (LHS) for one target point's +/// neighborhood: pairwise covariances between the `k` neighbors, plus +/// measurement error on the diagonal, plus (optionally) the +/// unbiasedness row/column and drift rows/columns. Mirrors the `lhs` +/// block in `local_kriging_numpy.calc_field_krige_local`. +/// +/// # Arguments +/// * `nbr_pos` - neighbor positions, shape `(dim, k)` +/// * `model` - covariance model +/// * `cond_err` - measurement error per neighbor, shape `(k,)` +/// * `unbiased` - whether to add the Ordinary-Kriging Lagrange row/column +/// * `drift_nbrs` - drift terms evaluated at the neighbor positions, +/// shape `(p, k)` (`p` may be 0) +/// +/// # Returns +/// the `(k + p + unbiased) x (k + p + unbiased)` LHS matrix +pub fn build_local_lhs( + nbr_pos: ArrayView2<'_, f64>, + model: &CovModel, + cond_err: ArrayView1<'_, f64>, + unbiased: bool, + drift_nbrs: ArrayView2<'_, f64>, +) -> Array2 { + let k = nbr_pos.ncols(); + let drift_no = drift_nbrs.nrows(); + let pad = drift_no + unbiased as usize; + let size = k + pad; + + let mut lhs = Array2::::zeros((size, size)); + for i in 0..k { + let pi = nbr_pos.column(i); + for j in 0..k { + let pj = nbr_pos.column(j); + let dist = euclidean_dist(pi, pj); + lhs[[i, j]] = model.covariance(dist); + } + lhs[[i, i]] += cond_err[i]; + } + if unbiased { + let u = k; + for i in 0..k { + lhs[[u, i]] = 1.0; + lhs[[i, u]] = 1.0; + } + } + if drift_no > 0 { + let offset = k + unbiased as usize; + for d in 0..drift_no { + for i in 0..k { + lhs[[offset + d, i]] = drift_nbrs[[d, i]]; + lhs[[i, offset + d]] = drift_nbrs[[d, i]]; + } + } + } + lhs +} + +/// Assemble the local kriging vector (RHS) for one target point: the +/// covariance (or, in `exact` mode, the nugget-aware covariance) between +/// each neighbor and the target point, plus the unbiasedness/drift +/// entries matching [`build_local_lhs`]'s padding. Mirrors the `rhs` +/// block in `local_kriging_numpy.calc_field_krige_local`. +/// +/// # Arguments +/// * `nbr_dists` - distance from each neighbor to the target point, +/// shape `(k,)` +/// * `model` - covariance model +/// * `exact` - use `model.cov_nugget` instead of `model.covariance`, so a +/// target point that coincides with a neighbor reproduces its value +/// exactly +/// * `unbiased` - must match [`build_local_lhs`] +/// * `drift_target` - drift terms evaluated at the target position, +/// shape `(p,)`, `p` must match [`build_local_lhs`]'s `drift_nbrs` +/// +/// # Returns +/// the RHS vector, same length as [`build_local_lhs`]'s matrix dimension +pub fn build_local_rhs( + nbr_dists: &[f64], + model: &CovModel, + exact: bool, + unbiased: bool, + drift_target: ArrayView1<'_, f64>, +) -> Array1 { + let k = nbr_dists.len(); + let drift_no = drift_target.len(); + let pad = drift_no + unbiased as usize; + + let mut rhs = Array1::::zeros(k + pad); + for (i, &dist) in nbr_dists.iter().enumerate() { + rhs[i] = if exact { + model.cov_nugget(dist) + } else { + model.covariance(dist) + }; + } + if unbiased { + rhs[k] = 1.0; + } + if drift_no > 0 { + let offset = k + unbiased as usize; + for d in 0..drift_no { + rhs[offset + d] = drift_target[d]; + } + } + rhs +} + +/// Tolerance below which an SVD singular value is treated as zero in the +/// least-squares fallback -- analogous to NumPy `lstsq`'s default `rcond` +/// cutoff. +const LSTSQ_EPS: f64 = 1e-12; + +/// Solve the dense linear system `a @ x = b`. Pure function (takes views, +/// returns a fresh `Array1`), so it is independently testable with plain +/// literal matrices, unrelated to kriging. +/// +/// Local kriging systems are small (a handful to a few hundred neighbors), +/// so this is not a performance bottleneck compared to the covariance +/// evaluations that build `a`. +/// +/// Tries `nalgebra::linalg::LU` (partial-pivoted Gaussian elimination) +/// first -- cheap, and sufficient for the well-conditioned case (nugget > +/// 0 keeps the diagonal away from zero). `LU::solve` only treats an +/// *exactly* zero pivot as singular (no tolerance), so it can, in +/// principle, still miss a numerically near-singular system; that +/// tradeoff is accepted here in favor of the fallback below actually +/// covering the case NumPy's `lstsq` was covering in the Python +/// prototype: if `LU` reports singular, this falls back to +/// `nalgebra::linalg::SVD::solve`, i.e. a least-squares solve via the +/// Moore-Penrose pseudo-inverse (singular values below [`LSTSQ_EPS`] are +/// treated as zero). That fallback finds *a* solution for consistent +/// rank-deficient systems (e.g. two neighbors at the same position) and +/// the residual-minimizing solution for inconsistent ones, rather than +/// erroring out -- so `Err(SingularSystem)` is now reserved for cases +/// where even that fails (in practice: essentially never, since both `u` +/// and `v` are requested from the SVD). +pub fn solve_dense( + a: ArrayView2<'_, f64>, + b: ArrayView1<'_, f64>, +) -> Result, SingularSystem> { + let n = b.len(); + assert_eq!(a.nrows(), n, "matrix/vector size mismatch"); + assert_eq!(a.ncols(), n, "matrix must be square"); + + // Convert `ndarray` views to `nalgebra` types, required by + // `nalgebra::linalg::LU`/`SVD`. + let a_na = DMatrix::from_fn(n, n, |r, c| a[[r, c]]); + let b_na = DVector::from_fn(n, |r, _| b[r]); + + // try cheap LU first, then fall back to SVD least-squares if it reports singular + if let Some(x) = a_na.clone().lu().solve(&b_na) { + return Ok(Array1::from_iter(x.iter().copied())); + } + + a_na.svd(true, true) + .solve(&b_na, LSTSQ_EPS) + .map(|x| Array1::from_iter(x.iter().copied())) + .map_err(|_| SingularSystem) +} + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct SingularSystem; + +fn euclidean_dist(a: ArrayView1<'_, f64>, b: ArrayView1<'_, f64>) -> f64 { + a.iter() + .zip(b.iter()) // (&f64, &f64) + .map(|(x, y)| (x - y) * (x - y)) // |(x, y)| is type for closure + .sum::() // sum over all dimensions + .sqrt() +} + +/// Evaluate local (moving-neighborhood) kriging at every target point. +/// +/// Rust port of `local_kriging_numpy.calc_field_krige_local`. Builds a +/// k-d tree over the conditioning positions once, then for each target +/// point: finds every conditioning point within `local_radius`, +/// assembles that neighborhood's local system via [`build_local_lhs`]/ +/// [`build_local_rhs`], solves it with [`solve_dense`], and accumulates +/// the field value and (raw, not yet `sill - error`) error variance. +/// +/// # Arguments +/// * `cond_pos` - isometrized conditioning positions, shape `(dim, n)` +/// * `cond_val` - centered conditioning values, shape `(n,)` +/// * `target_pos` - isometrized target positions, shape `(dim, m)` +/// * `model` - covariance model +/// * `cond_err` - measurement error per conditioning point, shape `(n,)` +/// * `drift_cond` - drift terms at conditioning positions, shape `(p, n)` +/// (`p` may be 0) +/// * `drift_target` - drift terms at target positions, shape `(p, m)` +/// * `unbiased` - whether to enforce the Ordinary-Kriging unbiasedness +/// condition +/// * `exact` - whether to use `model.cov_nugget` on the RHS diagonal case +/// (target coincides with a conditioning point) +/// * `local_radius` - search radius (isometrized distance); every +/// conditioning point within this distance of a target point enters +/// its local system -- the only neighborhood-selection criterion +/// * `num_threads` - reserved for a future `rayon` parallelization of the +/// per-target-point loop below; unused for now (see module docs) +/// +/// # Returns +/// `(field, error)`, both shape `(m,)`; `field` is the raw local-kriging +/// estimate, `error` is the raw `rhs . weights` (not yet post-processed +/// into a variance) +// 11 parameters, matching the Python signature it mirrors 1:1 (see module +// docs) -- not a case of a function that grew organically and should be +// split up. +#[allow(clippy::too_many_arguments)] +pub fn calc_field_krige_local( + cond_pos: ArrayView2<'_, f64>, + cond_val: ArrayView1<'_, f64>, + target_pos: ArrayView2<'_, f64>, + model: &CovModel, + cond_err: ArrayView1<'_, f64>, + drift_cond: ArrayView2<'_, f64>, + drift_target: ArrayView2<'_, f64>, + unbiased: bool, + exact: bool, + local_radius: f64, + num_threads: Option, +) -> Result<(Array1, Array1), LocalKrigeError> { + // Unused until the loop below becomes a `par_iter` (see module docs). + // let _ = num_threads; + + let dim = cond_pos.nrows(); + let cond_no = cond_pos.ncols(); + let pnt_cnt = target_pos.ncols(); + + // Create Kdtree for fast neighbor search. Each point is a column of `cond_pos`, + // need to transpose it to get a Vec for each point. + let mut tree: KdTree> = KdTree::new(dim); + for i in 0..cond_no { + tree.add(cond_pos.column(i).to_vec(), i) + .expect("cond_pos columns are finite and match `dim`"); + } + // Using squared_euclidean distance function with squared distance + // to avoid computing square roots unnecessarily. + // Kdtree crate docs recommend this + let radius_sq = local_radius * local_radius; + + // Closure + // unlike functions closures can access variables from the surrounding scope + // this makes parallelization easier and faster because all threads can access the same data without having to copy it + // j is the index of the target point in the target_pos array + let solve_one = |j: usize| -> Result<(f64, f64), LocalKrigeError> { + let target_point = target_pos.column(j).to_vec(); + + // tree.within() returns a Vec of tuples (squared_distance, index) of all points within the radius + let neighbors = tree + .within(&target_point, radius_sq, &squared_euclidean) + .expect("target point is finite and matches `dim`"); + if neighbors.is_empty() { + return Err(LocalKrigeError::NoNeighbors { target_index: j }); + } + + let nbr_indices: Vec = neighbors.iter().map(|&(_, idx)| *idx).collect(); + let nbr_dists: Vec = neighbors + .iter() + .map(|&(sq_dist, _)| sq_dist.sqrt()) + .collect(); + let k = nbr_indices.len(); + + // Select the neighbor data from the full conditioning arrays using the indices found by the k-d tree search + let nbr_pos: Array2 = cond_pos.select(Axis(1), &nbr_indices); + let nbr_cond_err: Array1 = cond_err.select(Axis(0), &nbr_indices); + let nbr_drift: Array2 = drift_cond.select(Axis(1), &nbr_indices); + let nbr_cond_val: Array1 = cond_val.select(Axis(0), &nbr_indices); + + let lhs: Array2 = build_local_lhs( + nbr_pos.view(), + model, + nbr_cond_err.view(), + unbiased, + nbr_drift.view(), + ); + let rhs: Array1 = + build_local_rhs(&nbr_dists, model, exact, unbiased, drift_target.column(j)); + + let weights: Array1 = solve_dense(lhs.view(), rhs.view()) + .map_err(|_| LocalKrigeError::SingularSystem { target_index: j })?; + + let pad = drift_target.nrows() + unbiased as usize; + let mut local_cond: Array1 = Array1::::zeros(k + pad); + for i in 0..k { + local_cond[i] = nbr_cond_val[i]; + } + + let field_j: f64 = local_cond.dot(&weights); + let error_j: f64 = rhs.dot(&weights); + Ok((field_j, error_j)) + }; + + // let (field, error): (Vec, Vec) = (0..pnt_cnt) // range over target points is the parameter for the closure solve_one + // .map(solve_one) + // .collect::, _>>()? + // .into_iter() + // .unzip(); + + let results: Vec<(f64, f64)> = rayon::ThreadPoolBuilder::new() + .num_threads(num_threads.unwrap_or(rayon::current_num_threads())) + .build() + .unwrap() + .install(|| { + (0..pnt_cnt) + .into_par_iter() + .map(solve_one) + .collect::, _>>() + })?; + + let (field, error): (Vec, Vec) = results.into_iter().unzip(); + + Ok((Array1::from(field), Array1::from(error))) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::{array, Array2}; + + fn gaussian(var: f64, len_scale: f64, nugget: f64) -> CovModel { + CovModel::build(crate::covmodel_spec::CovModelSpec::Gaussian { + var, + len_scale, + nugget, + }) + .unwrap() + } + + // --- solve_dense ----------------------------------------------------- + + #[test] + fn solve_dense_matches_hand_solved_system() { + // [2 1; 1 3] x = [3; 5] -> x = [4/5, 7/5] + let a = array![[2.0, 1.0], [1.0, 3.0]]; + let b = array![3.0, 5.0]; + let x = solve_dense(a.view(), b.view()).unwrap(); + assert!((x[0] - 0.8).abs() < 1e-10); + assert!((x[1] - 1.4).abs() < 1e-10); + } + + #[test] + fn solve_dense_needs_pivoting() { + // zero in the natural pivot position -> must swap rows to solve + let a = array![[0.0, 1.0], [1.0, 1.0]]; + let b = array![2.0, 3.0]; + let x = solve_dense(a.view(), b.view()).unwrap(); + assert!((x[0] - 1.0).abs() < 1e-10); + assert!((x[1] - 2.0).abs() < 1e-10); + } + + #[test] + fn solve_dense_lstsq_fallback_finds_min_norm_solution_for_consistent_singular_system() { + // row2 = 2 * row1, and b is consistent with that (b2 = 2*b1): the + // system x1 + 2*x2 = 1 has infinitely many exact solutions. LU's + // pivot search hits an exact zero (row2 becomes all-zero after + // eliminating row1) and reports singular; the SVD fallback must + // not error out here -- it picks the minimum-norm one among all + // exact solutions, which (for a single constraint a.x = c) is + // x = c * a / |a|^2 = (1, 2) / 5. + let a = array![[1.0, 2.0], [2.0, 4.0]]; + let b = array![1.0, 2.0]; + + let x = solve_dense(a.view(), b.view()).expect("lstsq fallback should not error out"); + + assert!((x[0] - 0.2).abs() < 1e-8); + assert!((x[1] - 0.4).abs() < 1e-8); + // it is a genuine solution of the (consistent) original system: + assert!((a.row(0).dot(&x) - b[0]).abs() < 1e-8); + assert!((a.row(1).dot(&x) - b[1]).abs() < 1e-8); + } + + #[test] + fn solve_dense_lstsq_fallback_minimizes_residual_for_inconsistent_singular_system() { + // Same singular `a`, but now b is *inconsistent* with row2 = 2 * + // row1 (b2 != 2*b1): no exact solution exists. This is the case + // NumPy's `lstsq` (and our SVD fallback) is actually for -- + // minimizing |a.x - b|^2. Since a = (1,2)(1,2)^T, a.x only depends + // on s = x1 + 2*x2, and minimizing |(s,2s) - (1,3)|^2 over s gives + // s = 1.4 (hand-derived below), with the minimum-norm x for that s + // being s * (1, 2) / 5 = (0.28, 0.56). + let a = array![[1.0, 2.0], [2.0, 4.0]]; + let b = array![1.0, 3.0]; + + let x = solve_dense(a.view(), b.view()).expect("lstsq fallback should not error out"); + + assert!((x[0] - 0.28).abs() < 1e-8); + assert!((x[1] - 0.56).abs() < 1e-8); + + // sanity: this residual must be strictly smaller than for a + // (non-least-squares) candidate exact solution of just row 1. + let exact_row0_only = array![1.0, 0.0]; // satisfies row 0 exactly, ignores row 1 + let residual = |v: &Array1| { + let r0 = a.row(0).dot(v) - b[0]; + let r1 = a.row(1).dot(v) - b[1]; + r0 * r0 + r1 * r1 + }; + assert!(residual(&x) < residual(&exact_row0_only)); + } + + // --- build_local_lhs / build_local_rhs -------------------------------- + + #[test] + fn lhs_matches_model_covariance_and_adds_cond_err_on_diagonal() { + let model = gaussian(1.5, 4.0, 0.1); + let nbr_pos = array![[0.0, 3.0], [0.0, 0.0]]; // 1D positions on the x-axis + let cond_err = array![0.1, 0.2]; + let drift = Array2::::zeros((0, 2)); + + let lhs = build_local_lhs(nbr_pos.view(), &model, cond_err.view(), false, drift.view()); + + assert_eq!(lhs.dim(), (2, 2)); + assert!((lhs[[0, 0]] - (model.covariance(0.0) + 0.1)).abs() < 1e-12); + assert!((lhs[[1, 1]] - (model.covariance(0.0) + 0.2)).abs() < 1e-12); + assert!((lhs[[0, 1]] - model.covariance(3.0)).abs() < 1e-12); + assert_eq!(lhs[[0, 1]], lhs[[1, 0]]); // symmetric + } + + #[test] + fn lhs_unbiased_adds_ones_row_and_column() { + let model = gaussian(1.0, 2.0, 0.0); + let nbr_pos = array![[0.0, 1.0]]; + let cond_err = array![0.0, 0.0]; + let drift = Array2::::zeros((0, 2)); + + let lhs = build_local_lhs(nbr_pos.view(), &model, cond_err.view(), true, drift.view()); + + assert_eq!(lhs.dim(), (3, 3)); + assert_eq!(lhs.row(2).to_vec(), vec![1.0, 1.0, 0.0]); + assert_eq!(lhs.column(2).to_vec(), vec![1.0, 1.0, 0.0]); + } + + #[test] + fn rhs_uses_cov_nugget_only_in_exact_mode() { + let model = gaussian(1.5, 4.0, 0.2); + let drift = Array1::::zeros(0); + + let rhs_normal = build_local_rhs(&[0.0, 3.0], &model, false, false, drift.view()); + let rhs_exact = build_local_rhs(&[0.0, 3.0], &model, true, false, drift.view()); + + // away from the origin both agree... + assert!((rhs_normal[1] - rhs_exact[1]).abs() < 1e-12); + // ...but at distance 0 `exact` must include the nugget (the sill). + assert!((rhs_normal[0] - model.covariance(0.0)).abs() < 1e-12); + assert!((rhs_exact[0] - model.sill()).abs() < 1e-12); + } + + // --- calc_field_krige_local ------------------------------------------- + + fn scattered_cond() -> (Array2, Array1) { + let cond_pos = array![[0.0, 4.0, 8.0, 2.0, 6.0], [0.0, 3.0, 1.0, 7.0, 5.0]]; + let cond_val = array![0.3, -0.5, 1.2, 0.1, -0.8]; + (cond_pos, cond_val) + } + + #[test] + fn exact_mode_reproduces_conditioning_value_at_conditioning_point() { + let model = gaussian(2.0, 3.0, 0.0); + let (cond_pos, cond_val) = scattered_cond(); + let cond_no = cond_val.len(); + let cond_err = Array1::::zeros(cond_no); + let drift_cond = Array2::::zeros((0, cond_no)); + + // target == the first conditioning point exactly + let target_pos = cond_pos.column(0).to_owned().insert_axis(Axis(1)); + let drift_target = Array2::::zeros((0, 1)); + + let (field, _error) = calc_field_krige_local( + cond_pos.view(), + cond_val.view(), + target_pos.view(), + &model, + cond_err.view(), + drift_cond.view(), + drift_target.view(), + true, + true, + 100.0, + None, + ) + .unwrap(); + + assert!((field[0] - cond_val[0]).abs() < 1e-6); + } + + #[test] + fn no_neighbors_within_radius_is_an_error() { + let model = gaussian(1.0, 2.0, 0.1); + let (cond_pos, cond_val) = scattered_cond(); + let cond_no = cond_val.len(); + let cond_err = Array1::::from_elem(cond_no, 0.1); + let drift_cond = Array2::::zeros((0, cond_no)); + + let target_pos = array![[1000.0], [1000.0]]; // far from every cond point + let drift_target = Array2::::zeros((0, 1)); + + let result = calc_field_krige_local( + cond_pos.view(), + cond_val.view(), + target_pos.view(), + &model, + cond_err.view(), + drift_cond.view(), + drift_target.view(), + true, + false, + 1.0, + None, + ); + + assert_eq!( + result, + Err(LocalKrigeError::NoNeighbors { target_index: 0 }) + ); + } + + #[test] + fn matches_manual_lhs_rhs_solve_when_radius_covers_everyone() { + let model = gaussian(1.5, 3.0, 0.1); + let (cond_pos, cond_val) = scattered_cond(); + let cond_no = cond_val.len(); + let cond_err = Array1::::from_elem(cond_no, 0.1); + let drift_cond = Array2::::zeros((0, cond_no)); + + let target_pos = array![[3.0], [3.0]]; + let drift_target = Array2::::zeros((0, 1)); + + let (field, error) = calc_field_krige_local( + cond_pos.view(), + cond_val.view(), + target_pos.view(), + &model, + cond_err.view(), + drift_cond.view(), + drift_target.view(), + true, + false, + 1e6, + None, + ) + .unwrap(); + + // Manually assemble the very same (full-neighborhood) system with + // build_local_lhs/build_local_rhs/solve_dense directly, bypassing + // the k-d tree entirely, and check the two paths agree. + let target = target_pos.column(0); + let dists: Vec = (0..cond_no) + .map(|i| euclidean_dist(cond_pos.column(i), target)) + .collect(); + let lhs = build_local_lhs( + cond_pos.view(), + &model, + cond_err.view(), + true, + drift_cond.view(), + ); + let rhs = build_local_rhs(&dists, &model, false, true, drift_target.column(0)); + let weights = solve_dense(lhs.view(), rhs.view()).unwrap(); + + let mut local_cond = Array1::::zeros(cond_no + 1); + for i in 0..cond_no { + local_cond[i] = cond_val[i]; + } + let expected_field = local_cond.dot(&weights); + let expected_error = rhs.dot(&weights); + + assert!((field[0] - expected_field).abs() < 1e-8); + assert!((error[0] - expected_error).abs() < 1e-8); + } + + // `cond_pos`/`target_pos` being `ArrayView2` only fixes the array's + // *rank* (a matrix: axis 0 = spatial dimension, axis 1 = point index). + // The actual spatial dimensionality is `cond_pos.nrows()`, read at + // runtime -- so the exact same function works unchanged for 1D and 3D + // positions, not just the 2D data used in every other test here. + + #[test] + fn works_with_1d_positions() { + let model = gaussian(1.0, 3.0, 0.1); + let cond_pos = array![[0.0, 2.0, 5.0, 8.0]]; // shape (1, 4) + let cond_val = array![0.1, 0.4, -0.2, 0.9]; + let cond_no = cond_val.len(); + let cond_err = Array1::::from_elem(cond_no, 0.1); + let drift_cond = Array2::::zeros((0, cond_no)); + let target_pos = array![[1.0, 6.0]]; // shape (1, 2) + let drift_target = Array2::::zeros((0, 2)); + + let (field, _error) = calc_field_krige_local( + cond_pos.view(), + cond_val.view(), + target_pos.view(), + &model, + cond_err.view(), + drift_cond.view(), + drift_target.view(), + true, + false, + 100.0, + None, + ) + .unwrap(); + + assert_eq!(field.len(), 2); + } + + #[test] + fn works_with_3d_positions() { + let model = gaussian(1.0, 3.0, 0.1); + let cond_pos = array![ + [0.0, 2.0, 5.0, 8.0], + [1.0, 0.0, 3.0, 2.0], + [4.0, 1.0, 0.0, 5.0], + ]; // shape (3, 4) + let cond_val = array![0.1, 0.4, -0.2, 0.9]; + let cond_no = cond_val.len(); + let cond_err = Array1::::from_elem(cond_no, 0.1); + let drift_cond = Array2::::zeros((0, cond_no)); + let target_pos = array![[1.0, 6.0], [0.5, 1.0], [2.0, 3.0]]; // shape (3, 2) + let drift_target = Array2::::zeros((0, 2)); + + let (field, _error) = calc_field_krige_local( + cond_pos.view(), + cond_val.view(), + target_pos.view(), + &model, + cond_err.view(), + drift_cond.view(), + drift_target.view(), + true, + false, + 100.0, + None, + ) + .unwrap(); + + assert_eq!(field.len(), 2); + } +}