diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 4276037..35c3d72 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -25,6 +25,8 @@ jobs: toolchain: stable override: true components: rustfmt, clippy + - name: Install Dependencies + run: sudo apt-get install -y build-essentials gfortran - uses: actions-rs/cargo@v1 name: fmt with: @@ -50,7 +52,7 @@ jobs: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: '3.10' + python-version: "3.10" architecture: x64 - name: Build Wheels uses: messense/maturin-action@v1 @@ -78,7 +80,7 @@ jobs: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: '3.10' + python-version: "3.10" architecture: ${{ matrix.target }} - name: Install Rust Toolchain uses: actions-rs/toolchain@v1 @@ -110,7 +112,7 @@ jobs: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: '3.10' + python-version: "3.10" - name: Install Rust toolchain uses: actions-rs/toolchain@v1 with: @@ -152,7 +154,7 @@ jobs: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: '3.10' + python-version: "3.10" architecture: x64 - name: Install dependencies run: | diff --git a/.gitignore b/.gitignore index 8130c3a..2e63061 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ debug/ target/ +*.png diff --git a/Cargo.lock b/Cargo.lock index f5308ef..90f86e2 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,12 @@ # It is not intended for manual editing. version = 3 +[[package]] +name = "adler" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f26201604c87b1e01bd3d98f8d5d9a8fcbb815e8cedb41ffccbeb4bf593a35fe" + [[package]] name = "aho-corasick" version = "1.1.3" @@ -11,6 +17,21 @@ dependencies = [ "memchr", ] +[[package]] +name = "android-tzdata" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e999941b234f3131b00bc13c22d06e8c5ff726d1b6318ac7eb276997bbb4fef0" + +[[package]] +name = "android_system_properties" +version = "0.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "819e7219dbd41043ac279b19830f2efc897156490d7fd6ea916720117ee66311" +dependencies = [ + "libc", +] + [[package]] name = "anes" version = "0.1.6" @@ -23,6 +44,21 @@ version = "1.0.7" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "038dfcf04a5feb68e9c60b21c9625a54c2c0616e79b72b0fd87075a056ae1d1b" +[[package]] +name = "anyhow" +version = "1.0.86" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b3d1d046238990b9cf5bcde22a3fb3584ee5cf65fb2765f454ed428c7a0063da" + +[[package]] +name = "approx" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f2a05fd1bd10b2527e20a2cd32d8873d115b8b39fe219ee25f42a8aca6ba278" +dependencies = [ + "num-traits", +] + [[package]] name = "approx" version = "0.5.1" @@ -38,6 +74,18 @@ version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" +[[package]] +name = "base64" +version = "0.22.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + [[package]] name = "bitflags" version = "2.6.0" @@ -50,18 +98,71 @@ version = "3.16.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "79296716171880943b8470b5f8d03aa55eb2e645a4874bdbb28adb49162e012c" +[[package]] +name = "bytemuck" +version = "1.14.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2ef034f05691a48569bd920a96c81b9d91bbad1ab5ac7c4616c1f6ef36cb79f" + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + [[package]] name = "cast" version = "0.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" +[[package]] +name = "cauchy" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ff11ddd2af3b5e80dd0297fee6e56ac038d9bdc549573cdb51bd6d2efe7f05e" +dependencies = [ + "num-complex", + "num-traits", + "rand", + "serde", +] + +[[package]] +name = "cblas-sys" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6feecd82cce51b0204cf063f0041d69f24ce83f680d87514b004248e7b0fa65" +dependencies = [ + "libc", +] + +[[package]] +name = "cc" +version = "1.0.104" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "74b6a57f98764a267ff415d50a25e6e166f3831a5071af4995296ea97d210490" + [[package]] name = "cfg-if" version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" +[[package]] +name = "chrono" +version = "0.4.38" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a21f936df1771bf62b77f047b726c4625ff2e8aa607c01ec06e5a05bd8463401" +dependencies = [ + "android-tzdata", + "iana-time-zone", + "js-sys", + "num-traits", + "wasm-bindgen", + "windows-targets", +] + [[package]] name = "ciborium" version = "0.2.2" @@ -114,6 +215,73 @@ version = "0.7.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4b82cf0babdbd58558212896d1a4272303a57bdb245c2bf1147185fb45640e70" +[[package]] +name = "color_quant" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d7b894f5411737b7867f4827955924d7c254fc9f4d91a6aad6b097804b1018b" + +[[package]] +name = "core-foundation" +version = "0.9.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91e195e091a93c46f7102ec7818a2aa394e1e1771c3ab4825963fa03e45afb8f" +dependencies = [ + "core-foundation-sys", + "libc", +] + +[[package]] +name = "core-foundation-sys" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "06ea2b9bc92be3c2baa9334a323ebca2d6f074ff852cd1d7b11064035cd3868f" + +[[package]] +name = "core-graphics" +version = "0.23.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "970a29baf4110c26fedbc7f82107d42c23f7e88e404c4577ed73fe99ff85a212" +dependencies = [ + "bitflags 1.3.2", + "core-foundation", + "core-graphics-types", + "foreign-types 0.5.0", + "libc", +] + +[[package]] +name = "core-graphics-types" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "45390e6114f68f718cc7a830514a96f903cccd70d02a8f6d9f643ac4ba45afaf" +dependencies = [ + "bitflags 1.3.2", + "core-foundation", + "libc", +] + +[[package]] +name = "core-text" +version = "20.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c9d2790b5c08465d49f8dc05c8bcae9fea467855947db39b0f8145c091aaced5" +dependencies = [ + "core-foundation", + "core-graphics", + "foreign-types 0.5.0", + "libc", +] + +[[package]] +name = "crc32fast" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b3855a8a784b474f333699ef2bbca9db2c4a1f6d9088a90a2d25b1eb53111eaa" +dependencies = [ + "cfg-if", +] + [[package]] name = "criterion" version = "0.5.1" @@ -181,12 +349,233 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7a81dae078cea95a014a339291cec439d2f232ebe854a9d672b796c6afafa9b7" +[[package]] +name = "cstr" +version = "0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68523903c8ae5aacfa32a0d9ae60cadeb764e1da14ee0d26b1f3089f13a54636" +dependencies = [ + "proc-macro2", + "quote", +] + +[[package]] +name = "dirs" +version = "3.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "30baa043103c9d0c2a57cf537cc2f35623889dc0d405e6c3cccfadbc81c71309" +dependencies = [ + "dirs-sys", +] + +[[package]] +name = "dirs-next" +version = "2.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b98cf8ebf19c3d1b223e151f99a4f9f0690dca41414773390fc824184ac833e1" +dependencies = [ + "cfg-if", + "dirs-sys-next", +] + +[[package]] +name = "dirs-sys" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1b1d1d91c932ef41c0f2663aa8b0ca0342d444d842c06914aa0a7e352d0bada6" +dependencies = [ + "libc", + "redox_users", + "winapi", +] + +[[package]] +name = "dirs-sys-next" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ebda144c4fe02d1f7ea1a7d9641b6fc6b580adcfa024ae48797ecdeb6825b4d" +dependencies = [ + "libc", + "redox_users", + "winapi", +] + +[[package]] +name = "dlib" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "330c60081dcc4c72131f8eb70510f1ac07223e5d4163db481a04a0befcffa412" +dependencies = [ + "libloading", +] + +[[package]] +name = "dwrote" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "439a1c2ba5611ad3ed731280541d36d2e9c4ac5e7fb818a27b604bdc5a6aa65b" +dependencies = [ + "lazy_static", + "libc", + "winapi", + "wio", +] + [[package]] name = "either" version = "1.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" +[[package]] +name = "errno" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "534c5cf6194dfab3db3242765c03bbe257cf92f22b38f6bc0c58d59108a820ba" +dependencies = [ + "libc", + "windows-sys", +] + +[[package]] +name = "fast_polynomial" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62eea6ee590b08a5f8b1139f4d6caee195b646d0c07e4b1808fbd5c4dea4829a" +dependencies = [ + "num-traits", +] + +[[package]] +name = "fastrand" +version = "2.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9fc0510504f03c51ada170672ac806f1f105a88aa97a5281117e1ddc3368e51a" + +[[package]] +name = "fdeflate" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4f9bfee30e4dedf0ab8b422f03af778d9612b63f502710fc500a334ebe2de645" +dependencies = [ + "simd-adler32", +] + +[[package]] +name = "filetime" +version = "0.2.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ee447700ac8aa0b2f2bd7bc4462ad686ba06baa6727ac149a2d6277f0d240fd" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall 0.4.1", + "windows-sys", +] + +[[package]] +name = "flate2" +version = "1.0.28" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "46303f565772937ffe1d394a4fac6f411c6013172fadde9dcdb1e147a086940e" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + +[[package]] +name = "float-ord" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8ce81f49ae8a0482e4c55ea62ebbd7e5a686af544c00b9d090bba3ff9be97b3d" + +[[package]] +name = "font-kit" +version = "0.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2845a73bbd781e691ab7c2a028c579727cd254942e8ced57ff73e0eafd60de87" +dependencies = [ + "bitflags 2.6.0", + "byteorder", + "core-foundation", + "core-graphics", + "core-text", + "dirs-next", + "dwrote", + "float-ord", + "freetype-sys", + "lazy_static", + "libc", + "log", + "pathfinder_geometry", + "pathfinder_simd", + "walkdir", + "winapi", + "yeslogic-fontconfig-sys", +] + +[[package]] +name = "foreign-types" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f6f339eb8adc052cd2ca78910fda869aefa38d22d5cb648e6485e4d3fc06f3b1" +dependencies = [ + "foreign-types-shared 0.1.1", +] + +[[package]] +name = "foreign-types" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d737d9aa519fb7b749cbc3b962edcf310a8dd1f4b67c91c4f83975dbdd17d965" +dependencies = [ + "foreign-types-macros", + "foreign-types-shared 0.3.1", +] + +[[package]] +name = "foreign-types-macros" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1a5c6c585bc94aaf2c7b51dd4c2ba22680844aba4c687be581871a6f518c5742" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.68", +] + +[[package]] +name = "foreign-types-shared" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "00b0228411908ca8685dba7fc2cdd70ec9990a6e753e89b6ac91a84c40fbaf4b" + +[[package]] +name = "foreign-types-shared" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa9a19cbb55df58761df49b23516a86d432839add4af60fc256da840f66ed35b" + +[[package]] +name = "form_urlencoded" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e13624c2627564efccf4934284bdd98cbaa14e79b0b5a141218e507b3a823456" +dependencies = [ + "percent-encoding", +] + +[[package]] +name = "freetype-sys" +version = "0.20.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0e7edc5b9669349acfda99533e9e0bcf26a51862ab43b08ee7745c55d28eb134" +dependencies = [ + "cc", + "libc", + "pkg-config", +] + [[package]] name = "getrandom" version = "0.2.15" @@ -198,17 +587,30 @@ dependencies = [ "wasi", ] +[[package]] +name = "gif" +version = "0.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "80792593675e051cf94a4b111980da2ba60d4a83e43e0048c5693baab3977045" +dependencies = [ + "color_quant", + "weezl", +] + [[package]] name = "gstools-core" version = "1.0.0" dependencies = [ - "approx", + "approx 0.5.1", "criterion", "ndarray", + "ndarray-linalg", "ndarray-rand", "numpy", + "plotters", "pyo3", "rayon", + "special", ] [[package]] @@ -233,6 +635,53 @@ version = "0.3.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d231dfb89cfffdbc30e7fc41579ed6066ad03abda9e567ccafae602b97ec5024" +[[package]] +name = "iana-time-zone" +version = "0.1.60" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e7ffbb5a1b541ea2561f8c41c087286cc091e21e556a4f09a8f6cbf17b69b141" +dependencies = [ + "android_system_properties", + "core-foundation-sys", + "iana-time-zone-haiku", + "js-sys", + "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 = "idna" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "634d9b1461af396cad843f47fdba5597a4f9e6ddd4bfb6ff5d85028c25cb12f6" +dependencies = [ + "unicode-bidi", + "unicode-normalization", +] + +[[package]] +name = "image" +version = "0.24.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "034bbe799d1909622a74d1193aa50147769440040ff36cb2baa947609b0a4e23" +dependencies = [ + "bytemuck", + "byteorder", + "color_quant", + "jpeg-decoder", + "num-traits", + "png", +] + [[package]] name = "indoc" version = "2.0.5" @@ -265,6 +714,12 @@ version = "1.0.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "49f1f14873335454500d59611f1cf4a4b0f786f9ac11f4312a78e4cf2566695b" +[[package]] +name = "jpeg-decoder" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f5d4a7da358eff58addd2877a45865158f0d78c911d43a5784ceb7bbf52833b0" + [[package]] name = "js-sys" version = "0.3.69" @@ -274,18 +729,95 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "katexit" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb1304c448ce2c207c2298a34bc476ce7ae47f63c23fa2b498583b26be9bc88c" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "lambert_w" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8918581f60bc58a272d0c752ccee83645cee54f25192fe66f12c2ef653ee80ad" +dependencies = [ + "fast_polynomial", + "libm", +] + +[[package]] +name = "lapack-sys" +version = "0.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "447f56c85fb410a7a3d36701b2153c1018b1d2b908c5fbaf01c1b04fac33bcbe" +dependencies = [ + "libc", +] + +[[package]] +name = "lax" +version = "0.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f96a229d9557112e574164f8024ce703625ad9f88a90964c1780809358e53da" +dependencies = [ + "cauchy", + "katexit", + "lapack-sys", + "num-traits", + "openblas-src", + "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.155" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "97b3888a4aecf77e811145cadf6eef5901f4782c53886191b2f693f24761847c" +[[package]] +name = "libloading" +version = "0.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0c2a198fb6b0eada2a8df47933734e6d35d350665a33a3593d7164fa52c75c19" +dependencies = [ + "cfg-if", + "windows-targets", +] + [[package]] name = "libm" version = "0.2.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" +[[package]] +name = "libredox" +version = "0.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85c833ca1e66078851dba29046874e38f08b2c883700aa29a03ddd3b23814ee8" +dependencies = [ + "bitflags 2.6.0", + "libc", + "redox_syscall 0.4.1", +] + +[[package]] +name = "linux-raw-sys" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78b3ae25bc7c8c38cec158d1f2757ee79e9b3740fbc7ccf0e59e4b08d793fa89" + [[package]] name = "lock_api" version = "0.4.12" @@ -327,13 +859,43 @@ dependencies = [ "autocfg", ] +[[package]] +name = "miniz_oxide" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9d811f3e15f28568be3407c8e7fdb6514c1cda3cb30683f15b6a1a1dc4ea14a7" +dependencies = [ + "adler", + "simd-adler32", +] + +[[package]] +name = "native-tls" +version = "0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a8614eb2c83d59d1c8cc974dd3f920198647674a0a035e1af1fa58707e317466" +dependencies = [ + "libc", + "log", + "openssl", + "openssl-probe", + "openssl-sys", + "schannel", + "security-framework", + "security-framework-sys", + "tempfile", +] + [[package]] name = "ndarray" version = "0.15.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" dependencies = [ - "approx", + "approx 0.4.0", + "approx 0.5.1", + "cblas-sys", + "libc", "matrixmultiply", "num-complex", "num-integer", @@ -342,6 +904,22 @@ dependencies = [ "rayon", ] +[[package]] +name = "ndarray-linalg" +version = "0.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b0e8dda0c941b64a85c5deb2b3e0144aca87aced64678adfc23eacea6d2cc42" +dependencies = [ + "cauchy", + "katexit", + "lax", + "ndarray", + "num-complex", + "num-traits", + "rand", + "thiserror", +] + [[package]] name = "ndarray-rand" version = "0.14.0" @@ -360,6 +938,8 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" dependencies = [ "num-traits", + "rand", + "serde", ] [[package]] @@ -408,6 +988,76 @@ version = "11.1.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0ab1bc2a289d34bd04a330323ac98a1b4bc82c9d9fcb1e66b63caa84da26b575" +[[package]] +name = "openblas-build" +version = "0.10.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d4b6b44095098cafc71915cfac3427135b6dd2ea85820a7d94a5871cb0d1e169" +dependencies = [ + "anyhow", + "flate2", + "native-tls", + "tar", + "thiserror", + "ureq", + "walkdir", +] + +[[package]] +name = "openblas-src" +version = "0.10.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa4958649f766a1013db4254a852cdf2836764869b6654fa117316905f537363" +dependencies = [ + "dirs", + "openblas-build", + "vcpkg", +] + +[[package]] +name = "openssl" +version = "0.10.66" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9529f4786b70a3e8c61e11179af17ab6188ad8d0ded78c5529441ed39d4bd9c1" +dependencies = [ + "bitflags 2.6.0", + "cfg-if", + "foreign-types 0.3.2", + "libc", + "once_cell", + "openssl-macros", + "openssl-sys", +] + +[[package]] +name = "openssl-macros" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a948666b637a0f465e8564c73e89d4dde00d72d4d473cc972f390fc3dcee7d9c" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.68", +] + +[[package]] +name = "openssl-probe" +version = "0.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ff011a302c396a5197692431fc1948019154afc178baf7d8e37367442a4601cf" + +[[package]] +name = "openssl-sys" +version = "0.9.103" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f9e8deee91df40a943c71b917e5874b951d32a802526c85721ce3b776c929d6" +dependencies = [ + "cc", + "libc", + "pkg-config", + "vcpkg", +] + [[package]] name = "parking_lot" version = "0.12.3" @@ -426,20 +1076,58 @@ checksum = "1e401f977ab385c9e4e3ab30627d6f26d00e2c73eef317493c4ec6d468726cf8" dependencies = [ "cfg-if", "libc", - "redox_syscall", + "redox_syscall 0.5.2", "smallvec", "windows-targets", ] +[[package]] +name = "pathfinder_geometry" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b7b7e7b4ea703700ce73ebf128e1450eb69c3a8329199ffbfb9b2a0418e5ad3" +dependencies = [ + "log", + "pathfinder_simd", +] + +[[package]] +name = "pathfinder_simd" +version = "0.5.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5cf07ef4804cfa9aea3b04a7bbdd5a40031dbb6b4f2cbaf2b011666c80c5b4f2" +dependencies = [ + "rustc_version", +] + +[[package]] +name = "percent-encoding" +version = "2.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3148f5046208a5d56bcfc03053e3ca6334e51da8dfb19b6cdc8b306fae3283e" + +[[package]] +name = "pkg-config" +version = "0.3.30" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d231b230927b5e4ad203db57bbcbee2802f6bce620b1e4a9024a07d94e2907ec" + [[package]] name = "plotters" version = "0.3.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a15b6eccb8484002195a3e44fe65a4ce8e93a625797a063735536fd59cb01cf3" dependencies = [ + "chrono", + "font-kit", + "image", + "lazy_static", "num-traits", + "pathfinder_geometry", "plotters-backend", + "plotters-bitmap", "plotters-svg", + "ttf-parser", "wasm-bindgen", "web-sys", ] @@ -450,6 +1138,17 @@ version = "0.3.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "414cec62c6634ae900ea1c56128dfe87cf63e7caece0852ec76aba307cebadb7" +[[package]] +name = "plotters-bitmap" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f7e7f6fb8302456d7c264a94dada86f76d76e1a03e2294ee86ca7da92983b0a6" +dependencies = [ + "gif", + "image", + "plotters-backend", +] + [[package]] name = "plotters-svg" version = "0.3.6" @@ -459,6 +1158,19 @@ dependencies = [ "plotters-backend", ] +[[package]] +name = "png" +version = "0.17.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78c2378060fb13acff3ba0325b83442c1d2c44fbb76df481160ddc1687cce160" +dependencies = [ + "bitflags 1.3.2", + "crc32fast", + "fdeflate", + "flate2", + "miniz_oxide", +] + [[package]] name = "portable-atomic" version = "1.6.0" @@ -527,7 +1239,7 @@ dependencies = [ "proc-macro2", "pyo3-macros-backend", "quote", - "syn", + "syn 2.0.68", ] [[package]] @@ -540,7 +1252,7 @@ dependencies = [ "proc-macro2", "pyo3-build-config", "quote", - "syn", + "syn 2.0.68", ] [[package]] @@ -618,13 +1330,33 @@ dependencies = [ "crossbeam-utils", ] +[[package]] +name = "redox_syscall" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4722d768eff46b75989dd134e5c353f0d6296e5aaa3132e776cbdb56be7731aa" +dependencies = [ + "bitflags 1.3.2", +] + [[package]] name = "redox_syscall" version = "0.5.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c82cf8cff14456045f55ec4241383baeff27af886adb72ffb2162f99911de0fd" dependencies = [ - "bitflags", + "bitflags 2.6.0", +] + +[[package]] +name = "redox_users" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a18479200779601e498ada4e8c1e1f50e3ee19deb0259c25825a98b5603b2cb4" +dependencies = [ + "getrandom", + "libredox", + "thiserror", ] [[package]] @@ -662,6 +1394,57 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" +[[package]] +name = "rustc_version" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bfa0f585226d2e68097d4f95d113b15b83a82e819ab25717ec0590d9584ef366" +dependencies = [ + "semver", +] + +[[package]] +name = "rustix" +version = "0.38.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "70dc5ec042f7a43c4a73241207cecc9873a06d45debb38b329f8541d85c2730f" +dependencies = [ + "bitflags 2.6.0", + "errno", + "libc", + "linux-raw-sys", + "windows-sys", +] + +[[package]] +name = "rustls-native-certs" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "04182dffc9091a404e0fc069ea5cd60e5b866c3adf881eff99a32d048242dffa" +dependencies = [ + "openssl-probe", + "rustls-pemfile", + "rustls-pki-types", + "schannel", + "security-framework", +] + +[[package]] +name = "rustls-pemfile" +version = "2.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "196fe16b00e106300d3e45ecfcb764fa292a535d7326a29a5875c579c7417425" +dependencies = [ + "base64", + "rustls-pki-types", +] + +[[package]] +name = "rustls-pki-types" +version = "1.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0a2ce646f8655401bb81e7927b812614bd5d91dbc968696be50603510fcaf0" + [[package]] name = "ryu" version = "1.0.18" @@ -677,12 +1460,50 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "schannel" +version = "0.1.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fbc91545643bcf3a0bbb6569265615222618bdf33ce4ffbbd13c4bbd4c093534" +dependencies = [ + "windows-sys", +] + [[package]] name = "scopeguard" version = "1.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" +[[package]] +name = "security-framework" +version = "2.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "897b2245f0b511c87893af39b033e5ca9cce68824c4d7e7630b5a1d339658d02" +dependencies = [ + "bitflags 2.6.0", + "core-foundation", + "core-foundation-sys", + "libc", + "security-framework-sys", +] + +[[package]] +name = "security-framework-sys" +version = "2.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "75da29fe9b9b08fe9d6b22b5b4bcbc75d8db3aa31e639aa56bb62e9d46bfceaf" +dependencies = [ + "core-foundation-sys", + "libc", +] + +[[package]] +name = "semver" +version = "1.0.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61697e0a1c7e512e84a621326239844a24d8207b4669b41bc18b32ea5cbf988b" + [[package]] name = "serde" version = "1.0.203" @@ -700,7 +1521,7 @@ checksum = "500cbc0ebeb6f46627f50f3f5811ccf6bf00643be300b4c3eabc0ef55dc5b5ba" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.68", ] [[package]] @@ -714,12 +1535,39 @@ dependencies = [ "serde", ] +[[package]] +name = "simd-adler32" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d66dc143e6b11c1eddc06d5c423cfc97062865baf299914ab64caa38182078fe" + [[package]] name = "smallvec" version = "1.13.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" +[[package]] +name = "special" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b2e14d17dd8acbc25b769fb726b81aef7f19f48d92bbc0652a97283799c8b9d7" +dependencies = [ + "lambert_w", + "libm", +] + +[[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.68" @@ -731,12 +1579,56 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "tar" +version = "0.4.41" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cb797dad5fb5b76fcf519e702f4a589483b5ef06567f160c392832c1f5e44909" +dependencies = [ + "filetime", + "libc", + "xattr", +] + [[package]] name = "target-lexicon" version = "0.12.14" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e1fc403891a21bcfb7c37834ba66a547a8f402146eba7265b5a6d88059c9ff2f" +[[package]] +name = "tempfile" +version = "3.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b8fcd239983515c23a32fb82099f97d0b11b8c72f654ed659363a95c3dad7a53" +dependencies = [ + "cfg-if", + "fastrand", + "once_cell", + "rustix", + "windows-sys", +] + +[[package]] +name = "thiserror" +version = "1.0.61" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c546c80d6be4bc6a00c0f01730c08df82eaa7a7a61f11d656526506112cc1709" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.61" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "46c3384250002a6d5af4d114f2845d37b57521033f30d5c3f46c4d70e1197533" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.68", +] + [[package]] name = "tinytemplate" version = "1.2.1" @@ -747,18 +1639,86 @@ dependencies = [ "serde_json", ] +[[package]] +name = "tinyvec" +version = "1.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "445e881f4f6d382d5f27c034e25eb92edd7c784ceab92a0937db7f2e9471b938" +dependencies = [ + "tinyvec_macros", +] + +[[package]] +name = "tinyvec_macros" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20" + +[[package]] +name = "ttf-parser" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "17f77d76d837a7830fe1d4f12b7b4ba4192c1888001c7164257e4bc6d21d96b4" + +[[package]] +name = "unicode-bidi" +version = "0.3.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08f95100a766bf4f8f28f90d77e0a5461bbdb219042e7679bebe79004fed8d75" + [[package]] name = "unicode-ident" version = "1.0.12" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" +[[package]] +name = "unicode-normalization" +version = "0.1.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a56d1686db2308d901306f92a263857ef59ea39678a5458e7cb17f01415101f5" +dependencies = [ + "tinyvec", +] + [[package]] name = "unindent" version = "0.2.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" +[[package]] +name = "ureq" +version = "2.10.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b74fc6b57825be3373f7054754755f03ac3a8f5d70015ccad699ba2029956f4a" +dependencies = [ + "base64", + "flate2", + "log", + "native-tls", + "once_cell", + "rustls-native-certs", + "url", +] + +[[package]] +name = "url" +version = "2.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22784dbdf76fdde8af1aeda5622b546b422b6fc585325248a2bf9f5e41e94d6c" +dependencies = [ + "form_urlencoded", + "idna", + "percent-encoding", +] + +[[package]] +name = "vcpkg" +version = "0.2.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "accd4ea62f7bb7a82fe23066fb0957d48ef677f6eeb8215f372f52e48bb32426" + [[package]] name = "walkdir" version = "2.5.0" @@ -796,7 +1756,7 @@ dependencies = [ "once_cell", "proc-macro2", "quote", - "syn", + "syn 2.0.68", "wasm-bindgen-shared", ] @@ -818,7 +1778,7 @@ checksum = "e94f17b526d0a461a191c78ea52bbce64071ed5c04c9ffe424dcb38f74171bb7" dependencies = [ "proc-macro2", "quote", - "syn", + "syn 2.0.68", "wasm-bindgen-backend", "wasm-bindgen-shared", ] @@ -839,6 +1799,28 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "weezl" +version = "0.1.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53a85b86a771b1c87058196170769dd264f66c0782acf1ae6cc51bfd64b39082" + +[[package]] +name = "winapi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + [[package]] name = "winapi-util" version = "0.1.8" @@ -848,6 +1830,21 @@ dependencies = [ "windows-sys", ] +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" + +[[package]] +name = "windows-core" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "33ab640c8d7e35bf8ba19b884ba838ceb4fba93a4e8c65a9059d08afcfc683d9" +dependencies = [ + "windows-targets", +] + [[package]] name = "windows-sys" version = "0.52.0" @@ -920,3 +1917,35 @@ name = "windows_x86_64_msvc" version = "0.52.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bec47e5bfd1bff0eeaf6d8b485cc1074891a197ab4225d504cb7a1ab88b02bf0" + +[[package]] +name = "wio" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5d129932f4644ac2396cb456385cbf9e63b5b30c6e8dc4820bdca4eb082037a5" +dependencies = [ + "winapi", +] + +[[package]] +name = "xattr" +version = "1.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8da84f1a25939b27f6820d92aed108f83ff920fdf11a7b19366c27c4cda81d4f" +dependencies = [ + "libc", + "linux-raw-sys", + "rustix", +] + +[[package]] +name = "yeslogic-fontconfig-sys" +version = "5.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ffb6b23999a8b1a997bf47c7bb4d19ad4029c3327bb3386ebe0a5ff584b33c7a" +dependencies = [ + "cstr", + "dlib", + "once_cell", + "pkg-config", +] diff --git a/Cargo.toml b/Cargo.toml index 56a6944..e3314af 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -32,7 +32,11 @@ codegen-units = 1 pyo3 = { version = "0.21", features = ["abi3-py310", "extension-module"] } numpy = "0.21" ndarray = { version = "0.15", features = ["rayon", "approx-0_5"] } +# "openblas-static" needs gcc, gfortran, and make installed +ndarray-linalg = { version = "0.16", features = ["openblas-static"] } rayon = "1.10" +special = "0.11" +plotters = "0.3" [dev-dependencies] criterion = { version = "0.5", features = ["html_reports"] } diff --git a/benches/gen_benchmark_inputs.py b/benches/gen_benchmark_inputs.py index b6d9f38..892b3f1 100755 --- a/benches/gen_benchmark_inputs.py +++ b/benches/gen_benchmark_inputs.py @@ -8,7 +8,7 @@ def gen_field_summate(path, seed): - pos_no = 10_000 + pos_no = 1_000 mode_no = 1_000 x = np.linspace(0.0, 10.0, pos_no) y = np.linspace(-5.0, 5.0, pos_no) @@ -58,7 +58,7 @@ def prepare_data(pos_no, cond_no): return krige_mat, k_vec, krige_cond - krige_mat, k_vec, krige_cond = prepare_data(pos_no=10_000, cond_no=500) + krige_mat, k_vec, krige_cond = prepare_data(pos_no=1_000, cond_no=300) np.savetxt(path / 'krige_krige_mat.txt', krige_mat) np.savetxt(path / 'krige_k_vec.txt', k_vec) diff --git a/benches/main.rs b/benches/main.rs index fe0cd38..82bad98 100644 --- a/benches/main.rs +++ b/benches/main.rs @@ -11,6 +11,7 @@ use ndarray_rand::{ }; use gstools_core::field::{summator, summator_fourier, summator_incompr}; +use gstools_core::krige::methods::cdist; use gstools_core::krige::{calculator_field_krige, calculator_field_krige_and_variance}; use gstools_core::variogram::{ variogram_directional, variogram_ma_structured, variogram_structured, variogram_unstructured, @@ -210,10 +211,27 @@ pub fn variogram_benchmark(c: &mut Criterion) { }); } +pub fn covmodel_benchmark(c: &mut Criterion) { + let a_ordered = arr2(&[[0.3, 1.9, 1.1, 3.3, 4.7, 5.3, 10.0, 12.8, 15.764, 20.0, 20.0]]); + c.bench_function("cdist ordered", |b| { + b.iter(|| cdist(a_ordered.view(), a_ordered.view())) + }); + let a_rand = arr2(&[[ + 10.3, 1.9, 11.1, 33.3, 4.7, 15.3, 10.0, 22.8, 5.764, 2.0, 2.1, + ]]); + c.bench_function("cdist ordered", |b| { + b.iter(|| cdist(a_rand.view(), a_rand.view())) + }); + c.bench_function("cdist mix", |b| { + b.iter(|| cdist(a_ordered.view(), a_rand.view())) + }); +} + criterion_group!( benches, field_benchmark, krige_benchmark, - variogram_benchmark + variogram_benchmark, + covmodel_benchmark, ); criterion_main!(benches); diff --git a/src/cov_model.rs b/src/cov_model.rs new file mode 100644 index 0000000..c625df4 --- /dev/null +++ b/src/cov_model.rs @@ -0,0 +1,207 @@ +//! Provide Covariance Models + +use plotters::prelude::*; +use special::Gamma; + +use plotters::{backend::BitMapBackend, drawing::IntoDrawingArea}; + +/// Minimum functionality for a covariance model. +pub trait CovModel { + /// Isotropic correlation function of the model + fn correlation(&self, h: f64) -> f64; + /// Getter for the variance + fn var(&self) -> f64; + /// Getter for the length scale + fn len_scale(&self) -> f64; +} + +/// Spectral density of the model +pub trait SpectralDensity { + /// Spectral density of the model + fn spectral_density(&self, k: f64) -> f64; +} + +/// Isotropic variogram of the model +pub trait Variogram { + /// Isotropic variogram of the model + fn variogram(&self, r: f64) -> f64; +} + +/// Isotropic covariance of the model +pub trait Covariance { + /// Isotropic covariance of the model + fn covariance(&self, r: f64) -> f64; +} + +/// Base variables of a covariance model +pub struct BaseCovModel { + /// spatial dimension + pub dim: usize, + /// variance of the model, not including the nugget + pub var: f64, + /// length scale of the model + pub len_scale: f64, +} +/// Gaussian covariance model +pub struct Gaussian { + /// covariance base variables + pub base: BaseCovModel, + /// rescale factor resulting in integral length scale + pub rescale: f64, +} +/// Stable covariance model +pub struct Stable { + /// covariance base variables + pub base: BaseCovModel, + /// shape parameter, valid values ɑ ∊ (0, 2] + pub alpha: f64, + /// rescale factor resulting in integral length scale + pub rescale: f64, +} + +impl Default for BaseCovModel { + fn default() -> Self { + Self { + dim: 3, + var: 1.0, + len_scale: 1.0, + } + } +} +impl Default for Gaussian { + fn default() -> Self { + Self { + base: BaseCovModel { + dim: 3, + var: 1.0, + len_scale: 1.0, + }, + rescale: std::f64::consts::PI.sqrt() / 2.0, + } + } +} +impl Default for Stable { + fn default() -> Self { + Self { + base: BaseCovModel::default(), + alpha: 1.5, + rescale: 1.0, + } + } +} + +impl CovModel for Gaussian { + fn correlation(&self, h: f64) -> f64 { + (-(h / self.len_scale()).powi(2)).exp() + } + fn var(&self) -> f64 { + self.base.var + } + fn len_scale(&self) -> f64 { + self.base.len_scale / self.rescale + } +} +impl CovModel for Stable { + #![allow(unstable_name_collisions)] + fn correlation(&self, h: f64) -> f64 { + (-(h / self.len_scale()).powf(self.alpha) * (1.0 + 1.0 / self.alpha).gamma()).exp() + } + fn var(&self) -> f64 { + self.base.var + } + fn len_scale(&self) -> f64 { + self.base.len_scale / self.rescale + } +} + +impl SpectralDensity for Gaussian { + fn spectral_density(&self, k: f64) -> f64 { + // TODO use std::f64::consts::FRAC_1_SQRT_PI, once it's stable + self.base.len_scale / 2.0 / std::f64::consts::PI.sqrt().powi(self.base.dim as i32) + * (-(k * self.base.len_scale / 2.0).powi(2)).exp() + } +} + +impl Variogram for Gaussian { + fn variogram(&self, r: f64) -> f64 { + // TODO nugget missing, add `+ self.nugget` + self.base.var * (1.0 - self.correlation(r)) + } +} + +impl Variogram for Stable { + fn variogram(&self, r: f64) -> f64 { + // TODO nugget missing, add `+ self.nugget` + self.base.var * (1.0 - self.correlation(r)) + } +} + +impl Covariance for Gaussian { + fn covariance(&self, r: f64) -> f64 { + self.base.var * self.correlation(r) + } +} + +impl Covariance for Stable { + fn covariance(&self, r: f64) -> f64 { + self.base.var * self.correlation(r) + } +} + +/// Plot the correlation function of a given correlation model. +pub fn plot_cor(model: Box) { + let root_drawing_area = BitMapBackend::new("cor.png", (1024, 768)).into_drawing_area(); + + root_drawing_area.fill(&WHITE).unwrap(); + + let x_max = model.len_scale() * 3.0; + let mut chart = ChartBuilder::on(&root_drawing_area) + .build_cartesian_2d(0.0..x_max, 0.0..1.2) + .unwrap(); + + chart + .draw_series(LineSeries::new( + (0..1000) + .map(|x| x as f64 / 1000.0 * x_max) + .map(|x| (x, model.correlation(x))), + &RED, + )) + .unwrap(); +} + +#[cfg(test)] +mod tests { + use super::*; + + use approx::assert_ulps_eq; + + #[test] + fn test_plotting() { + let gau = Box::new(Gaussian::default()); + plot_cor(gau); + } + #[test] + fn test_covmodels() { + let gau = Gaussian { + base: BaseCovModel { + dim: 1, + var: 0.5, + len_scale: 2.0, + }, + ..Default::default() + }; + assert_ulps_eq!(gau.correlation(0.0), 1.0, max_ulps = 6,); + assert_ulps_eq!(gau.correlation(0.1), 0.9980384309875864, max_ulps = 6,); + assert_ulps_eq!(gau.correlation(1.0), 0.8217249580338771, max_ulps = 6,); + assert_ulps_eq!(gau.correlation(10.0), 2.96925699656481e-09, max_ulps = 6,); + assert_ulps_eq!(gau.correlation(-0.5), 0.9520979267837046, max_ulps = 6,); + + assert_ulps_eq!(gau.covariance(0.0), 0.5, max_ulps = 6,); + assert_ulps_eq!(gau.covariance(0.1), 0.4990192154937932, max_ulps = 6,); + assert_ulps_eq!(gau.covariance(1.0), 0.41086247901693856, max_ulps = 6,); + assert_ulps_eq!(gau.covariance(10.0), 1.4846284982824e-09, max_ulps = 6,); + assert_ulps_eq!(gau.covariance(-0.5), 0.4760489633918523, max_ulps = 6,); + + // TODO test Stable model + } +} diff --git a/src/krige/methods.rs b/src/krige/methods.rs new file mode 100644 index 0000000..e2020fe --- /dev/null +++ b/src/krige/methods.rs @@ -0,0 +1,405 @@ +//! Different kriging algorithms implemented here. So far, only +//! simple kriging is implemented. + +use ndarray::{Array, Array1, Array2, ArrayView1, ArrayView2, Zip}; +use ndarray_linalg::{norm::Norm, solve::Inverse}; + +use crate::calculator_field_krige_and_variance; +use crate::cov_model::{CovModel, Covariance}; + +/// Compute the L2 distance between each pair of the two input arrays. +/// +/// # Arguments +/// +/// * `a` - first input array +/// * `b` - second input array +/// +/// # Returns +/// +/// * `res` - the pairwise distance matrix +pub fn cdist(a: ArrayView2<'_, f64>, b: ArrayView2<'_, f64>) -> Array2 { + let mut res = Array::::zeros((a.dim().1, b.dim().1)); + // TODO benchmark if parallelization helps here + Zip::from(res.rows_mut()) + .and(a.columns()) + .for_each(|rr, a| { + Zip::from(rr).and(b.columns()).for_each(|r, b| { + *r = (&b - &a).norm_l2(); + }) + }); + res +} + +/// Generate the kriging field with simple kriging. +pub fn simple( + pos: ArrayView2<'_, f64>, + cov_model: &T, + cond_pos: ArrayView2<'_, f64>, + cond_val: ArrayView1<'_, f64>, + mean: f64, +) -> (Array1, Array1) +where + T: CovModel + Covariance, +{ + // TODO calculate pseudo inverse, instead of inverse + let krig_mat = match cdist(cond_pos, cond_pos) + .mapv(|d| cov_model.covariance(d)) + .inv() + { + Ok(inv) => inv, + Err(error) => { + panic!("{error:?}\nKrige matrix does not have an inverse. Pseudo inverse needs to be implemented.") + } + }; + + // calculate the RHS of kriging system + let krig_vecs = cdist(cond_pos, pos).mapv(|d| cov_model.covariance(d)); + let krige_cond = cond_val.mapv(|v| v - mean); + let (field, error) = calculator_field_krige_and_variance( + krig_mat.view(), + krig_vecs.view(), + krige_cond.view(), + None, + ); + let field = field.mapv(|f| f + mean); + let error = error.mapv(|e| (cov_model.var() - e).max(0.0)); + (field, error) +} + +#[cfg(test)] +mod tests { + use crate::{ + cov_model::{BaseCovModel, Gaussian}, + krige::calculator_field_krige, + }; + + use super::*; + + use approx::assert_ulps_eq; + use ndarray::{arr1, arr2, Array2}; + + struct Setup { + krig_mat: Array2, + krig_vecs: Array2, + cond: Array1, + } + + impl Setup { + fn new() -> Self { + Self { + krig_mat: arr2(&[ + [ + 5.00000000068981e-01, + -5.87287095364834e-06, + 7.82325812566282e-12, + ], + [ + -5.87287095378827e-06, + 5.00000000070158e-01, + -7.67370103394336e-07, + ], + [ + 7.82331319334681e-12, + -7.67370103410243e-07, + 5.00000000001178e-01, + ], + ]), + krig_vecs: arr2(&[ + [ + 3.00650970845165e-01, + 7.92958674144233e-11, + 7.34102993092809e-02, + 1.10371060304999e-08, + 2.00114256042442e-01, + 7.23018134159345e-03, + ], + [ + 5.51416575736629e-09, + 4.79656668238205e-09, + 3.91247964853073e-03, + 3.59846942149471e-11, + 2.10720573114332e-10, + 4.83625846265317e-04, + ], + [ + 7.08796598544206e-13, + 1.09700007286403e-01, + 2.46322359027701e-05, + 1.75889992745405e-07, + 3.05671083940413e-17, + 2.38513785599550e-11, + ], + ]), + cond: arr1(&[ + -1.27755407195723e+00, + 1.15554040655641e+00, + 8.47374235895458e-01, + ]), + } + } + } + + #[test] + fn test_calculator_field_krige_and_variance() { + let setup = Setup::new(); + + let (kr_field, kr_error) = calculator_field_krige_and_variance( + setup.krig_mat.view(), + setup.krig_vecs.view(), + setup.cond.view(), + None, + ); + assert_ulps_eq!( + kr_field, + arr1(&[ + -0.19205097317842723, + 0.04647838537175125, + -0.04462233428403452, + 0.0000000674926344864219, + -0.12782974926973434, + -0.0043390949462510245 + ]), + max_ulps = 6, + ); + assert_ulps_eq!( + kr_error, + arr1(&[ + 0.04519550314128594, + 0.006017045799331816, + 0.0027021867008690937, + 0.000000000000015529554261898964, + 0.020022857738471924, + 0.00002625466702800745 + ]), + max_ulps = 6, + ); + } + + #[test] + fn test_calculator_field_krige() { + let setup = Setup::new(); + + assert_ulps_eq!( + calculator_field_krige( + setup.krig_mat.view(), + setup.krig_vecs.view(), + setup.cond.view(), + None + ), + arr1(&[ + -0.19205097317842723, + 0.04647838537175125, + -0.04462233428403452, + 0.0000000674926344864219, + -0.12782974926973434, + -0.0043390949462510245 + ]), + max_ulps = 6, + ); + } + + #[test] + fn test_linalg_1d() { + let a = arr2(&[[0.3, 1.9, 1.1, 3.3, 4.7]]); + assert_ulps_eq!( + cdist(a.view(), a.view()), + arr2(&[ + [0.0, 1.6, 0.8, 3.0, 4.4], + [1.6, 0.0, 0.8, 1.4, 2.8], + [0.8, 0.8, 0.0, 2.2, 3.6], + [3., 1.4, 2.2, 0.0, 1.4], + [4.4, 2.8, 3.6, 1.4, 0.0], + ]), + max_ulps = 6, + ); + let b = Array::linspace(0.0, 15.0, 4); + let b = b.broadcast((1, b.len())).unwrap(); + assert_ulps_eq!( + cdist(a.view(), b.view()), + arr2(&[ + [0.3, 4.7, 9.7, 14.7], + [1.9, 3.1, 8.1, 13.1], + [1.1, 3.9, 8.9, 13.9], + [3.3, 1.7, 6.7, 11.7], + [4.7, 0.3, 5.3, 10.3], + ]), + max_ulps = 6, + ); + } + #[test] + fn test_linalg_2d() { + let a = arr2(&[[0.3, 1.9, 1.1, 3.3, 4.7], [1.2, 0.6, 3.2, 4.4, 3.8]]); + assert_ulps_eq!( + cdist(a.view(), a.view()), + arr2(&[ + [ + 0.0, + 1.708800749063506, + 2.154065922853802, + 4.386342439892262, + 5.110772935672255 + ], + [ + 1.708800749063506, + 0.0, + 2.7202941017470885, + 4.049691346263318, + 4.252058325093859 + ], + [ + 2.154065922853802, + 2.7202941017470885, + 0.0, + 2.5059928172283334, + 3.6496575181789317 + ], + [ + 4.386342439892262, + 4.049691346263318, + 2.5059928172283334, + 0.0, + 1.5231546211727822 + ], + [ + 5.110772935672255, + 4.252058325093859, + 3.6496575181789317, + 1.5231546211727822, + 0.0 + ] + ]), + max_ulps = 6, + ); + let b = arr2(&[[0., 0., 5., 5.], [-5., 0., 5., -5.]]); + assert_ulps_eq!( + cdist(a.view(), b.view()), + arr2(&[ + [ + 6.2072538211353985, + 1.2369316876852983, + 6.04400529450463, + 7.780102827083971 + ], + [ + 5.913543776789007, + 1.9924858845171274, + 5.3823786563191565, + 6.400781202322104 + ], + [ + 8.27345151674922, + 3.3837848631377265, + 4.295346318982906, + 9.080198235721507 + ], + [ + 9.962429422585638, + 5.5, + 1.8027756377319946, + 9.552486587271401 + ], + [ + 9.976472322419385, + 6.04400529450463, + 1.2369316876852983, + 8.805112151472008 + ] + ]), + max_ulps = 6, + ); + } + #[test] + fn test_krige_1d() { + let pos = Array::linspace(0.0, 15.0, 4); + let pos = pos.broadcast((1, pos.len())).unwrap(); + + let model = Gaussian { + base: BaseCovModel { + dim: 1, + var: 0.5, + len_scale: 2.0, + }, + ..Default::default() + }; + let cond_pos = arr2(&[[0.3, 1.9, 1.1, 3.3, 4.7]]); + let cond_val = arr1(&[0.47, 0.56, 0.74, 1.47, 1.74]); + let (field, error) = simple(pos.view(), &model, cond_pos.view(), cond_val.view(), 1.0); + assert_ulps_eq!( + field, + arr1(&[ + 0.18792146005726795, + 1.485600187143261, + 0.9892592457925125, + 0.9999999972832743 + ]), + max_ulps = 6, + ); + assert_ulps_eq!( + error, + arr1(&[ + 0.0007866097958594831, + 0.0036412749665615807, + 0.499975563693493, + 0.5 + ]), + max_ulps = 6, + ); + } + #[test] + fn test_krige_2d() { + let model = Gaussian { + base: BaseCovModel { + dim: 2, + var: 0.8, + len_scale: 3.0, + }, + ..Default::default() + }; + // TODO maybe use meshgrid? + // https://github.com/rust-ndarray/ndarray/issues/1355 + let pos = arr2(&[ + [0., 0., 0., 5., 5., 5., 10., 10., 10., 15., 15., 15.], + [-5., 0., 5., -5., 0., 5., -5., 0., 5., -5., 0., 5.], + ]); + let cond_pos = arr2(&[[0.3, 1.9, 1.1, 3.3, 4.7], [1.2, 0.6, 3.2, 4.4, 3.8]]); + let cond_val = arr1(&[0.47, 0.56, 0.74, 1.47, 1.74]); + let (field, error) = simple(pos.view(), &model, cond_pos.view(), cond_val.view(), 0.0); + + assert_ulps_eq!( + field, + arr1(&[ + 0.009557358693436949, + 0.3277615223971525, + 0.3145520862286027, + 0.0021019598010470313, + 0.4966425910922256, + 1.5162561066738178, + 0.00015835319007678262, + 0.04079589391922381, + 0.12775437308455703, + 1.8320948846857e-7, + 4.501362280436092e-5, + 0.00014020995740099762 + ]), + max_ulps = 6, + ); + assert_ulps_eq!( + error, + arr1(&[ + 0.797396133570957, + 0.10565469106661851, + 0.36630191151733443, + 0.7987743258653218, + 0.5324374361056077, + 0.1548205035388157, + 0.7999998997165715, + 0.7987431018679746, + 0.7901204798243383, + 0.7999999999999668, + 0.7999999980743688, + 0.7999999824714752 + ]), + max_ulps = 6, + ); + } +} diff --git a/src/krige.rs b/src/krige/mod.rs similarity index 52% rename from src/krige.rs rename to src/krige/mod.rs index bb2f182..af6deea 100644 --- a/src/krige.rs +++ b/src/krige/mod.rs @@ -1,8 +1,10 @@ -//! Perform the kriging matrix operations. +//! Provides different kriging functionality use ndarray::{Array1, ArrayView1, ArrayView2, Zip}; use rayon::prelude::*; +pub mod methods; + /// Calculate the interpolated field and also return the variance. /// /// # Arguments @@ -116,131 +118,3 @@ pub fn calculator_field_krige( }) }) } - -#[cfg(test)] -mod tests { - use super::*; - - use approx::assert_ulps_eq; - use ndarray::{arr1, arr2, Array2}; - - struct Setup { - krig_mat: Array2, - krig_vecs: Array2, - cond: Array1, - } - - impl Setup { - fn new() -> Self { - Self { - krig_mat: arr2(&[ - [ - 5.00000000068981e-01, - -5.87287095364834e-06, - 7.82325812566282e-12, - ], - [ - -5.87287095378827e-06, - 5.00000000070158e-01, - -7.67370103394336e-07, - ], - [ - 7.82331319334681e-12, - -7.67370103410243e-07, - 5.00000000001178e-01, - ], - ]), - krig_vecs: arr2(&[ - [ - 3.00650970845165e-01, - 7.92958674144233e-11, - 7.34102993092809e-02, - 1.10371060304999e-08, - 2.00114256042442e-01, - 7.23018134159345e-03, - ], - [ - 5.51416575736629e-09, - 4.79656668238205e-09, - 3.91247964853073e-03, - 3.59846942149471e-11, - 2.10720573114332e-10, - 4.83625846265317e-04, - ], - [ - 7.08796598544206e-13, - 1.09700007286403e-01, - 2.46322359027701e-05, - 1.75889992745405e-07, - 3.05671083940413e-17, - 2.38513785599550e-11, - ], - ]), - cond: arr1(&[ - -1.27755407195723e+00, - 1.15554040655641e+00, - 8.47374235895458e-01, - ]), - } - } - } - - #[test] - fn test_calculator_field_krige_and_variance() { - let setup = Setup::new(); - - let (kr_field, kr_error) = calculator_field_krige_and_variance( - setup.krig_mat.view(), - setup.krig_vecs.view(), - setup.cond.view(), - None, - ); - assert_ulps_eq!( - kr_field, - arr1(&[ - -0.19205097317842723, - 0.04647838537175125, - -0.04462233428403452, - 0.0000000674926344864219, - -0.12782974926973434, - -0.0043390949462510245 - ]), - max_ulps = 6, - ); - assert_ulps_eq!( - kr_error, - arr1(&[ - 0.04519550314128594, - 0.006017045799331816, - 0.0027021867008690937, - 0.000000000000015529554261898964, - 0.020022857738471924, - 0.00002625466702800745 - ]), - max_ulps = 6, - ); - } - - #[test] - fn test_calculator_field_krige() { - let setup = Setup::new(); - - assert_ulps_eq!( - calculator_field_krige( - setup.krig_mat.view(), - setup.krig_vecs.view(), - setup.cond.view(), - None - ), - arr1(&[ - -0.19205097317842723, - 0.04647838537175125, - -0.04462233428403452, - 0.0000000674926344864219, - -0.12782974926973434, - -0.0043390949462510245 - ]), - max_ulps = 6, - ); - } -} diff --git a/src/lib.rs b/src/lib.rs index c1a62e1..5bbd14f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -21,6 +21,7 @@ use variogram::{ variogram_directional, variogram_ma_structured, variogram_structured, variogram_unstructured, }; +pub mod cov_model; pub mod field; pub mod krige; mod short_vec;