commit aaa212b9a9a7c53d273c1a4e0e86ed7add7839b7
Author: Pascal Engélibert <tuxmain@zettascript.org>
Date: Thu Nov 16 23:28:16 2023 +0100
Initial commit
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..ea8c4bf
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+/target
diff --git a/Cargo.lock b/Cargo.lock
new file mode 100644
index 0000000..c061f88
--- /dev/null
+++ b/Cargo.lock
@@ -0,0 +1,1199 @@
+# This file is automatically @generated by Cargo.
+# 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 = "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 = "approx"
+version = "0.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "autocfg"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa"
+
+[[package]]
+name = "bitflags"
+version = "1.3.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a"
+
+[[package]]
+name = "bitflags"
+version = "2.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "327762f6e5a765692301e5bb513e0d9fef63be86bbc14528052b1cd3e6f03e07"
+
+[[package]]
+name = "bumpalo"
+version = "3.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7f30e7476521f6f8af1a1c4c0b8cc94f0bee37d91763d0ca2665f299b6cd8aec"
+
+[[package]]
+name = "bytemuck"
+version = "1.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "374d28ec25809ee0e23827c2ab573d729e293f281dfe393500e7ad618baa61c6"
+
+[[package]]
+name = "byteorder"
+version = "1.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b"
+
+[[package]]
+name = "cc"
+version = "1.0.83"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f1174fb0b6ec23863f8b971027804a42614e347eafb0a95bf0b12cdae21fc4d0"
+dependencies = [
+ "libc",
+]
+
+[[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.31"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7f2c685bad3eb3d45a01354cedb7d5faa66194d1d58ba6e267a8de788f79db38"
+dependencies = [
+ "android-tzdata",
+ "iana-time-zone",
+ "js-sys",
+ "num-traits",
+ "wasm-bindgen",
+ "windows-targets",
+]
+
+[[package]]
+name = "cmake"
+version = "0.1.50"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a31c789563b815f77f4250caee12365734369f942439b7defd71e18a48197130"
+dependencies = [
+ "cc",
+]
+
+[[package]]
+name = "color_quant"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3d7b894f5411737b7867f4827955924d7c254fc9f4d91a6aad6b097804b1018b"
+
+[[package]]
+name = "const-cstr"
+version = "0.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ed3d0b5ff30645a68f35ece8cea4556ca14ef8a1651455f789a099a0513532a6"
+
+[[package]]
+name = "core-foundation"
+version = "0.9.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "194a7a9e6de53fa55116934067c844d9d749312f75c6f6d0980e8c252f8c2146"
+dependencies = [
+ "core-foundation-sys",
+ "libc",
+]
+
+[[package]]
+name = "core-foundation-sys"
+version = "0.8.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e496a50fda8aacccc86d7529e2c1e0892dbd0f898a6b5645b5561b89c3210efa"
+
+[[package]]
+name = "core-graphics"
+version = "0.22.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2581bbab3b8ffc6fcbd550bf46c355135d16e9ff2a6ea032ad6b9bf1d7efe4fb"
+dependencies = [
+ "bitflags 1.3.2",
+ "core-foundation",
+ "core-graphics-types",
+ "foreign-types",
+ "libc",
+]
+
+[[package]]
+name = "core-graphics-types"
+version = "0.1.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2bb142d41022986c1d8ff29103a1411c8a3dfad3552f87a4f8dc50d61d4f4e33"
+dependencies = [
+ "bitflags 1.3.2",
+ "core-foundation",
+ "libc",
+]
+
+[[package]]
+name = "core-text"
+version = "19.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "99d74ada66e07c1cefa18f8abfba765b486f250de2e4a999e5727fc0dd4b4a25"
+dependencies = [
+ "core-foundation",
+ "core-graphics",
+ "foreign-types",
+ "libc",
+]
+
+[[package]]
+name = "crc32fast"
+version = "1.3.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b540bd8bc810d3885c6ea91e2018302f68baba2129ab3e88f32389ee9370880d"
+dependencies = [
+ "cfg-if",
+]
+
+[[package]]
+name = "darling"
+version = "0.14.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7b750cb3417fd1b327431a470f388520309479ab0bf5e323505daf0290cd3850"
+dependencies = [
+ "darling_core",
+ "darling_macro",
+]
+
+[[package]]
+name = "darling_core"
+version = "0.14.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "109c1ca6e6b7f82cc233a97004ea8ed7ca123a9af07a8230878fcfda9b158bf0"
+dependencies = [
+ "fnv",
+ "ident_case",
+ "proc-macro2",
+ "quote",
+ "strsim",
+ "syn 1.0.109",
+]
+
+[[package]]
+name = "darling_macro"
+version = "0.14.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a4aab4dbc9f7611d8b55048a3a16d2d010c2c8334e46304b40ac1cc14bf3b48e"
+dependencies = [
+ "darling_core",
+ "quote",
+ "syn 1.0.109",
+]
+
+[[package]]
+name = "derive_builder"
+version = "0.12.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8d67778784b508018359cbc8696edb3db78160bab2c2a28ba7f56ef6932997f8"
+dependencies = [
+ "derive_builder_macro",
+]
+
+[[package]]
+name = "derive_builder_core"
+version = "0.12.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c11bdc11a0c47bc7d37d582b5285da6849c96681023680b906673c5707af7b0f"
+dependencies = [
+ "darling",
+ "proc-macro2",
+ "quote",
+ "syn 1.0.109",
+]
+
+[[package]]
+name = "derive_builder_macro"
+version = "0.12.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ebcda35c7a396850a55ffeac740804b40ffec779b98fffbb1738f4033f0ee79e"
+dependencies = [
+ "derive_builder_core",
+ "syn 1.0.109",
+]
+
+[[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-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 = "epidemics"
+version = "0.1.0"
+dependencies = [
+ "nalgebra",
+ "rand",
+ "rustmodel",
+]
+
+[[package]]
+name = "fast-srgb8"
+version = "1.0.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dd2e7510819d6fbf51a5545c8f922716ecfb14df168a3242f7d33e0239efe6a1"
+
+[[package]]
+name = "fdeflate"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "64d6dafc854908ff5da46ff3f8f473c6984119a2876a383a860246dd7841a868"
+dependencies = [
+ "simd-adler32",
+]
+
+[[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.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7bad48618fdb549078c333a7a8528acb57af271d0433bdecd523eb620628364e"
+
+[[package]]
+name = "fnv"
+version = "1.0.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1"
+
+[[package]]
+name = "font-kit"
+version = "0.11.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "21fe28504d371085fae9ac7a3450f0b289ab71e07c8e57baa3fb68b9e57d6ce5"
+dependencies = [
+ "bitflags 1.3.2",
+ "byteorder",
+ "core-foundation",
+ "core-graphics",
+ "core-text",
+ "dirs-next",
+ "dwrote",
+ "float-ord",
+ "freetype",
+ "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",
+]
+
+[[package]]
+name = "foreign-types-shared"
+version = "0.1.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "00b0228411908ca8685dba7fc2cdd70ec9990a6e753e89b6ac91a84c40fbaf4b"
+
+[[package]]
+name = "freetype"
+version = "0.7.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bee38378a9e3db1cc693b4f88d166ae375338a0ff75cb8263e1c601d51f35dc6"
+dependencies = [
+ "freetype-sys",
+ "libc",
+]
+
+[[package]]
+name = "freetype-sys"
+version = "0.13.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a37d4011c0cc628dfa766fcc195454f4b068d7afdc2adfd28861191d866e731a"
+dependencies = [
+ "cmake",
+ "libc",
+ "pkg-config",
+]
+
+[[package]]
+name = "getrandom"
+version = "0.2.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fe9006bed769170c11f845cf00c7c1e9092aeb3f268e007c3e760ac68008070f"
+dependencies = [
+ "cfg-if",
+ "libc",
+ "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 = "iana-time-zone"
+version = "0.1.58"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8326b86b6cff230b97d0d312a6c40a60726df3332e721f72a1b035f451663b20"
+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 = "ident_case"
+version = "1.0.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39"
+
+[[package]]
+name = "image"
+version = "0.24.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6f3dfdbdd72063086ff443e297b61695500514b1e41095b6fb9a5ab48a70a711"
+dependencies = [
+ "bytemuck",
+ "byteorder",
+ "color_quant",
+ "jpeg-decoder",
+ "num-rational",
+ "num-traits",
+ "png",
+]
+
+[[package]]
+name = "jpeg-decoder"
+version = "0.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bc0000e42512c92e31c2252315bda326620a4e034105e900c98ec492fa077b3e"
+
+[[package]]
+name = "js-sys"
+version = "0.3.65"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "54c0c35952f67de54bb584e9fd912b3023117cbafc0a77d8f3dee1fb5f572fe8"
+dependencies = [
+ "wasm-bindgen",
+]
+
+[[package]]
+name = "lazy_static"
+version = "1.4.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
+
+[[package]]
+name = "libc"
+version = "0.2.150"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "89d92a4743f9a61002fae18374ed11e7973f530cb3a3255fb354818118b2203c"
+
+[[package]]
+name = "libloading"
+version = "0.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c571b676ddfc9a8c12f1f3d3085a7b163966a8fd8098a90640953ce5f6170161"
+dependencies = [
+ "cfg-if",
+ "windows-sys",
+]
+
+[[package]]
+name = "libredox"
+version = "0.0.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "85c833ca1e66078851dba29046874e38f08b2c883700aa29a03ddd3b23814ee8"
+dependencies = [
+ "bitflags 2.4.1",
+ "libc",
+ "redox_syscall",
+]
+
+[[package]]
+name = "log"
+version = "0.4.20"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b5e6163cb8c49088c2c36f57875e58ccd8c87c7427f7fbd50ea6710b2f3f2e8f"
+
+[[package]]
+name = "matrixmultiply"
+version = "0.3.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
+[[package]]
+name = "miniz_oxide"
+version = "0.7.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e7810e0be55b428ada41041c41f32c9f1a42817901b4ccf45fa3d4b6561e74c7"
+dependencies = [
+ "adler",
+ "simd-adler32",
+]
+
+[[package]]
+name = "nalgebra"
+version = "0.32.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "307ed9b18cc2423f29e83f84fd23a8e73628727990181f18641a8b5dc2ab1caa"
+dependencies = [
+ "approx",
+ "matrixmultiply",
+ "nalgebra-macros",
+ "num-complex",
+ "num-rational",
+ "num-traits",
+ "simba",
+ "typenum",
+]
+
+[[package]]
+name = "nalgebra-macros"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "91761aed67d03ad966ef783ae962ef9bbaca728d2dd7ceb7939ec110fffad998"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 1.0.109",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.4.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1ba157ca0885411de85d6ca030ba7e2a83a28636056c7c699b07c8b6f7383214"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.45"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9"
+dependencies = [
+ "autocfg",
+ "num-traits",
+]
+
+[[package]]
+name = "num-rational"
+version = "0.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.17"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "39e3200413f237f41ab11ad6d161bc7239c84dcb631773ccd7de3dfe4b5c267c"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "once_cell"
+version = "1.18.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d"
+
+[[package]]
+name = "palette"
+version = "0.7.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b2e2f34147767aa758aa649415b50a69eeb46a67f9dc7db8011eeb3d84b351dc"
+dependencies = [
+ "approx",
+ "fast-srgb8",
+ "palette_derive",
+ "phf",
+]
+
+[[package]]
+name = "palette_derive"
+version = "0.7.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b7db010ec5ff3d4385e4f133916faacd9dad0f6a09394c92d825b3aed310fa0a"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.39",
+]
+
+[[package]]
+name = "paste"
+version = "1.0.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c"
+
+[[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.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0444332826c70dc47be74a7c6a5fc44e23a7905ad6858d4162b658320455ef93"
+dependencies = [
+ "rustc_version",
+]
+
+[[package]]
+name = "phf"
+version = "0.11.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ade2d8b8f33c7333b51bcf0428d37e217e9f32192ae4772156f65063b8ce03dc"
+dependencies = [
+ "phf_macros",
+ "phf_shared",
+]
+
+[[package]]
+name = "phf_generator"
+version = "0.11.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "48e4cc64c2ad9ebe670cb8fd69dd50ae301650392e81c05f9bfcb2d5bdbc24b0"
+dependencies = [
+ "phf_shared",
+ "rand",
+]
+
+[[package]]
+name = "phf_macros"
+version = "0.11.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3444646e286606587e49f3bcf1679b8cef1dc2c5ecc29ddacaffc305180d464b"
+dependencies = [
+ "phf_generator",
+ "phf_shared",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.39",
+]
+
+[[package]]
+name = "phf_shared"
+version = "0.11.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "90fcb95eef784c2ac79119d1dd819e162b5da872ce6f3c3abe1e8ca1c082f72b"
+dependencies = [
+ "siphasher",
+]
+
+[[package]]
+name = "pkg-config"
+version = "0.3.27"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "26072860ba924cbfa98ea39c8c19b4dd6a4a25423dbdf219c1eca91aa0cf6964"
+
+[[package]]
+name = "plotters"
+version = "0.3.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d2c224ba00d7cadd4d5c660deaf2098e5e80e07846537c51f9cfa4be50c1fd45"
+dependencies = [
+ "chrono",
+ "font-kit",
+ "image",
+ "lazy_static",
+ "num-traits",
+ "pathfinder_geometry",
+ "plotters-backend",
+ "plotters-bitmap",
+ "plotters-svg",
+ "ttf-parser",
+ "wasm-bindgen",
+ "web-sys",
+]
+
+[[package]]
+name = "plotters-backend"
+version = "0.3.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9e76628b4d3a7581389a35d5b6e2139607ad7c75b17aed325f210aa91f4a9609"
+
+[[package]]
+name = "plotters-bitmap"
+version = "0.3.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0cebbe1f70205299abc69e8b295035bb52a6a70ee35474ad10011f0a4efb8543"
+dependencies = [
+ "gif",
+ "image",
+ "plotters-backend",
+]
+
+[[package]]
+name = "plotters-svg"
+version = "0.3.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "38f6d39893cca0701371e3c27294f09797214b86f1fb951b89ade8ec04e2abab"
+dependencies = [
+ "plotters-backend",
+]
+
+[[package]]
+name = "png"
+version = "0.17.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dd75bf2d8dd3702b9707cdbc56a5b9ef42cec752eb8b3bafc01234558442aa64"
+dependencies = [
+ "bitflags 1.3.2",
+ "crc32fast",
+ "fdeflate",
+ "flate2",
+ "miniz_oxide",
+]
+
+[[package]]
+name = "ppv-lite86"
+version = "0.2.17"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de"
+
+[[package]]
+name = "proc-macro2"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "134c189feb4956b20f6f547d2cf727d4c0fe06722b20a0eec87ed445a97f92da"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "quote"
+version = "1.0.33"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae"
+dependencies = [
+ "proc-macro2",
+]
+
+[[package]]
+name = "rand"
+version = "0.8.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404"
+dependencies = [
+ "libc",
+ "rand_chacha",
+ "rand_core",
+]
+
+[[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",
+]
+
+[[package]]
+name = "rand_core"
+version = "0.6.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
+dependencies = [
+ "getrandom",
+]
+
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
+[[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_users"
+version = "0.4.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a18479200779601e498ada4e8c1e1f50e3ee19deb0259c25825a98b5603b2cb4"
+dependencies = [
+ "getrandom",
+ "libredox",
+ "thiserror",
+]
+
+[[package]]
+name = "rustc_version"
+version = "0.4.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bfa0f585226d2e68097d4f95d113b15b83a82e819ab25717ec0590d9584ef366"
+dependencies = [
+ "semver",
+]
+
+[[package]]
+name = "rustmodel"
+version = "0.1.0"
+dependencies = [
+ "derive_builder",
+ "nalgebra",
+ "num-traits",
+ "palette",
+ "plotters",
+]
+
+[[package]]
+name = "safe_arch"
+version = "0.7.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354"
+dependencies = [
+ "bytemuck",
+]
+
+[[package]]
+name = "same-file"
+version = "1.0.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502"
+dependencies = [
+ "winapi-util",
+]
+
+[[package]]
+name = "semver"
+version = "1.0.20"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "836fa6a3e1e547f9a2c4040802ec865b5d85f4014efe00555d7090a3dcaa1090"
+
+[[package]]
+name = "simba"
+version = "0.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "061507c94fc6ab4ba1c9a0305018408e312e17c041eb63bef8aa726fa33aceae"
+dependencies = [
+ "approx",
+ "num-complex",
+ "num-traits",
+ "paste",
+ "wide",
+]
+
+[[package]]
+name = "simd-adler32"
+version = "0.3.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d66dc143e6b11c1eddc06d5c423cfc97062865baf299914ab64caa38182078fe"
+
+[[package]]
+name = "siphasher"
+version = "0.3.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "38b58827f4464d87d377d175e90bf58eb00fd8716ff0a62f80356b5e61555d0d"
+
+[[package]]
+name = "strsim"
+version = "0.10.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73473c0e59e6d5812c5dfe2a064a6444949f089e20eec9a2e5506596494e4623"
+
+[[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.39"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "23e78b90f2fcf45d3e842032ce32e3f2d1545ba6636271dcbf24fa306d87be7a"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "thiserror"
+version = "1.0.50"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f9a7210f5c9a7156bb50aa36aed4c95afb51df0df00713949448cf9e97d382d2"
+dependencies = [
+ "thiserror-impl",
+]
+
+[[package]]
+name = "thiserror-impl"
+version = "1.0.50"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "266b2e40bc00e5a6c09c3584011e08b06f123c00362c92b975ba9843aaaa14b8"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.39",
+]
+
+[[package]]
+name = "ttf-parser"
+version = "0.17.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "375812fa44dab6df41c195cd2f7fecb488f6c09fbaafb62807488cefab642bff"
+
+[[package]]
+name = "typenum"
+version = "1.17.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825"
+
+[[package]]
+name = "unicode-ident"
+version = "1.0.12"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b"
+
+[[package]]
+name = "walkdir"
+version = "2.4.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d71d857dc86794ca4c280d616f7da00d2dbfd8cd788846559a6813e6aa4b54ee"
+dependencies = [
+ "same-file",
+ "winapi-util",
+]
+
+[[package]]
+name = "wasi"
+version = "0.11.0+wasi-snapshot-preview1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423"
+
+[[package]]
+name = "wasm-bindgen"
+version = "0.2.88"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7daec296f25a1bae309c0cd5c29c4b260e510e6d813c286b19eaadf409d40fce"
+dependencies = [
+ "cfg-if",
+ "wasm-bindgen-macro",
+]
+
+[[package]]
+name = "wasm-bindgen-backend"
+version = "0.2.88"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e397f4664c0e4e428e8313a469aaa58310d302159845980fd23b0f22a847f217"
+dependencies = [
+ "bumpalo",
+ "log",
+ "once_cell",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.39",
+ "wasm-bindgen-shared",
+]
+
+[[package]]
+name = "wasm-bindgen-macro"
+version = "0.2.88"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5961017b3b08ad5f3fe39f1e79877f8ee7c23c5e5fd5eb80de95abc41f1f16b2"
+dependencies = [
+ "quote",
+ "wasm-bindgen-macro-support",
+]
+
+[[package]]
+name = "wasm-bindgen-macro-support"
+version = "0.2.88"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c5353b8dab669f5e10f5bd76df26a9360c748f054f862ff5f3f8aae0c7fb3907"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.39",
+ "wasm-bindgen-backend",
+ "wasm-bindgen-shared",
+]
+
+[[package]]
+name = "wasm-bindgen-shared"
+version = "0.2.88"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0d046c5d029ba91a1ed14da14dca44b68bf2f124cfbaf741c54151fdb3e0750b"
+
+[[package]]
+name = "web-sys"
+version = "0.3.65"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5db499c5f66323272151db0e666cd34f78617522fb0c1604d31a27c50c206a85"
+dependencies = [
+ "js-sys",
+ "wasm-bindgen",
+]
+
+[[package]]
+name = "weezl"
+version = "0.1.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9193164d4de03a926d909d3bc7c30543cecb35400c02114792c2cae20d5e2dbb"
+
+[[package]]
+name = "wide"
+version = "0.7.13"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c68938b57b33da363195412cfc5fc37c9ed49aa9cfe2156fde64b8d2c9498242"
+dependencies = [
+ "bytemuck",
+ "safe_arch",
+]
+
+[[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.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f29e6f9198ba0d26b4c9f07dbe6f9ed633e1f3d5b8b414090084349e46a52596"
+dependencies = [
+ "winapi",
+]
+
+[[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.51.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f1f8cf84f35d2db49a46868f947758c7a1138116f7fac3bc844f43ade1292e64"
+dependencies = [
+ "windows-targets",
+]
+
+[[package]]
+name = "windows-sys"
+version = "0.48.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "677d2418bec65e3338edb076e806bc1ec15693c5d0104683f2efe857f61056a9"
+dependencies = [
+ "windows-targets",
+]
+
+[[package]]
+name = "windows-targets"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9a2fa6e2155d7247be68c096456083145c183cbbbc2764150dda45a87197940c"
+dependencies = [
+ "windows_aarch64_gnullvm",
+ "windows_aarch64_msvc",
+ "windows_i686_gnu",
+ "windows_i686_msvc",
+ "windows_x86_64_gnu",
+ "windows_x86_64_gnullvm",
+ "windows_x86_64_msvc",
+]
+
+[[package]]
+name = "windows_aarch64_gnullvm"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2b38e32f0abccf9987a4e3079dfb67dcd799fb61361e53e2882c3cbaf0d905d8"
+
+[[package]]
+name = "windows_aarch64_msvc"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dc35310971f3b2dbbf3f0690a219f40e2d9afcf64f9ab7cc1be722937c26b4bc"
+
+[[package]]
+name = "windows_i686_gnu"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a75915e7def60c94dcef72200b9a8e58e5091744960da64ec734a6c6e9b3743e"
+
+[[package]]
+name = "windows_i686_msvc"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8f55c233f70c4b27f66c523580f78f1004e8b5a8b659e05a4eb49d4166cca406"
+
+[[package]]
+name = "windows_x86_64_gnu"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "53d40abd2583d23e4718fddf1ebec84dbff8381c07cae67ff7768bbf19c6718e"
+
+[[package]]
+name = "windows_x86_64_gnullvm"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0b7b52767868a23d5bab768e390dc5f5c55825b6d30b86c844ff2dc7414044cc"
+
+[[package]]
+name = "windows_x86_64_msvc"
+version = "0.48.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ed94fce61571a4006852b7389a063ab983c02eb1bb37b47f8272ce92d06d9538"
+
+[[package]]
+name = "wio"
+version = "0.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5d129932f4644ac2396cb456385cbf9e63b5b30c6e8dc4820bdca4eb082037a5"
+dependencies = [
+ "winapi",
+]
+
+[[package]]
+name = "yeslogic-fontconfig-sys"
+version = "3.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f2bbd69036d397ebbff671b1b8e4d918610c181c5a16073b96f984a38d08c386"
+dependencies = [
+ "const-cstr",
+ "dlib",
+ "once_cell",
+ "pkg-config",
+]
diff --git a/Cargo.toml b/Cargo.toml
new file mode 100644
index 0000000..5706882
--- /dev/null
+++ b/Cargo.toml
@@ -0,0 +1,29 @@
+[package]
+name = "rustmodel"
+version = "0.1.0"
+edition = "2021"
+description = "simple ODE solving framework"
+authors = ["tuxmain <tuxmain@zettascript.org>"]
+repository = "https://git.txmn.tk/tuxmain/rustmodel"
+license = "AGPL-3.0-only"
+
+[dependencies]
+derive_builder = { version = "0.12.0", optional = true }
+nalgebra = "0.32.3"
+num-traits = "0.2.17"
+palette = { version = "0.7.3", optional = true }
+plotters = { version = "0.3.5", optional = true }
+#plotters = { git = "https://github.com/plotters-rs/plotters" }
+#rand = {version = "0.8.5", optional = true }
+#rayon = "1.8.0"
+#sdl2 = "0.35.2"
+
+[features]
+default = ["plot"]
+
+# opti is about numerical optimization methods, not computing speed!
+#opti = ["nalgebra/rand", "rand"]
+plot = ["derive_builder", "palette", "plotters"]
+
+[workspace]
+members = ["examples/epidemics"]
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..be3f7b2
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,661 @@
+ GNU AFFERO GENERAL PUBLIC LICENSE
+ Version 3, 19 November 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU Affero General Public License is a free, copyleft license for
+software and other kinds of works, specifically designed to ensure
+cooperation with the community in the case of network server software.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+our General Public Licenses are intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ Developers that use our General Public Licenses protect your rights
+with two steps: (1) assert copyright on the software, and (2) offer
+you this License which gives you legal permission to copy, distribute
+and/or modify the software.
+
+ A secondary benefit of defending all users' freedom is that
+improvements made in alternate versions of the program, if they
+receive widespread use, become available for other developers to
+incorporate. Many developers of free software are heartened and
+encouraged by the resulting cooperation. However, in the case of
+software used on network servers, this result may fail to come about.
+The GNU General Public License permits making a modified version and
+letting the public access it on a server without ever releasing its
+source code to the public.
+
+ The GNU Affero General Public License is designed specifically to
+ensure that, in such cases, the modified source code becomes available
+to the community. It requires the operator of a network server to
+provide the source code of the modified version running there to the
+users of that server. Therefore, public use of a modified version, on
+a publicly accessible server, gives the public access to the source
+code of the modified version.
+
+ An older license, called the Affero General Public License and
+published by Affero, was designed to accomplish similar goals. This is
+a different license, not a version of the Affero GPL, but Affero has
+released a new version of the Affero GPL which permits relicensing under
+this license.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU Affero General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Remote Network Interaction; Use with the GNU General Public License.
+
+ Notwithstanding any other provision of this License, if you modify the
+Program, your modified version must prominently offer all users
+interacting with it remotely through a computer network (if your version
+supports such interaction) an opportunity to receive the Corresponding
+Source of your version by providing access to the Corresponding Source
+from a network server at no charge, through some standard or customary
+means of facilitating copying of software. This Corresponding Source
+shall include the Corresponding Source for any work covered by version 3
+of the GNU General Public License that is incorporated pursuant to the
+following paragraph.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the work with which it is combined will remain governed by version
+3 of the GNU General Public License.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU Affero General Public License from time to time. Such new versions
+will be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU Affero General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU Affero General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU Affero General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Affero General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Affero General Public License for more details.
+
+ You should have received a copy of the GNU Affero General Public License
+ along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If your software can interact with users remotely through a computer
+network, you should also make sure that it provides a way for users to
+get its source. For example, if your program is a web application, its
+interface could display a "Source" link that leads users to an archive
+of the code. There are many ways you could offer source, and different
+solutions will be better for different programs; see section 13 for the
+specific requirements.
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU AGPL, see
+<https://www.gnu.org/licenses/>.
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..adcc076
--- /dev/null
+++ b/README.md
@@ -0,0 +1,17 @@
+Simple framework for solving first-order ODEs
+
+## Example
+
+Simulate an epidemics using the compartmental model SIR (Susceptible-Infected-Removed). This command prints tab-separated columns Susceptible and Infected, then writes a plot at `target/sir.png`.
+
+```bash
+cargo run -p epidemics
+```
+
+## License
+
+GNU AGPL v3, CopyLeft 2021-2023 Pascal Engélibert [(why copyleft?)](https://txmn.tk/blog/why-copyleft/)
+
+This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, version 3 of the License.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details.
+You should have received a copy of the GNU Affero General Public License along with this program. If not, see https://www.gnu.org/licenses/.
diff --git a/examples/epidemics/Cargo.toml b/examples/epidemics/Cargo.toml
new file mode 100644
index 0000000..7fdf992
--- /dev/null
+++ b/examples/epidemics/Cargo.toml
@@ -0,0 +1,9 @@
+[package]
+name = "epidemics"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+nalgebra = "0.32.3"
+rand = "0.8.5"
+rustmodel = { path = "../..", features = ["plot"] }
diff --git a/examples/epidemics/src/main.rs b/examples/epidemics/src/main.rs
new file mode 100644
index 0000000..795599b
--- /dev/null
+++ b/examples/epidemics/src/main.rs
@@ -0,0 +1,107 @@
+use nalgebra::vector;
+use rustmodel::prelude::*;
+
+fn main() {
+ // Start with 1% infected
+ let x0 = vector![0.99, 0.01];
+ let dt = 0.1;
+ let nsamples: usize = 400;
+
+ let settings = sir::SirSettings {
+ beta: 0.6,
+ gamma: 0.1,
+ pop: 1.0,
+ };
+
+ let model = sir::Sir {
+ s: settings.clone(),
+ };
+ let solver = ImplicitEuler {
+ dt,
+ tol: 0.000001,
+ niters: 100,
+ };
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+ let mut x = x0;
+ for _ in 1..nsamples {
+ x = solver.f(&model, x);
+ xlist.push(x);
+ }
+
+ xlist.iter().for_each(|x| println!("{}\t{}", x[0], x[1]));
+
+ PlotBuilder::default()
+ .dt(dt)
+ .title(Some(String::from("Epidemics Compartmental Model")))
+ .y_min(Some(0.0))
+ .y_max(Some(1.0))
+ .x_label(Some(String::from("Time")))
+ .build()
+ .unwrap()
+ .plot(
+ "target/sir.png",
+ &xlist,
+ &[("Susceptible", colors::BLUE), ("Infected", colors::RED)],
+ );
+}
+
+/// SIR model without vital dynamics
+pub mod sir {
+ use nalgebra::{base::*, matrix, vector};
+ use rustmodel::prelude::*;
+
+ #[derive(Clone, Debug)]
+ pub struct SirSettings<T> {
+ /// Transmission probability
+ pub beta: T,
+ /// Removal probability
+ pub gamma: T,
+ /// Total population
+ pub pop: T,
+ }
+
+ impl<T> Settings for SirSettings<T> {}
+
+ #[derive(Clone)]
+ pub struct Sir<T> {
+ pub s: SirSettings<T>,
+ }
+
+ impl Model<f64, SirSettings<f64>, 2> for Sir<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ vector![
+ -self.s.beta * x[0] * x[1] / self.s.pop,
+ self.s.beta * x[0] * x[1] / self.s.pop - self.s.gamma * x[1]
+ ]
+ }
+ fn df(&self, x: Vector2<f64>) -> Matrix2<f64> {
+ matrix![
+ -self.s.beta*x[1]/self.s.pop, -self.s.beta*x[0]/self.s.pop;
+ self.s.beta*x[1]/self.s.pop, self.s.beta*x[0]/self.s.pop-self.s.gamma
+ ]
+ }
+ fn get_settings(&self) -> &SirSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut SirSettings<f64> {
+ &mut self.s
+ }
+ }
+
+ impl<T> From<SirSettings<T>> for Vect<T, 3> {
+ fn from(s: SirSettings<T>) -> Self {
+ vector![s.beta, s.gamma, s.pop]
+ }
+ }
+
+ impl<T: Scalar + Copy> From<Vect<T, 3>> for SirSettings<T> {
+ fn from(v: Vect<T, 3>) -> Self {
+ Self {
+ beta: v[0],
+ gamma: v[1],
+ pop: v[2],
+ }
+ }
+ }
+}
diff --git a/examples/old/charts.rs b/examples/old/charts.rs
new file mode 100644
index 0000000..ce0ce2f
--- /dev/null
+++ b/examples/old/charts.rs
@@ -0,0 +1,608 @@
+use crate::{
+ model::{self, Model},
+ opti::GradientDescentOptimizer,
+ solver::Solver,
+ utils::*,
+};
+
+use plotters::prelude::*;
+use rayon::prelude::*;
+
+const CHART_SIZE: (u32, u32) = (1000, 800); //(500, 400);
+const CHART_SIZE_OBJ: (u32, u32) = (960, 960); //(480, 480);
+
+pub fn draw_chart(filename: &str, title: Option<&str>, pop: f64, xlist: &[Vect<f64, 2>], dt: f64) {
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, CHART_SIZE).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .y_label_area_size(30)
+ .x_label_area_size(30)
+ .build_cartesian_2d(0.0f64..xlist.len() as f64 * dt, 0.0f64..1.)
+ .unwrap();
+
+ chart.configure_mesh().x_desc("Time").draw().unwrap();
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist.iter().enumerate().map(|(i, x)| (i as f64 * dt, x[0])),
+ BLUE,
+ ))
+ .unwrap()
+ .label("Susceptible")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], BLUE));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist.iter().enumerate().map(|(i, x)| (i as f64 * dt, x[1])),
+ RED,
+ ))
+ .unwrap()
+ .label("Infected")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RED));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt, pop - x[0] - x[1])),
+ GREEN,
+ ))
+ .unwrap()
+ .label("Removed")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], GREEN));
+
+ chart
+ .configure_series_labels()
+ .background_style(WHITE.mix(0.8))
+ .border_style(BLACK)
+ .draw()
+ .unwrap();
+}
+
+pub fn draw_error_chart(filename: &str, title: Option<&str>, xlist: &[f64]) {
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, CHART_SIZE).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .y_label_area_size(50)
+ .x_label_area_size(30)
+ .build_cartesian_2d(0..xlist.len(), (0.0f64..max(xlist)).log_scale())
+ .unwrap();
+
+ let printer = plotters::data::float::FloatPrettyPrinter {
+ allow_scientific: true,
+ min_decimal: 0,
+ max_decimal: 2,
+ };
+ chart
+ .configure_mesh()
+ .x_desc("Iterations")
+ .y_desc("Mean error")
+ .y_label_formatter(&|y| printer.print(*y))
+ .draw()
+ .unwrap();
+
+ chart
+ .draw_series(LineSeries::new(xlist.iter().copied().enumerate(), BLACK))
+ .unwrap();
+}
+
+pub fn draw_error_chart2(
+ filename: &str,
+ title: Option<&str>,
+ xlist_batch: &[f64],
+ xlist_sto: &[f64],
+) {
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, CHART_SIZE).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .y_label_area_size(50)
+ .x_label_area_size(30)
+ .build_cartesian_2d(
+ 0..xlist_batch.len().max(xlist_sto.len()),
+ (0.0f64..max(xlist_batch).max(max(xlist_sto))).log_scale(),
+ )
+ .unwrap();
+
+ let printer = plotters::data::float::FloatPrettyPrinter {
+ allow_scientific: true,
+ min_decimal: 0,
+ max_decimal: 2,
+ };
+ chart
+ .configure_mesh()
+ .x_desc("Iterations")
+ .y_desc("Mean error")
+ .y_label_formatter(&|y| printer.print(*y))
+ .draw()
+ .unwrap();
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_batch.iter().copied().enumerate(),
+ BLACK,
+ ))
+ .unwrap()
+ .label("Batch")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], BLACK));
+
+ chart
+ .draw_series(LineSeries::new(xlist_sto.iter().copied().enumerate(), RED))
+ .unwrap()
+ .label("Stochastic")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RED));
+
+ chart
+ .configure_series_labels()
+ .background_style(WHITE.mix(0.8))
+ .border_style(BLACK)
+ .position(SeriesLabelPosition::MiddleLeft)
+ .draw()
+ .unwrap();
+}
+
+pub fn plot_objective<
+ R: Solver<f64, model::sir::SirSettings<f64>, model::sir::Sir<f64>, 2> + Clone + Sync,
+>(
+ filename: &str,
+ title: Option<&str>,
+ optimizer: GradientDescentOptimizer<
+ f64,
+ model::sir::SirSettings<f64>,
+ model::sir::Sir<f64>,
+ R,
+ 2,
+ 3,
+ >,
+ ylist_true: &[Vect<f64, 2>],
+ path_batch: Option<&[(f64, f64)]>,
+ path_sto: Option<&[(f64, f64)]>,
+) {
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, CHART_SIZE_OBJ).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .x_label_area_size(30)
+ .y_label_area_size(40)
+ .build_cartesian_2d(0.0f64..1., 0.0f64..1.)
+ .unwrap();
+
+ chart
+ .configure_mesh()
+ .x_desc("beta")
+ .y_desc("gamma")
+ .draw()
+ .unwrap();
+
+ let area = chart.plotting_area();
+
+ let range = area.get_pixel_range();
+ let (pw, ph) = (range.0.end - range.0.start, range.1.end - range.1.start);
+ let (xr, yr) = (chart.x_range(), chart.y_range());
+ let step = (
+ (xr.end - xr.start) / pw as f64,
+ (yr.end - yr.start) / ph as f64,
+ );
+
+ let mut min = f64::MAX;
+ let mut max = f64::MIN;
+ let vals: Vec<(f64, f64, f64)> = (0..pw * ph)
+ .into_par_iter()
+ .map(|i| {
+ let (x, y) = (
+ xr.start + step.0 * (i % pw) as f64,
+ yr.start + step.1 * (i / pw) as f64,
+ );
+ let mut optimizer = optimizer.clone();
+ let s = optimizer.model.get_settings_mut();
+ s.beta = x;
+ s.gamma = y;
+ let val = optimizer.objective_batch(&optimizer.model, ylist_true);
+ (x, y, val)
+ })
+ .collect();
+ vals.iter().for_each(|(_, _, val)| {
+ if *val > max {
+ max = *val;
+ }
+ if *val < min {
+ min = *val;
+ }
+ });
+ let ampl = 0.825 / (max - min);
+
+ for (x, y, c) in vals {
+ area.draw_pixel((x, y), &HSLColor((c - min) * ampl, 1.0, 0.5))
+ .unwrap();
+ }
+
+ if let Some(path_sto) = path_sto {
+ chart
+ .draw_series(std::iter::once(PathElement::new(
+ path_sto,
+ RGBColor(128, 128, 128),
+ )))
+ .unwrap()
+ .label("Stochastic")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RGBColor(128, 128, 128)));
+ }
+
+ if let Some(path_batch) = path_batch {
+ chart
+ .draw_series(std::iter::once(PathElement::new(path_batch, BLACK)))
+ .unwrap()
+ .label("Batch")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], BLACK));
+ }
+
+ chart
+ .configure_series_labels()
+ .background_style(WHITE.mix(0.8))
+ .border_style(BLACK)
+ .position(SeriesLabelPosition::UpperRight)
+ .draw()
+ .unwrap();
+}
+
+pub fn draw_comparison_chart(
+ filename: &str,
+ title: Option<&str>,
+ s: &model::sir::SirSettings<f64>,
+ xlist_explicit: &[Vect<f64, 2>],
+ xlist_implicit: &[Vect<f64, 2>],
+ xlist_true: &[Vect<f64, 2>],
+ dt_explicit: f64,
+ dt_implicit: f64,
+ dt_true: f64,
+) {
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, CHART_SIZE).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .y_label_area_size(30)
+ .x_label_area_size(30)
+ .build_cartesian_2d(0.0f64..xlist_true.len() as f64 * dt_true, 0.0f64..1.)
+ .unwrap();
+
+ chart.configure_mesh().x_desc("Time").draw().unwrap();
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_explicit
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_explicit, x[0])),
+ ShapeStyle::from(BLUE).stroke_width(3),
+ ))
+ .unwrap()
+ .label("Susceptible (explicit)")
+ .legend(|(x, y)| {
+ PathElement::new(
+ vec![(x, y), (x + 20, y)],
+ ShapeStyle::from(BLUE).stroke_width(3),
+ )
+ });
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_explicit
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_explicit, x[1])),
+ ShapeStyle::from(RED).stroke_width(3),
+ ))
+ .unwrap()
+ .label("Infected (explicit)")
+ .legend(|(x, y)| {
+ PathElement::new(
+ vec![(x, y), (x + 20, y)],
+ ShapeStyle::from(RED).stroke_width(3),
+ )
+ });
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_explicit
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_explicit, s.pop - x[0] - x[1])),
+ ShapeStyle::from(GREEN).stroke_width(3),
+ ))
+ .unwrap()
+ .label("Removed (explicit)")
+ .legend(|(x, y)| {
+ PathElement::new(
+ vec![(x, y), (x + 20, y)],
+ ShapeStyle::from(GREEN).stroke_width(3),
+ )
+ });
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_implicit
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_implicit, x[0])),
+ BLUE,
+ ))
+ .unwrap()
+ .label("Susceptible (implicit)")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], BLUE));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_implicit
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_implicit, x[1])),
+ RED,
+ ))
+ .unwrap()
+ .label("Infected (implicit)")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RED));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_implicit
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_implicit, s.pop - x[0] - x[1])),
+ GREEN,
+ ))
+ .unwrap()
+ .label("Removed (implicit)")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], GREEN));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_true
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_true, x[0])),
+ RGBColor(0, 0, 128),
+ ))
+ .unwrap()
+ .label("Susceptible (true)")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RGBColor(0, 0, 128)));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_true
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_true, x[1])),
+ RGBColor(128, 0, 0),
+ ))
+ .unwrap()
+ .label("Infected (true)")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RGBColor(128, 0, 0)));
+
+ chart
+ .draw_series(LineSeries::new(
+ xlist_true
+ .iter()
+ .enumerate()
+ .map(|(i, x)| (i as f64 * dt_true, s.pop - x[0] - x[1])),
+ RGBColor(0, 128, 0),
+ ))
+ .unwrap()
+ .label("Removed (true)")
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], RGBColor(0, 128, 0)));
+
+ chart
+ .configure_series_labels()
+ .background_style(WHITE.mix(0.8))
+ .border_style(BLACK)
+ .draw()
+ .unwrap();
+}
+
+pub fn draw_bike_chart(filename: &str, title: Option<&str>, xlists: &[(&str, &[f64])], dt: f64) {
+ let max_x = xlists
+ .iter()
+ .map(|(_, xlist)| xlist.len())
+ .max()
+ .expect("at least one series expected");
+ let max_y = *xlists
+ .iter()
+ .map(|(_, xlist)| {
+ xlist
+ .iter()
+ .max_by(|a, b| a.total_cmp(b))
+ .expect("at least one sample per series expected")
+ })
+ .max_by(|a, b| a.total_cmp(b))
+ .unwrap();
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, (640, 480)).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .margin_top(12)
+ .y_label_area_size(30)
+ .x_label_area_size(30)
+ .build_cartesian_2d(0.0f64..max_x as f64 * dt, 0.0f64..max_y)
+ .unwrap();
+
+ chart.configure_mesh().x_desc("Temps (s)").draw().unwrap();
+
+ for (list_i, (label, xlist)) in xlists.into_iter().enumerate() {
+ chart
+ .draw_series(LineSeries::new(
+ xlist.iter().enumerate().map(|(i, x)| (i as f64 * dt, *x)),
+ Palette100::pick(list_i + 1).stroke_width(2),
+ ))
+ .unwrap()
+ .label(*label)
+ .legend(move |(x, y)| {
+ PathElement::new(
+ vec![(x, y), (x + 20, y)],
+ Palette100::pick(list_i + 1).stroke_width(2),
+ )
+ });
+ }
+
+ chart
+ .configure_series_labels()
+ .border_style(BLACK)
+ .background_style(WHITE.mix(0.8))
+ .label_font(("Libertinus Serif", 20))
+ .draw()
+ .unwrap();
+}
+
+pub fn draw_bike_chart2(
+ filename: &str,
+ title: Option<&str>,
+ xlists1: &[(&str, &[f64])],
+ xlists2: &[(&str, &[f64])],
+ dt: f64,
+) {
+ let max_x = xlists1
+ .iter()
+ .chain(xlists2.iter())
+ .map(|(_, xlist)| xlist.len())
+ .max()
+ .expect("at least one series expected");
+ let max_y1 = *xlists1
+ .iter()
+ .map(|(_, xlist)| {
+ xlist
+ .iter()
+ .max_by(|a, b| a.total_cmp(b))
+ .expect("at least one sample per series expected")
+ })
+ .max_by(|a, b| a.total_cmp(b))
+ .unwrap();
+ let max_y2 = *xlists2
+ .iter()
+ .map(|(_, xlist)| {
+ xlist
+ .iter()
+ .max_by(|a, b| a.total_cmp(b))
+ .expect("at least one sample per series expected")
+ })
+ .max_by(|a, b| a.total_cmp(b))
+ .unwrap();
+ let filepath = format!("target/{}.png", filename);
+ let root = BitMapBackend::new(&filepath, (640, 480)).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::Name("cantarell"), 28.0, FontStyle::Normal),
+ );
+ }
+ let mut chart = chart
+ .margin_right(12)
+ .margin_top(12)
+ .y_label_area_size(50)
+ .x_label_area_size(30)
+ .right_y_label_area_size(50)
+ .build_cartesian_2d(0.0f64..max_x as f64 * dt, 0.0f64..max_y1)
+ .unwrap()
+ .set_secondary_coord(0.0f64..max_x as f64 * dt, 0.0f64..max_y2);
+
+ chart
+ .configure_mesh()
+ .x_desc("Temps (s)")
+ .y_desc("Vitesse (m/s)")
+ .draw()
+ .unwrap();
+ chart
+ .configure_secondary_axes()
+ .y_desc("Freinage")
+ .draw()
+ .unwrap();
+
+ for (list_i, (label, xlist)) in xlists1.into_iter().enumerate() {
+ chart
+ .draw_series(LineSeries::new(
+ xlist.iter().enumerate().map(|(i, x)| (i as f64 * dt, *x)),
+ Palette100::pick(list_i + 1).stroke_width(2),
+ ))
+ .unwrap()
+ .label(*label)
+ .legend(move |(x, y)| {
+ PathElement::new(
+ vec![(x, y), (x + 20, y)],
+ Palette100::pick(list_i + 1).stroke_width(2),
+ )
+ });
+ }
+
+ for (list_i, (label, xlist)) in (xlists1.len()..).zip(xlists2.into_iter()) {
+ chart
+ .draw_secondary_series(LineSeries::new(
+ xlist.iter().enumerate().map(|(i, x)| (i as f64 * dt, *x)),
+ Palette100::pick(list_i + 1).stroke_width(2),
+ ))
+ .unwrap()
+ .label(*label)
+ .legend(move |(x, y)| {
+ PathElement::new(
+ vec![(x, y), (x + 20, y)],
+ Palette100::pick(list_i + 1).stroke_width(2),
+ )
+ });
+ }
+
+ chart
+ .configure_series_labels()
+ .border_style(BLACK)
+ .background_style(WHITE.mix(0.8))
+ .label_font(("Libertinus Serif", 20))
+ .draw()
+ .unwrap();
+}
diff --git a/examples/old/live.rs b/examples/old/live.rs
new file mode 100644
index 0000000..91ec93e
--- /dev/null
+++ b/examples/old/live.rs
@@ -0,0 +1,91 @@
+use crate::{
+ model::{Coloring, Model, Settings},
+ solver::Solver,
+ space::Space,
+};
+
+use sdl2::{event::Event, keyboard::Keycode, pixels::PixelFormatEnum};
+use std::{
+ sync::{Arc, RwLock, RwLockReadGuard},
+ time::Duration,
+};
+
+const FRAMEDUR: u64 = 30;
+
+pub fn run<
+ T: Copy,
+ P,
+ S: Settings,
+ M: Model<T, S, D> + Coloring<T, P, [u8; 3], D>,
+ V: Solver<T, S, M, D>,
+ const D: usize,
+>(
+ space: Arc<RwLock<Space<T, S, M, V, D>>>,
+) {
+ let sdl_context = sdl2::init().unwrap();
+ let video_subsystem = sdl_context.video().unwrap();
+
+ let size = space.read().unwrap().size;
+
+ let window = video_subsystem
+ .window("Model", size.0 as u32, size.1 as u32)
+ .resizable()
+ .position_centered()
+ .opengl()
+ .build()
+ .unwrap();
+
+ let mut canvas = window.into_canvas().build().unwrap();
+ let texture_creator = canvas.texture_creator();
+
+ let mut event_pump = sdl_context.event_pump().unwrap();
+ let interval = Duration::from_millis(FRAMEDUR);
+
+ let mut image = vec![0u8; size.0 * size.1 * 3];
+
+ 'running: loop {
+ for event in event_pump.poll_iter() {
+ match event {
+ Event::Quit { .. }
+ | Event::KeyDown {
+ keycode: Some(Keycode::Escape),
+ ..
+ } => {
+ break 'running;
+ }
+ _ => {}
+ }
+ }
+
+ render(space.read().unwrap(), &mut image);
+
+ let mut texture = texture_creator
+ .create_texture_streaming(PixelFormatEnum::RGB24, size.0 as u32, size.1 as u32)
+ .unwrap();
+ texture.update(None, &image, size.0 * 3).unwrap();
+ canvas.copy(&texture, None, None).unwrap();
+
+ canvas.present();
+
+ std::thread::sleep(interval);
+ }
+}
+
+fn render<
+ T: Copy,
+ P,
+ S: Settings,
+ M: Model<T, S, D> + Coloring<T, P, [u8; 3], D>,
+ V: Solver<T, S, M, D>,
+ const D: usize,
+>(
+ space: RwLockReadGuard<Space<T, S, M, V, D>>,
+ image: &mut Vec<u8>,
+) {
+ image.resize(space.points.len() * 3, 0);
+ if let Some(pre) = M::prepare(space.points.iter().map(|point| point.pos)) {
+ for (i, point) in space.points.iter().enumerate() {
+ image[i * 3..i * 3 + 3].copy_from_slice(&M::color(&pre, point.pos));
+ }
+ }
+}
diff --git a/examples/old/main.rs b/examples/old/main.rs
new file mode 100644
index 0000000..5136deb
--- /dev/null
+++ b/examples/old/main.rs
@@ -0,0 +1,676 @@
+mod charts;
+mod live;
+mod model;
+mod opti;
+mod solver;
+mod space;
+mod utils;
+
+use model::Model;
+use opti::GradientDescentOptimizer;
+use solver::*;
+
+use nalgebra::vector;
+use rand::Rng;
+use std::{
+ collections::HashMap,
+ sync::{Arc, RwLock},
+ thread,
+ time::Duration,
+};
+
+fn main() {
+ bike();
+ //lyfe();
+ //stage();
+}
+
+#[allow(dead_code)]
+fn lyfe() {
+ let size = (800, 800);
+ let diffusion = vector![0.2, 0.2];
+
+ let model = model::constrained::Constrained {
+ s: model::constrained::ConstrainedSettings {
+ model: model::lyfe::Lyfe {
+ s: model::lyfe::LyfeSettings {
+ da: 0.2,
+ db: 0.3,
+ f: 0.03,
+ r: 0.061,
+ },
+ },
+ constraint: model::MinMaxConstraint {
+ min: [0.0, 0.0],
+ max: [1.0, 1.0],
+ },
+ _p: Default::default(),
+ },
+ };
+ /*let solver = ImplicitEulerSolver {
+ dt: 0.1,
+ tol: 0.000001,
+ niters: 100,
+ };*/
+ let solver = ExplicitEulerSolver { dt: 0.01 };
+
+ let mut space = space::Space {
+ model,
+ solver,
+ old_points: vec![
+ space::Point {
+ pos: vector![0.0, 0.0],
+ diffusion,
+ };
+ size.0 * size.1
+ ],
+ points: vec![
+ space::Point {
+ pos: vector![0.0, 0.0],
+ diffusion,
+ };
+ size.0 * size.1
+ ],
+ size,
+ sources: HashMap::new(),
+ time: 0.0,
+ _p: Default::default(),
+ };
+
+ let mut rng = rand::thread_rng();
+ for _ in 0..100 {
+ space.points[rng.gen_range(0..space.old_points.len())].pos[0] = 0.8;
+ }
+ //space.points[size.0 * size.1 / 2 + size.0 / 2].pos[0] = 0.1;
+ //space.points[size.0 * size.1 / 2 + size.0 / 2 + 100].pos[0] = 0.05;
+
+ let space = Arc::new(RwLock::new(space));
+
+ thread::spawn({
+ let space = space.clone();
+ let interval = Duration::from_millis(1);
+ move || loop {
+ space.write().unwrap().simulate(0.1);
+ std::thread::sleep(interval);
+ }
+ });
+
+ thread::spawn(move || live::run(space)).join().unwrap();
+}
+
+#[allow(dead_code)]
+fn giraffe() {
+ let size = (800, 800);
+ let diffusion = vector![0.1, 0.1];
+
+ let model = model::constrained::Constrained {
+ s: model::constrained::ConstrainedSettings {
+ model: model::giraffe::Giraffe {
+ s: model::giraffe::GiraffeSettings {
+ a_a: 0.7,
+ a_b: 0.2,
+ b_a: -0.5,
+ b_b: 0.1,
+ },
+ },
+ constraint: model::MinMaxConstraint {
+ min: [0.0, 0.0],
+ max: [1.0, 1.0],
+ },
+ _p: Default::default(),
+ },
+ };
+ /*let solver = ImplicitEulerSolver {
+ dt: 0.1,
+ tol: 0.000001,
+ niters: 100,
+ };*/
+ let solver = ExplicitEulerSolver { dt: 0.1 };
+
+ let mut space = space::Space {
+ model,
+ solver,
+ old_points: vec![
+ space::Point {
+ pos: vector![0.0, 0.0],
+ diffusion,
+ };
+ size.0 * size.1
+ ],
+ points: vec![
+ space::Point {
+ pos: vector![0.0, 0.0],
+ diffusion,
+ };
+ size.0 * size.1
+ ],
+ size,
+ sources: HashMap::new(),
+ time: 0.0,
+ _p: Default::default(),
+ };
+
+ let mut rng = rand::thread_rng();
+ for _ in 0..100 {
+ space.points[rng.gen_range(0..space.old_points.len())].pos[0] = 0.5;
+ }
+ //space.points[size.0 * size.1 / 2 + size.0 / 2].pos[0] = 0.1;
+ //space.points[size.0 * size.1 / 2 + size.0 / 2 + 100].pos[0] = 0.05;
+
+ let space = Arc::new(RwLock::new(space));
+
+ thread::spawn({
+ let space = space.clone();
+ let interval = Duration::from_millis(1);
+ move || loop {
+ space.write().unwrap().simulate(0.1);
+ std::thread::sleep(interval);
+ }
+ });
+
+ thread::spawn(move || live::run(space)).join().unwrap();
+}
+
+#[allow(dead_code)]
+fn stage() {
+ let mut rng = rand::thread_rng();
+
+ // ---- Initialization
+
+ let x0 = vector![0.99, 0.01];
+ let dt = 0.1;
+ let nsamples: usize = 400;
+ let nsamples_partial: usize = 40;
+
+ // ---- True data generation
+
+ let settings_true = model::sir::SirSettings {
+ beta: 0.6,
+ gamma: 0.1,
+ pop: 1.0,
+ };
+
+ let model = model::sir::Sir {
+ s: settings_true.clone(),
+ };
+ let solver = ImplicitEulerSolver {
+ dt: 0.1,
+ tol: 0.000001,
+ niters: 100,
+ };
+ let mut xlist_true = Vec::with_capacity(nsamples);
+ xlist_true.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = solver.f(&model, x);
+ xlist_true.push(x);
+ }
+
+ // ---- Calibration
+
+ let mut optimizer = GradientDescentOptimizer::new(model, solver);
+ let settings_random = model::sir::SirSettings {
+ beta: rng.gen(),
+ gamma: rng.gen(),
+ pop: 1.0,
+ };
+ *optimizer.model.get_settings_mut() = model::sir::SirSettings {
+ beta: 0.38960491052564317,
+ gamma: 0.6549130899826807,
+ pop: 1.0,
+ }; //settings_random.clone();
+ let mut optimizer_sto = optimizer.clone();
+
+ let mut xlist_random = Vec::with_capacity(nsamples);
+ xlist_random.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = optimizer.solver.f(&optimizer.model, x);
+ xlist_random.push(x);
+ }
+
+ // Batch
+
+ let mut path = Vec::new();
+ let mut error = Vec::new();
+ for rate in [1.0, 0.1, 0.01, 0.001] {
+ let (path_, error_) = &mut optimizer.calibrate_batch_record(
+ &xlist_true[..nsamples_partial],
+ 0.00001,
+ rate,
+ 1000,
+ 0..2,
+ );
+ path.append(path_);
+ error.append(error_);
+ }
+
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = optimizer.solver.f(&optimizer.model, x);
+ xlist.push(x);
+ }
+
+ // Stochastic
+
+ let mut path_sto = Vec::new();
+ let mut error_sto = Vec::new();
+ for rate in [1.0, 0.1, 0.01, 0.001] {
+ let (path_, error_) = &mut optimizer_sto.calibrate_stochastic_record(
+ &xlist_true[..nsamples_partial],
+ 0.00001,
+ rate,
+ 10,
+ 0..2,
+ );
+ path_sto.append(path_);
+ error_sto.append(error_);
+ }
+
+ let mut xlist_sto = Vec::with_capacity(nsamples);
+ xlist_sto.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = optimizer_sto.solver.f(&optimizer_sto.model, x);
+ xlist_sto.push(x);
+ }
+
+ // ---- Printing
+
+ println!("Random settings:\n{:?}", settings_random);
+ println!("Calibrated settings:\n{:?}", optimizer.model.get_settings());
+ println!("True settings:\n{:?}", settings_true);
+
+ // ---- Drawing
+
+ // Main plots
+
+ charts::draw_chart("sir_true", None, settings_true.pop, &xlist_true, dt);
+ charts::draw_chart("sir_random", None, settings_random.pop, &xlist_random, dt);
+ charts::draw_chart(
+ "sir_calibrated",
+ None,
+ optimizer.model.get_settings().pop,
+ &xlist,
+ dt,
+ );
+ charts::draw_error_chart2("error", None, &error, &error_sto);
+
+ charts::plot_objective(
+ "obj_partial",
+ None,
+ optimizer.clone(),
+ &xlist_true[..nsamples_partial],
+ Some(&path.iter().map(|v| (v[0], v[1])).collect::<Vec<_>>()),
+ Some(&path_sto.iter().map(|v| (v[0], v[1])).collect::<Vec<_>>()),
+ );
+ charts::plot_objective(
+ "obj",
+ None,
+ optimizer.clone(),
+ &xlist_true,
+ Some(&path.iter().map(|v| (v[0], v[1])).collect::<Vec<_>>()),
+ Some(&path_sto.iter().map(|v| (v[0], v[1])).collect::<Vec<_>>()),
+ );
+
+ // Implicit/explicit Euler comparison
+
+ {
+ let dur = 40f64;
+ let settings = model::sir::SirSettings {
+ beta: 0.999,
+ gamma: 0.5,
+ pop: 1.0,
+ };
+ let model = model::sir::Sir {
+ s: settings.clone(),
+ };
+ let solver_explicit = ExplicitEulerSolver { dt: 1.0 };
+ let solver_implicit = ImplicitEulerSolver {
+ dt: 1.0,
+ tol: 0.000001,
+ niters: 100,
+ };
+ let solver_true = ImplicitEulerSolver {
+ dt: 0.001,
+ tol: 0.000001,
+ niters: 100,
+ };
+ let nsamples_explicit = (dur / solver_explicit.dt) as usize;
+ let nsamples_implicit = (dur / solver_implicit.dt) as usize;
+ let nsamples_true = (dur / solver_true.dt) as usize;
+ let mut xlist_explicit = Vec::with_capacity(nsamples_explicit);
+ xlist_explicit.push(x0);
+ let mut xlist_implicit = Vec::with_capacity(nsamples_implicit);
+ xlist_implicit.push(x0);
+ let mut xlist_true = Vec::with_capacity(nsamples_true);
+ xlist_true.push(x0);
+
+ let mut x = x0;
+ for _ in 1..nsamples_explicit {
+ x = solver_explicit.f(&model, x);
+ xlist_explicit.push(x);
+ }
+
+ x = x0;
+ for _ in 1..nsamples_implicit {
+ x = solver_implicit.f(&model, x);
+ xlist_implicit.push(x);
+ }
+
+ x = x0;
+ for _ in 1..nsamples_true {
+ x = solver_true.f(&model, x);
+ xlist_true.push(x);
+ }
+
+ charts::draw_comparison_chart(
+ "comp_euler",
+ None,
+ &settings,
+ &xlist_explicit,
+ &xlist_implicit,
+ &xlist_true,
+ solver_explicit.dt,
+ solver_implicit.dt,
+ solver_true.dt,
+ );
+ }
+
+ // SIRV charts
+
+ {
+ let nsamples = 1000;
+ let settings = model::sirv::SirvSettings {
+ beta: 0.8,
+ gamma: 0.2,
+ lambda: 0.025,
+ mu: 0.02,
+ pop: 1.0,
+ };
+ let model = model::sirv::Sirv { s: settings };
+
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+
+ let mut x = x0;
+ for _ in 1..nsamples {
+ x = optimizer.solver.f(&model, x);
+ xlist.push(x);
+ }
+
+ charts::draw_chart("sirv", None, model.get_settings().pop, &xlist, dt);
+ }
+}
+
+#[allow(dead_code)]
+fn bike() {
+ let mut rng = rand::thread_rng();
+
+ // ---- Initialization
+
+ let x0 = vector![0.0, 60. / 3.6];
+ let dt = 0.1;
+ let nsamples: usize = 2000;
+
+ // ---- Data generation
+
+ let settings_true = model::bike::BikeSettings::<f64> {
+ cx: 0.25,
+ g: 9.81,
+ m: 70.0,
+ th: 0.11,
+ };
+ println!(
+ "true: A={} ; B={}",
+ -settings_true.cx / settings_true.m,
+ settings_true.g * settings_true.th.sin() - 80. / settings_true.m
+ );
+
+ let model = model::bike::Bike {
+ s: settings_true.clone(),
+ };
+ /*let solver = ImplicitEulerSolver {
+ dt,
+ tol: 0.000001,
+ niters: 100,
+ };*/
+ let solver = ExplicitEulerSolver { dt };
+ let mut xlist_true = Vec::with_capacity(nsamples);
+ xlist_true.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = solver.f(&model, x);
+ if x[1] < 0. {
+ x[1] = 0.;
+ }
+ xlist_true.push(x);
+ }
+
+ // -- Alternative settings
+
+ // Greater theta
+
+ let mut settings_greater_th = settings_true.clone();
+ settings_greater_th.th = 0.12;
+ println!(
+ "gtth: A={} ; B={}",
+ -settings_greater_th.cx / settings_greater_th.m,
+ settings_greater_th.g * settings_greater_th.th.sin() - 80. / settings_greater_th.m
+ );
+
+ let xlist_greater_th = {
+ let model = model::bike::Bike {
+ s: settings_greater_th.clone(),
+ };
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = solver.f(&model, x);
+ if x[1] < 0. {
+ x[1] = 0.;
+ }
+ xlist.push(x);
+ }
+ xlist
+ };
+
+ // Optimal braking
+
+ let mut settings_opti1 = model::bike2::BikeSettings {
+ cx: settings_true.cx,
+ g: settings_true.g,
+ m: settings_true.m,
+ th: std::f64::consts::PI / 180. * 15.,
+ b: |x, v, s| {
+ let mu = 0.1;
+ let gx = 0.65;
+ let gy = 1.05;
+ let magic = 1.0;
+ (
+ (
+ s.m * s.g * (mu * s.th.cos() + s.th.sin()) - s.cx * v * v,
+ 0.,
+ ),
+ (-2. * s.cx * v, 0.),
+ )
+ },
+ };
+
+ let (xlist_opti1, blist_opti1) = {
+ let model = model::bike2::Bike {
+ s: settings_opti1.clone(),
+ };
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+ let mut blist = Vec::with_capacity(nsamples);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ blist.push((settings_opti1.b)(x[0], x[1], &settings_opti1).0 .0);
+ x = solver.f(&model, x);
+ if x[1] < 0. {
+ x[1] = 0.;
+ break;
+ }
+ xlist.push(x);
+ }
+ (xlist, blist)
+ };
+
+ let mut settings_opti2 = model::bike2::BikeSettings {
+ cx: settings_true.cx,
+ g: settings_true.g,
+ m: settings_true.m,
+ th: std::f64::consts::PI / 180. * 15.,
+ b: |x, v, s| {
+ let mu = 0.1;
+ let gx = 0.65;
+ let gy = 1.05;
+ let magic = 1.0;
+ (
+ (
+ s.m * s.g * (mu * s.th.cos() + s.th.sin()) - s.cx * v * v,
+ s.m * s.g * s.th.cos() * (mu + gx / gy),
+ ),
+ (-2. * s.cx * v, 0.),
+ )
+ },
+ };
+
+ let xlist_opti2 = {
+ let model = model::bike2::Bike {
+ s: settings_opti2.clone(),
+ };
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+ let mut x = x0;
+ for _ in 0..nsamples - 1 {
+ x = solver.f(&model, x);
+ if x[1] < 0. {
+ x[1] = 0.;
+ break;
+ }
+ xlist.push(x);
+ }
+ xlist
+ };
+
+ // -- ODE solution
+
+ let xlist_ode = {
+ let settings = settings_true.clone();
+ let a = -settings.cx / settings.m;
+ let b = settings.g * settings.th.sin() - 80. / settings.m;
+ /*let r = (-a*b).sqrt();
+ assert!(!r.is_nan());
+ let alpha = (1.+x0[1]*a/r)/(1.-x0[1]*a/r);*/
+ let r = (a * b).sqrt();
+ assert!(!r.is_nan());
+ let alpha = -r / (a * x0[1]);
+ println!("alpha: {alpha}");
+ let stop = (1. / alpha).atan() / r;
+ println!("Stop: {stop}");
+
+ let bmax =
+ settings.m * settings.g * (settings.th.cos() + settings.th.sin()) - settings.cx * 27.;
+ println!("bmax: {bmax}");
+
+ let mut xlist = Vec::with_capacity(nsamples);
+ xlist.push(x0);
+ let mut x = x0;
+ for t in 0..nsamples - 1 {
+ let t = t as f64 * dt;
+ //dbg!((beta*c*(-t*c).exp()-alpha*c*(t*c).exp()));
+ //dbg!((alpha*(t*c).exp()+beta*(-t*c).exp()));
+ //let v = (beta*c*(-t*c).exp()-alpha*c*(t*c).exp())/a/(alpha*(t*c).exp()+beta*(-t*c).exp());
+ //let v = r/a*(alpha*(-t*r).exp()-(t*r).exp())/(alpha*(-t*r).exp()+(t*r).exp());
+ let mut v = r / a * (alpha * (t * r).sin() - (t * r).cos())
+ / (alpha * (t * r).cos() + (t * r).sin());
+ if v.is_nan() {
+ panic!("NaN");
+ }
+ v = v.max(0.).min(100.);
+ x = vector![0., v];
+ xlist.push(x);
+ if t > stop {
+ break;
+ }
+ }
+ xlist
+ };
+
+ // ---- Drawing
+
+ // Main plots
+
+ charts::draw_bike_chart(
+ "bike_x",
+ None,
+ &[(
+ "x (m)",
+ &xlist_true.iter().map(|x| x[0]).collect::<Vec<_>>(),
+ )],
+ dt,
+ );
+ charts::draw_bike_chart(
+ "bike_v",
+ None,
+ &[(
+ "v (m/s)",
+ &xlist_true.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ )],
+ dt,
+ );
+ charts::draw_bike_chart2(
+ "bike_opti1",
+ None,
+ &[(
+ "v (m/s)",
+ &xlist_opti1.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ )],
+ &[("b", &blist_opti1)],
+ dt,
+ );
+ charts::draw_bike_chart(
+ "bike_opti2",
+ None,
+ &[
+ (
+ "v (m/s) (arrière)",
+ &xlist_opti1.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ ),
+ (
+ "v (m/s) (arrière+avant)",
+ &xlist_opti2.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ ),
+ ],
+ dt,
+ );
+ charts::draw_bike_chart(
+ "bike_var_theta",
+ None,
+ &[
+ (
+ &format!("v (θ={}, B<0)", settings_true.th),
+ &xlist_true.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ ),
+ (
+ &format!("v (θ={}, B>0)", settings_greater_th.th),
+ &xlist_greater_th.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ ),
+ ],
+ dt,
+ );
+ charts::draw_bike_chart(
+ "bike_ode_v",
+ None,
+ &[(
+ "v (m/s)",
+ &xlist_ode.iter().map(|x| x[1]).collect::<Vec<_>>(),
+ )],
+ dt,
+ );
+}
diff --git a/examples/old/model.rs b/examples/old/model.rs
new file mode 100644
index 0000000..0472e24
--- /dev/null
+++ b/examples/old/model.rs
@@ -0,0 +1,595 @@
+use crate::utils::*;
+
+use nalgebra::{base::*, matrix, vector};
+use std::marker::PhantomData;
+
+pub trait Settings {}
+
+pub trait Model<T, S: Settings, const D: usize> {
+ /// Returns f(x)
+ fn f(&self, x: Vect<T, D>) -> Vect<T, D>;
+ /// Returns df(x)/dx
+ fn df(&self, x: Vect<T, D>) -> Mat<T, D, D>;
+ fn get_settings(&self) -> &S;
+ fn get_settings_mut(&mut self) -> &mut S;
+}
+
+pub trait Coloring<T, P, C, const D: usize> {
+ fn prepare<I: Iterator<Item = Vect<T, D>>>(iter: I) -> Option<P>;
+ fn color(pre: &P, val: Vect<T, D>) -> C;
+}
+
+pub struct MinMaxColoringInfo<T, const D: usize> {
+ pub min: [T; D],
+ pub max: [T; D],
+}
+
+pub trait Constraint<T, const D: usize> {
+ fn constrain_mut(&self, val: &mut Vect<T, D>);
+ fn constrain(&self, mut val: Vect<T, D>) -> Vect<T, D> {
+ self.constrain_mut(&mut val);
+ val
+ }
+}
+
+pub struct MinMaxConstraint<T, const D: usize> {
+ pub min: [T; D],
+ pub max: [T; D],
+}
+
+impl<T: Copy + PartialOrd, const D: usize> Constraint<T, D> for MinMaxConstraint<T, D> {
+ fn constrain_mut(&self, val: &mut Vect<T, D>) {
+ for ((comp, min), max) in val.iter_mut().zip(self.min.iter()).zip(self.max.iter()) {
+ *comp = num_traits::clamp(*comp, *min, *max)
+ }
+ }
+}
+
+// Models
+
+pub mod constrained {
+ use super::*;
+
+ #[derive(Clone, Debug)]
+ pub struct ConstrainedSettings<T, O, M, S> {
+ pub constraint: O,
+ pub model: M,
+ pub _p: PhantomData<(T, S)>,
+ }
+
+ impl<T, O, M, S> Settings for ConstrainedSettings<T, O, M, S> {}
+
+ #[derive(Clone)]
+ pub struct Constrained<T, O, M, S> {
+ pub s: ConstrainedSettings<T, O, M, S>,
+ }
+
+ impl<T, O: Constraint<T, D>, M: Model<T, S, D>, S: Settings, const D: usize>
+ Model<T, ConstrainedSettings<T, O, M, S>, D> for Constrained<T, O, M, S>
+ {
+ fn f(&self, x: Vect<T, D>) -> Vect<T, D> {
+ self.s.constraint.constrain(self.s.model.f(x))
+ }
+ fn df(&self, x: Vect<T, D>) -> Mat<T, D, D> {
+ self.s.model.df(x)
+ }
+ fn get_settings(&self) -> &ConstrainedSettings<T, O, M, S> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut ConstrainedSettings<T, O, M, S> {
+ &mut self.s
+ }
+ }
+
+ impl<T, O, M: Coloring<T, P, C, D>, S, P, C, const D: usize> Coloring<T, P, C, D>
+ for Constrained<T, O, M, S>
+ {
+ fn prepare<I: Iterator<Item = Vect<T, D>>>(iter: I) -> Option<P> {
+ M::prepare(iter)
+ }
+ fn color(pre: &P, val: Vect<T, D>) -> C {
+ M::color(pre, val)
+ }
+ }
+}
+
+/// SIR model without vital dynamics
+pub mod sir {
+ use super::*;
+
+ #[derive(Clone, Debug)]
+ pub struct SirSettings<T> {
+ /// Transmission probability
+ pub beta: T,
+ /// Removal probability
+ pub gamma: T,
+ pub pop: T,
+ }
+
+ impl<T> Settings for SirSettings<T> {}
+
+ #[derive(Clone)]
+ pub struct Sir<T> {
+ pub s: SirSettings<T>,
+ }
+
+ impl Model<f64, SirSettings<f64>, 2> for Sir<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ vector![
+ -self.s.beta * x[0] * x[1] / self.s.pop,
+ self.s.beta * x[0] * x[1] / self.s.pop - self.s.gamma * x[1]
+ ]
+ }
+ fn df(&self, x: Vector2<f64>) -> Matrix2<f64> {
+ matrix![
+ -self.s.beta*x[1]/self.s.pop, -self.s.beta*x[0]/self.s.pop;
+ self.s.beta*x[1]/self.s.pop, self.s.beta*x[0]/self.s.pop-self.s.gamma
+ ]
+ }
+ fn get_settings(&self) -> &SirSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut SirSettings<f64> {
+ &mut self.s
+ }
+ }
+
+ impl<T> From<SirSettings<T>> for Vect<T, 3> {
+ fn from(s: SirSettings<T>) -> Self {
+ vector![s.beta, s.gamma, s.pop]
+ }
+ }
+
+ impl<T: Scalar + Copy> From<Vect<T, 3>> for SirSettings<T> {
+ fn from(v: Vect<T, 3>) -> Self {
+ Self {
+ beta: v[0],
+ gamma: v[1],
+ pop: v[2],
+ }
+ }
+ }
+}
+
+/// SIR model with vital dynamics, constant population
+pub mod sirv {
+ use super::*;
+ #[derive(Clone, Debug)]
+ pub struct SirvSettings<T> {
+ /// Transmission probability
+ pub beta: T,
+ /// Removal probability
+ pub gamma: T,
+ /// Birth rate
+ pub lambda: T,
+ /// Death rate
+ pub mu: T,
+ pub pop: T,
+ }
+
+ impl<T> Settings for SirvSettings<T> {}
+
+ #[derive(Clone)]
+ pub struct Sirv<T> {
+ pub s: SirvSettings<T>,
+ }
+
+ impl Model<f64, SirvSettings<f64>, 2> for Sirv<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ vector![
+ self.s.lambda - self.s.beta * x[0] * x[1] / self.s.pop - self.s.mu * x[0],
+ self.s.beta * x[0] * x[1] / self.s.pop - self.s.gamma * x[1] - self.s.mu * x[1]
+ ]
+ }
+ fn df(&self, x: Vector2<f64>) -> Matrix2<f64> {
+ matrix![
+ -self.s.beta*x[1]/self.s.pop - self.s.mu, -self.s.beta*x[0]/self.s.pop;
+ self.s.beta*x[1]/self.s.pop, self.s.beta*x[0]/self.s.pop-self.s.gamma - self.s.mu
+ ]
+ }
+ fn get_settings(&self) -> &SirvSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut SirvSettings<f64> {
+ &mut self.s
+ }
+ }
+
+ impl<T> From<SirvSettings<T>> for Vect<T, 5> {
+ fn from(s: SirvSettings<T>) -> Self {
+ vector![s.beta, s.gamma, s.lambda, s.mu, s.pop]
+ }
+ }
+
+ impl<T: Scalar + Copy> From<Vect<T, 5>> for SirvSettings<T> {
+ fn from(v: Vect<T, 5>) -> Self {
+ Self {
+ beta: v[0],
+ gamma: v[1],
+ lambda: v[2],
+ mu: v[3],
+ pop: v[4],
+ }
+ }
+ }
+}
+
+/// Giraffe
+pub mod giraffe {
+ use super::*;
+
+ #[derive(Clone, Debug)]
+ pub struct GiraffeSettings<T> {
+ pub a_a: T,
+ pub a_b: T,
+ pub b_a: T,
+ pub b_b: T,
+ }
+
+ impl<T> Settings for GiraffeSettings<T> {}
+
+ #[derive(Clone)]
+ pub struct Giraffe<T> {
+ pub s: GiraffeSettings<T>,
+ }
+
+ impl Model<f64, GiraffeSettings<f64>, 2> for Giraffe<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ vector![
+ self.s.a_a * x[0] + self.s.b_a * x[1],
+ self.s.a_b * x[0] + self.s.b_b * x[1]
+ ]
+ }
+ fn df(&self, _x: Vector2<f64>) -> Matrix2<f64> {
+ matrix![
+ self.s.a_a, self.s.b_a;
+ self.s.a_b, self.s.b_b
+ ]
+ }
+ fn get_settings(&self) -> &GiraffeSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut GiraffeSettings<f64> {
+ &mut self.s
+ }
+ }
+
+ impl<T> From<GiraffeSettings<T>> for Vect<T, 4> {
+ fn from(s: GiraffeSettings<T>) -> Self {
+ vector![s.a_a, s.a_b, s.b_a, s.b_b]
+ }
+ }
+
+ impl<T: Scalar + Copy> From<Vect<T, 4>> for GiraffeSettings<T> {
+ fn from(v: Vect<T, 4>) -> Self {
+ Self {
+ a_a: v[0],
+ a_b: v[1],
+ b_a: v[2],
+ b_b: v[3],
+ }
+ }
+ }
+
+ /*impl Coloring<f64, MinMaxColoringInfo<f64, 2>, [u8; 3], 2> for Giraffe<f64> {
+ fn prepare<I: Iterator<Item = Vect<f64, 2>>>(
+ iter: I,
+ ) -> Option<MinMaxColoringInfo<f64, 2>> {
+ let mut r = MinMaxColoringInfo {
+ min: [1.0, 1.0],
+ max: [0.0, 0.0],
+ };
+ for val in iter {
+ if val[0] < r.min[0] {
+ r.min[0] = val[0];
+ }
+ if val[0] > r.max[0] {
+ r.max[0] = val[0];
+ }
+ if val[1] < r.min[1] {
+ r.min[1] = val[1];
+ }
+ if val[1] > r.max[1] {
+ r.max[1] = val[1];
+ }
+ }
+ if r.min[0] == r.max[0] || r.min[1] == r.max[1] {
+ None
+ } else {
+ Some(r)
+ }
+ }
+ fn color(pre: &MinMaxColoringInfo<f64, 2>, val: Vect<f64, 2>) -> [u8; 3] {
+ [
+ ((val[1] - pre.min[1]).abs() / (pre.max[1] - pre.min[1]).abs() * 255.0) as u8,
+ 0,
+ ((val[0] - pre.min[0]).abs() / (pre.max[0] - pre.min[0]).abs() * 255.0) as u8,
+ ]
+ }*/
+
+ impl Coloring<f64, (), [u8; 3], 2> for Giraffe<f64> {
+ fn prepare<I: Iterator<Item = Vect<f64, 2>>>(_iter: I) -> Option<()> {
+ Some(())
+ }
+ fn color(_pre: &(), val: Vect<f64, 2>) -> [u8; 3] {
+ [
+ (val[0] * 255.0) as u8,
+ ((val[0] + val[1]) * 127.5) as u8,
+ 255 - (val[1] * 255.0) as u8,
+ ]
+ }
+ }
+}
+
+/// https://arxiv.org/abs/2210.05227
+pub mod lyfe {
+ use super::*;
+
+ #[derive(Clone, Debug)]
+ pub struct LyfeSettings<T> {
+ pub da: T,
+ pub db: T,
+ pub f: T,
+ pub r: T,
+ }
+
+ impl<T> Settings for LyfeSettings<T> {}
+
+ #[derive(Clone)]
+ pub struct Lyfe<T> {
+ pub s: LyfeSettings<T>,
+ }
+
+ impl Model<f64, LyfeSettings<f64>, 2> for Lyfe<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ vector![
+ x[0] * x[1].powi(2) + self.s.f * (1.0 - x[0]),
+ x[0] * x[1].powi(2) - (self.s.f + self.s.r) * x[1]
+ ]
+ }
+ fn df(&self, x: Vector2<f64>) -> Matrix2<f64> {
+ matrix![
+ x[1].powi(2)-self.s.f, 2.0*x[0]*x[1];
+ x[1].powi(2), 2.0*x[0]*x[1] - (self.s.f+self.s.r)
+ ]
+ }
+ fn get_settings(&self) -> &LyfeSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut LyfeSettings<f64> {
+ &mut self.s
+ }
+ }
+
+ impl<T> From<LyfeSettings<T>> for Vect<T, 4> {
+ fn from(s: LyfeSettings<T>) -> Self {
+ vector![s.da, s.db, s.f, s.r]
+ }
+ }
+
+ impl<T: Scalar + Copy> From<Vect<T, 4>> for LyfeSettings<T> {
+ fn from(v: Vect<T, 4>) -> Self {
+ Self {
+ da: v[0],
+ db: v[1],
+ f: v[2],
+ r: v[3],
+ }
+ }
+ }
+
+ /*impl Coloring<f64, MinMaxColoringInfo<f64, 2>, [u8; 3], 2> for Lyfe<f64> {
+ fn prepare<I: Iterator<Item = Vect<f64, 2>>>(
+ iter: I,
+ ) -> Option<MinMaxColoringInfo<f64, 2>> {
+ let mut r = MinMaxColoringInfo {
+ min: [1.0, 1.0],
+ max: [0.0, 0.0],
+ };
+ for val in iter {
+ if val[0] < r.min[0] {
+ r.min[0] = val[0];
+ }
+ if val[0] > r.max[0] {
+ r.max[0] = val[0];
+ }
+ if val[1] < r.min[1] {
+ r.min[1] = val[1];
+ }
+ if val[1] > r.max[1] {
+ r.max[1] = val[1];
+ }
+ }
+ if r.min[0] == r.max[0] || r.min[1] == r.max[1] {
+ None
+ } else {
+ Some(r)
+ }
+ }
+ fn color(pre: &MinMaxColoringInfo<f64, 2>, val: Vect<f64, 2>) -> [u8; 3] {
+ [
+ ((val[1] - pre.min[1]).abs() / (pre.max[1] - pre.min[1]).abs() * 255.0) as u8,
+ 0,
+ ((val[0] - pre.min[0]).abs() / (pre.max[0] - pre.min[0]).abs() * 255.0) as u8,
+ ]
+ }*/
+
+ impl Coloring<f64, (), [u8; 3], 2> for Lyfe<f64> {
+ fn prepare<I: Iterator<Item = Vect<f64, 2>>>(_iter: I) -> Option<()> {
+ Some(())
+ }
+ fn color(_pre: &(), val: Vect<f64, 2>) -> [u8; 3] {
+ [
+ (val[0] * 255.0) as u8,
+ ((val[0] + val[1]) * 127.5) as u8,
+ 255 - (val[1] * 255.0) as u8,
+ ]
+ }
+ }
+}
+
+pub mod bike {
+ use super::*;
+
+ #[derive(Clone, Debug)]
+ pub struct BikeSettings<T> {
+ /// Pénétration dans l'air
+ pub cx: T,
+ /// Gravité
+ pub g: T,
+ /// Masse
+ pub m: T,
+ /// Pente
+ pub th: T,
+ }
+
+ impl<T> Settings for BikeSettings<T> {}
+
+ fn b(x: f64, v: f64) -> f64 {
+ 80.
+ }
+
+ fn db(x: f64, v: f64) -> Vector2<f64> {
+ vector![0., 0.]
+ }
+
+ #[derive(Clone)]
+ pub struct Bike<T> {
+ pub s: BikeSettings<T>,
+ }
+
+ impl Model<f64, BikeSettings<f64>, 2> for Bike<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ vector![
+ x[1],
+ self.s.g * self.s.th.sin() - (self.s.cx * x[1].powi(2) + b(x[0], x[1])) / self.s.m
+ ]
+ }
+ fn df(&self, x: Vector2<f64>) -> Matrix2<f64> {
+ let dbx = db(x[0], x[1]);
+ matrix![
+ 0., 1.;
+ -dbx[0]/self.s.m, -(2.*self.s.cx*x[1]+dbx[1])/self.s.m
+ ]
+ }
+ fn get_settings(&self) -> &BikeSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut BikeSettings<f64> {
+ &mut self.s
+ }
+ }
+
+ impl<T> From<BikeSettings<T>> for Vect<T, 4> {
+ fn from(s: BikeSettings<T>) -> Self {
+ vector![s.cx, s.g, s.m, s.th]
+ }
+ }
+
+ impl<T: Scalar + Copy> From<Vect<T, 4>> for BikeSettings<T> {
+ fn from(v: Vect<T, 4>) -> Self {
+ Self {
+ cx: v[0],
+ g: v[1],
+ m: v[2],
+ th: v[3],
+ }
+ }
+ }
+
+ /*impl Coloring<f64, MinMaxColoringInfo<f64, 2>, [u8; 3], 2> for Lyfe<f64> {
+ fn prepare<I: Iterator<Item = Vect<f64, 2>>>(
+ iter: I,
+ ) -> Option<MinMaxColoringInfo<f64, 2>> {
+ let mut r = MinMaxColoringInfo {
+ min: [1.0, 1.0],
+ max: [0.0, 0.0],
+ };
+ for val in iter {
+ if val[0] < r.min[0] {
+ r.min[0] = val[0];
+ }
+ if val[0] > r.max[0] {
+ r.max[0] = val[0];
+ }
+ if val[1] < r.min[1] {
+ r.min[1] = val[1];
+ }
+ if val[1] > r.max[1] {
+ r.max[1] = val[1];
+ }
+ }
+ if r.min[0] == r.max[0] || r.min[1] == r.max[1] {
+ None
+ } else {
+ Some(r)
+ }
+ }
+ fn color(pre: &MinMaxColoringInfo<f64, 2>, val: Vect<f64, 2>) -> [u8; 3] {
+ [
+ ((val[1] - pre.min[1]).abs() / (pre.max[1] - pre.min[1]).abs() * 255.0) as u8,
+ 0,
+ ((val[0] - pre.min[0]).abs() / (pre.max[0] - pre.min[0]).abs() * 255.0) as u8,
+ ]
+ }*/
+
+ impl Coloring<f64, (), [u8; 3], 2> for Bike<f64> {
+ fn prepare<I: Iterator<Item = Vect<f64, 2>>>(_iter: I) -> Option<()> {
+ Some(())
+ }
+ fn color(_pre: &(), val: Vect<f64, 2>) -> [u8; 3] {
+ [
+ (val[0] * 255.0) as u8,
+ ((val[0] + val[1]) * 127.5) as u8,
+ 255 - (val[1] * 255.0) as u8,
+ ]
+ }
+ }
+}
+
+pub mod bike2 {
+ use super::*;
+
+ #[derive(Clone, Debug)]
+ pub struct BikeSettings<T> {
+ /// Pénétration dans l'air
+ pub cx: T,
+ /// Gravité
+ pub g: T,
+ /// Masse
+ pub m: T,
+ /// Pente
+ pub th: T,
+ /// (x, v, settings) -> ((b1, b2), (b1', b2'))
+ pub b: fn(T, T, &BikeSettings<T>) -> ((T, T), (T, T)),
+ }
+ impl<T> Settings for BikeSettings<T> {}
+
+ #[derive(Clone)]
+ pub struct Bike<T> {
+ pub s: BikeSettings<T>,
+ }
+
+ impl Model<f64, BikeSettings<f64>, 2> for Bike<f64> {
+ fn f(&self, x: Vector2<f64>) -> Vector2<f64> {
+ let ((b1, b2), (_db1, _db2)) = (self.s.b)(x[0], x[1], &self.s);
+ vector![
+ x[1],
+ self.s.g * self.s.th.sin() - (self.s.cx * x[1].powi(2) + b1 + b2) / self.s.m
+ ]
+ }
+ fn df(&self, x: Vector2<f64>) -> Matrix2<f64> {
+ let ((_b1, _b2), (db1, db2)) = (self.s.b)(x[0], x[1], &self.s);
+ matrix![
+ 0., 1.;
+ 0., -(2.*self.s.cx*x[1]+db1+db2)/self.s.m
+ ]
+ }
+ fn get_settings(&self) -> &BikeSettings<f64> {
+ &self.s
+ }
+ fn get_settings_mut(&mut self) -> &mut BikeSettings<f64> {
+ &mut self.s
+ }
+ }
+}
diff --git a/examples/old/opti.rs b/examples/old/opti.rs
new file mode 100644
index 0000000..d47f83c
--- /dev/null
+++ b/examples/old/opti.rs
@@ -0,0 +1,399 @@
+use crate::{model::*, solver::*, utils::*};
+
+use nalgebra::{base::*, ComplexField};
+use num_traits::{FromPrimitive, Zero};
+use std::marker::PhantomData;
+
+pub fn newton<
+ T: Copy + Scalar + ComplexField<RealField = T> + PartialOrd,
+ S: Settings,
+ const D: usize,
+>(
+ model: &impl Model<T, S, D>,
+ x0: Vect<T, D>,
+ dt: T,
+ tol: T,
+ niters: usize,
+) -> Vect<T, D>
+where
+ Const<D>: ToTypenum + DimMin<Const<D>, Output = Const<D>>,
+{
+ let mut x = x0;
+
+ for _ in 0..niters {
+ if let Some(m) = (Mat::<T, D, D>::identity() - model.df(x) * dt).try_inverse() {
+ let dx = m * (x - x0 - model.f(x) * dt);
+ if dx.norm() < tol {
+ break;
+ }
+ x -= dx;
+ } else {
+ break;
+ }
+ }
+ x
+}
+
+/// Slower version using a linear system.
+pub fn _newton_slow<
+ T: Copy + Scalar + ComplexField<RealField = T> + PartialOrd,
+ S: Settings,
+ const D: usize,
+>(
+ model: &impl Model<T, S, D>,
+ x0: Vect<T, D>,
+ dt: T,
+ tol: T,
+ niters: usize,
+) -> Vect<T, D>
+where
+ Const<D>: ToTypenum + DimMin<Const<D>, Output = Const<D>>,
+{
+ let mut x = x0;
+
+ for _ in 0..niters {
+ let fi = model.f(x);
+ let dfi = model.df(x);
+ let g = x0 + fi * dt - x;
+ let dgdx = dfi * dt - Mat::<T, D, D>::identity();
+ if let Some(dx) = dgdx.lu().solve(&g) {
+ if dx.norm() < tol {
+ break;
+ }
+ x -= dx;
+ } else {
+ break;
+ }
+ }
+ x
+}
+
+#[derive(Clone)]
+pub struct GradientDescentOptimizer<
+ T,
+ S: Settings,
+ M: Model<T, S, D> + Clone,
+ R: Solver<T, S, M, D>,
+ const D: usize,
+ const DS: usize,
+> {
+ pub model: M,
+ pub solver: R,
+ _ph: PhantomData<(T, S)>,
+}
+
+impl<
+ T: Copy + ComplexField<RealField = T> + FromPrimitive,
+ S: Settings + Clone + Into<Vect<T, DS>> + From<Vect<T, DS>>,
+ M: Model<T, S, D> + Clone,
+ R: Solver<T, S, M, D>,
+ const D: usize,
+ const DS: usize,
+ > GradientDescentOptimizer<T, S, M, R, D, DS>
+where
+ Vect<T, DS>: std::ops::DivAssign<T> + std::ops::Mul<T, Output = Vect<T, DS>>,
+{
+ pub fn new(model: M, solver: R) -> Self {
+ Self {
+ model,
+ solver,
+ _ph: PhantomData {},
+ }
+ }
+
+ /// Distance between f(x) and y_true, that we want to minimize
+ pub fn objective(&self, model: &M, x: Vect<T, D>, y_true: Vect<T, D>) -> T {
+ (self.solver.f(model, x) - y_true).norm()
+ }
+
+ /// Return gradient of the objective function
+ /// (opposite direction for Settings in order to make f(x) closer to y_true)
+ ///
+ /// `free_settings` gives the indices of the settings we can change.
+ /// For example, if all the settings can change, set it to `0..DS`.
+ pub fn objective_gradient(
+ &self,
+ x: Vect<T, D>,
+ y_true: Vect<T, D>,
+ ep: T,
+ free_settings: impl Iterator<Item = usize>,
+ ) -> Vect<T, DS> {
+ let diff = self.objective(&self.model, x, y_true);
+ let mut model = self.model.clone();
+ let s: Vect<T, DS> = model.get_settings().clone().into();
+ let mut si = s;
+ let mut ds = Vect::<T, DS>::zero();
+ for i in free_settings {
+ si[i] += ep;
+ *model.get_settings_mut() = si.into();
+ ds[i] = (self.objective(&model, x, y_true) - diff) / ep;
+ si[i] = s[i];
+ }
+ ds
+ }
+
+ pub fn objective_gradient_batch(
+ &self,
+ ylist_true: &[Vect<T, D>],
+ ep: T,
+ free_settings: impl Iterator<Item = usize> + Clone,
+ ) -> Vect<T, DS> {
+ let nsamples = T::from_usize(ylist_true.len() - 1).unwrap();
+
+ let mut ds_batch = Vect::<T, DS>::zero();
+ for (x, y_true) in ylist_true.iter().zip(ylist_true.iter().skip(1)) {
+ let ds = self.objective_gradient(*x, *y_true, ep, free_settings.clone());
+ ds_batch += ds;
+ }
+ ds_batch / nsamples
+ }
+
+ /// Calibrate settings using batch GD
+ pub fn calibrate_batch(
+ &mut self,
+ ylist_true: &[Vect<T, D>],
+ ep: T,
+ rate: T,
+ niters: usize,
+ free_settings: impl Iterator<Item = usize> + Clone,
+ ) {
+ for _ in 0..niters {
+ let ds_batch = self.objective_gradient_batch(ylist_true, ep, free_settings.clone());
+
+ // ds_batch is the mean of the gradient of the objective function for each sample
+ let mut s: Vect<T, DS> = self.model.get_settings().clone().into();
+ s -= ds_batch * rate;
+ *self.model.get_settings_mut() = s.into();
+ }
+ }
+
+ /// Calibrate settings using stochastic GD
+ pub fn calibrate_stochastic(
+ &mut self,
+ ylist_true: &[Vect<T, D>],
+ ep: T,
+ rate: T,
+ niters: usize,
+ free_settings: impl Iterator<Item = usize> + Clone,
+ ) {
+ for _ in 0..niters {
+ for (x, y_true) in ylist_true.iter().zip(ylist_true.iter().skip(1)) {
+ let ds = self.objective_gradient(*x, *y_true, ep, free_settings.clone());
+
+ let mut s: Vect<T, DS> = self.model.get_settings().clone().into();
+ s -= ds * rate;
+ *self.model.get_settings_mut() = s.into();
+ }
+ }
+ }
+
+ /// Calibrate settings using batch GD and record path and error
+ ///
+ /// Returns (path, error)
+ pub fn calibrate_batch_record(
+ &mut self,
+ ylist_true: &[Vect<T, D>],
+ ep: T,
+ rate: T,
+ niters: usize,
+ free_settings: impl Iterator<Item = usize> + Clone,
+ ) -> (Vec<Vect<T, DS>>, Vec<T>)
+ where
+ T: PartialOrd,
+ {
+ let mut path = Vec::with_capacity(niters + 1);
+ path.push(self.model.get_settings().clone().into());
+ let mut errorlist = Vec::with_capacity(niters + 1);
+ errorlist.push(self.objective_batch(&self.model, ylist_true));
+
+ for _ in 0..niters {
+ let ds_batch = self.objective_gradient_batch(ylist_true, ep, free_settings.clone());
+
+ // ds_batch is the mean of the gradient of the objective function for each sample
+ let mut s: Vect<T, DS> = self.model.get_settings().clone().into();
+ s -= ds_batch * rate;
+ path.push(s);
+ *self.model.get_settings_mut() = s.into();
+
+ let error = self.objective_batch(&self.model, ylist_true);
+ if error > *errorlist.last().unwrap() {
+ path.pop();
+ *self.model.get_settings_mut() = (*path.last().unwrap()).into();
+ break;
+ }
+ errorlist.push(error);
+ }
+ (path, errorlist)
+ }
+
+ /// Calibrate settings using stochastic GD and record path and error
+ ///
+ /// Returns (path, error)
+ pub fn calibrate_stochastic_record(
+ &mut self,
+ ylist_true: &[Vect<T, D>],
+ ep: T,
+ rate: T,
+ niters: usize,
+ free_settings: impl Iterator<Item = usize> + Clone,
+ ) -> (Vec<Vect<T, DS>>, Vec<T>)
+ where
+ T: PartialOrd,
+ {
+ let mut path = Vec::with_capacity(niters * (niters - 1) + 1);
+ path.push(self.model.get_settings().clone().into());
+ let mut errorlist = Vec::with_capacity(niters * (niters - 1) + 1);
+ errorlist.push(self.objective_batch(&self.model, ylist_true));
+
+ for _ in 0..niters {
+ for (x, y_true) in ylist_true.iter().zip(ylist_true.iter().skip(1)) {
+ let ds = self.objective_gradient(*x, *y_true, ep, free_settings.clone());
+
+ let mut s: Vect<T, DS> = self.model.get_settings().clone().into();
+ s -= ds * rate;
+ path.push(s);
+ *self.model.get_settings_mut() = s.into();
+
+ let error = self.objective_batch(&self.model, ylist_true);
+ if error > *errorlist.last().unwrap() {
+ path.pop();
+ *self.model.get_settings_mut() = (*path.last().unwrap()).into();
+ continue;
+ }
+ errorlist.push(error);
+ }
+ }
+ (path, errorlist)
+ }
+
+ /// Mean of the objective function on all the samples
+ pub fn objective_batch(&self, model: &M, ylist_true: &[Vect<T, D>]) -> T {
+ let nsamples = T::from_usize(ylist_true.len() - 1).unwrap();
+
+ let mut obj_batch = T::zero();
+ for (x, y_true) in ylist_true.iter().zip(ylist_true.iter().skip(1)) {
+ obj_batch += self.objective(model, *x, *y_true);
+ }
+ obj_batch / nsamples
+ }
+}
+
+#[cfg(test)]
+mod test {
+ use super::*;
+ use nalgebra::vector;
+ use rand::Rng;
+
+ /// This test can fail, but should succeed most of the time
+ #[test]
+ fn test_objective_gradient_convergence() {
+ let mut rng = rand::thread_rng();
+
+ let ep0 = 0.01;
+ let nep = 20;
+ let ntests = 1000;
+
+ let mut fail_spikes = 0;
+ let mut fail_d = 0;
+
+ for _ in 0..ntests {
+ let model = sir::Sir {
+ s: sir::SirSettings {
+ beta: rng.gen(),
+ gamma: rng.gen(),
+ pop: 1.0,
+ },
+ };
+ let solver = ImplicitEulerSolver {
+ dt: 0.1,
+ tol: 0.000001,
+ niters: 100,
+ };
+ let optimizer = GradientDescentOptimizer::new(model, solver);
+ let x = rng.gen();
+ let y_true = rng.gen();
+
+ let mut g = optimizer.objective_gradient(x, y_true, ep0, 0..2);
+ let mut d = f64::MAX;
+ let mut spikes = 5;
+ for ep in (1..nep).map(|i| ep0 / 2.0.powi(i)) {
+ let ng = optimizer.objective_gradient(x, y_true, dbg!(ep), 0..2);
+ let nd = (dbg!(ng) - g).norm();
+ if nd.is_zero() {
+ break;
+ }
+ // Allow obj' having a local minimum between s and s+ep
+ if dbg!(nd) >= dbg!(d) {
+ if spikes == 0 {
+ fail_spikes += 1;
+ break;
+ }
+ spikes -= 1;
+ }
+ g = ng;
+ d = nd;
+ }
+ // d should be very small
+ if d > 10.0 * ep0 / 2.0.powi(nep - 1) {
+ fail_d += 1;
+ }
+ }
+
+ let prop_fail_spikes = fail_spikes as f64 / ntests as f64;
+ let prop_fail_d = fail_d as f64 / ntests as f64;
+ println!("Fail spikes: {} %", prop_fail_spikes * 100.0);
+ println!("Fail d: {} %", prop_fail_d * 100.0);
+ assert!(prop_fail_spikes < 0.015);
+ assert!(prop_fail_d < 0.0015);
+ }
+
+ // TODO fix
+ //#[test]
+ fn _test_objective_gradient_direction() {
+ let mut rng = rand::thread_rng();
+
+ let niters: usize = 1000;
+ let x0 = vector![0.99, 0.01];
+ let dt = 0.1;
+
+ // Generate "true" data
+ let settings_true = sir::SirSettings {
+ beta: rng.gen(),
+ gamma: rng.gen(),
+ pop: 1.0,
+ };
+ let model = sir::Sir {
+ s: dbg!(settings_true),
+ };
+ let solver = ImplicitEulerSolver {
+ dt,
+ tol: 0.000001,
+ niters: 100,
+ };
+ let mut optimizer = GradientDescentOptimizer::new(model, solver);
+ let mut xlist_true = Vec::with_capacity(niters + 1);
+ xlist_true.push(x0);
+ let mut x = x0;
+ for _ in 0..niters {
+ x = optimizer.solver.f(&optimizer.model, x);
+ xlist_true.push(x);
+ }
+
+ // Start with random settings
+ *optimizer.model.get_settings_mut() = dbg!(sir::SirSettings {
+ beta: rng.gen(),
+ gamma: rng.gen(),
+ pop: 1.0,
+ });
+
+ // Compute descent direction
+ let dir = dbg!(optimizer.objective_gradient(x0, xlist_true[1], 0.0000001, 0..2));
+
+ // Check that this direction leads to smaller error
+ let s: Vect<f64, 3> = optimizer.model.get_settings().clone().into();
+ let y = optimizer.model.f(x0);
+ *optimizer.model.get_settings_mut() = (s - dir).into(); // Apply direction
+ let y_new = optimizer.model.f(x0);
+ assert!((y - xlist_true[1]).norm() >= (y_new - xlist_true[1]).norm());
+ }
+}
diff --git a/examples/old/solver.rs b/examples/old/solver.rs
new file mode 100644
index 0000000..c16d97e
--- /dev/null
+++ b/examples/old/solver.rs
@@ -0,0 +1,41 @@
+use crate::{model::*, opti::*, utils::*};
+
+use nalgebra::{base::*, ComplexField};
+
+pub trait Solver<T, S: Settings, M: Model<T, S, D>, const D: usize> {
+ fn f(&self, model: &M, x: Vect<T, D>) -> Vect<T, D>;
+}
+
+#[derive(Clone)]
+pub struct ExplicitEulerSolver<T: Copy + ComplexField<RealField = T>> {
+ pub dt: T,
+}
+
+impl<T: Copy + ComplexField<RealField = T>, S: Settings, M: Model<T, S, D>, const D: usize>
+ Solver<T, S, M, D> for ExplicitEulerSolver<T>
+{
+ fn f(&self, model: &M, x: Vect<T, D>) -> Vect<T, D> {
+ model.f(x) * self.dt + x
+ }
+}
+
+#[derive(Clone)]
+pub struct ImplicitEulerSolver<T: Clone + ComplexField<RealField = T>> {
+ pub dt: T,
+ pub tol: T,
+ pub niters: usize,
+}
+
+impl<
+ T: Copy + ComplexField<RealField = T> + PartialOrd,
+ S: Settings,
+ M: Model<T, S, D>,
+ const D: usize,
+ > Solver<T, S, M, D> for ImplicitEulerSolver<T>
+where
+ Const<D>: ToTypenum + DimMin<Const<D>, Output = Const<D>>,
+{
+ fn f(&self, model: &M, x: Vect<T, D>) -> Vect<T, D> {
+ newton(model, x, self.dt, self.tol, self.niters)
+ }
+}
diff --git a/examples/old/space.rs b/examples/old/space.rs
new file mode 100644
index 0000000..1bd6a5f
--- /dev/null
+++ b/examples/old/space.rs
@@ -0,0 +1,94 @@
+use crate::{
+ model::{Model, Settings},
+ solver::Solver,
+ utils::*,
+};
+
+use nalgebra::ComplexField;
+use std::{collections::HashMap, marker::PhantomData};
+
+#[derive(Clone)]
+pub struct Point<T: Copy, const D: usize> {
+ pub pos: Vect<T, D>,
+ pub diffusion: Vect<T, D>,
+}
+
+#[derive(Clone)]
+pub struct Source<T: Copy, const D: usize> {
+ pub amplitude: T,
+ pub fonction: fn(T) -> Vect<T, D>,
+ pub freq: T,
+ pub phase: T,
+}
+
+pub struct Space<T: Copy, S: Settings, M: Model<T, S, D>, V: Solver<T, S, M, D>, const D: usize> {
+ pub model: M,
+ pub old_points: Vec<Point<T, D>>,
+ pub points: Vec<Point<T, D>>,
+ pub size: (usize, usize),
+ pub solver: V,
+ pub sources: HashMap<usize, Source<T, D>>,
+ pub time: T,
+ pub _p: PhantomData<S>,
+}
+
+impl<
+ T: Copy + ComplexField<RealField = T>,
+ S: Settings,
+ M: Model<T, S, D>,
+ V: Solver<T, S, M, D>,
+ const D: usize,
+ > Space<T, S, M, V, D>
+{
+ pub fn simulate(&mut self, delta_time: T) {
+ std::mem::swap(&mut self.old_points, &mut self.points);
+ for (i, (point, old_point)) in self
+ .points
+ .iter_mut()
+ .zip(self.old_points.iter())
+ .enumerate()
+ {
+ if let Some(source) = self.sources.get(&i) {
+ let t = self.time * source.freq + source.phase;
+ //Point{pos:t.sin()*source.amplitude, speed: t.cos()*source.amplitude}
+ point.pos = (source.fonction)(t) * source.amplitude;
+ } else {
+ *point = old_point.clone()
+ };
+ point.pos = self.solver.f(&self.model, point.pos);
+ }
+
+ let mut i = 0usize;
+ for y in 0..self.size.1 {
+ for x in 0..self.size.0 {
+ let old_point = self.old_points[i].clone();
+ let point = &mut self.points[i];
+
+ if y > 0 {
+ point.pos += (self.old_points[i - self.size.0].pos - old_point.pos)
+ .component_mul(&point.diffusion)
+ * delta_time;
+ }
+ if y < self.size.1 - 1 {
+ point.pos += (self.old_points[i + self.size.0].pos - old_point.pos)
+ .component_mul(&point.diffusion)
+ * delta_time;
+ }
+ if x > 0 {
+ point.pos += (self.old_points[i - 1].pos - old_point.pos)
+ .component_mul(&point.diffusion)
+ * delta_time;
+ }
+ if x < self.size.0 - 1 {
+ point.pos += (self.old_points[i + 1].pos - old_point.pos)
+ .component_mul(&point.diffusion)
+ * delta_time;
+ }
+
+ i += 1;
+ }
+ }
+
+ self.time += delta_time;
+ }
+}
diff --git a/examples/old/utils.rs b/examples/old/utils.rs
new file mode 100644
index 0000000..8e051eb
--- /dev/null
+++ b/examples/old/utils.rs
@@ -0,0 +1,16 @@
+use nalgebra::base::*;
+
+/// `<Type, Rows, Columns>`
+pub type Mat<T, const R: usize, const C: usize> =
+ Matrix<T, Const<R>, Const<C>, ArrayStorage<T, R, C>>;
+pub type Vect<T, const R: usize> = Matrix<T, Const<R>, U1, ArrayStorage<T, R, 1>>;
+
+pub fn max(l: &[f64]) -> f64 {
+ let mut m = l[0];
+ for &i in &l[1..] {
+ if i > m {
+ m = i;
+ }
+ }
+ m
+}
diff --git a/rustfmt.toml b/rustfmt.toml
new file mode 100644
index 0000000..843c6fe
--- /dev/null
+++ b/rustfmt.toml
@@ -0,0 +1,9 @@
+hard_tabs = true
+newline_style = "unix"
+imports_granularity = "Crate"
+
+unstable_features = true
+format_code_in_doc_comments = true
+format_macro_bodies = true
+format_macro_matchers = true
+format_strings = true
diff --git a/src/lib.rs b/src/lib.rs
new file mode 100644
index 0000000..76bc5e1
--- /dev/null
+++ b/src/lib.rs
@@ -0,0 +1,20 @@
+pub mod model;
+#[cfg(feature = "plot")]
+pub mod plot;
+pub mod solver;
+mod util;
+
+pub mod prelude {
+ #[cfg(feature = "plot")]
+ pub use crate::plot::{Plot, PlotBuilder, PlotBuilderError};
+ pub use crate::{
+ model::{Model, Settings},
+ solver::{
+ euler::{ExplicitEuler, ImplicitEuler},
+ Solver,
+ },
+ util::{Mat, Vect},
+ };
+ #[cfg(feature = "plot")]
+ pub use plotters::style::colors;
+}
diff --git a/src/model.rs b/src/model.rs
new file mode 100644
index 0000000..b2c4d6f
--- /dev/null
+++ b/src/model.rs
@@ -0,0 +1,12 @@
+use crate::util::*;
+
+pub trait Settings {}
+
+pub trait Model<T, S: Settings, const D: usize> {
+ /// Returns f(x)
+ fn f(&self, x: Vect<T, D>) -> Vect<T, D>;
+ /// Returns df(x)/dx
+ fn df(&self, x: Vect<T, D>) -> Mat<T, D, D>;
+ fn get_settings(&self) -> &S;
+ fn get_settings_mut(&mut self) -> &mut S;
+}
diff --git a/src/plot.rs b/src/plot.rs
new file mode 100644
index 0000000..a8fde8a
--- /dev/null
+++ b/src/plot.rs
@@ -0,0 +1,107 @@
+use std::cmp::Ordering;
+
+use crate::util::*;
+
+use nalgebra::Matrix;
+use plotters::prelude::*;
+
+#[derive(Clone, Debug, derive_builder::Builder)]
+#[builder(setter(into), default)]
+pub struct Plot {
+ dt: f64,
+ size: (u32, u32),
+ title: Option<String>,
+ x_label: Option<String>,
+ x_max: Option<f64>,
+ x_min: Option<f64>,
+ y_max: Option<f64>,
+ y_min: Option<f64>,
+}
+
+impl Default for Plot {
+ fn default() -> Self {
+ Self {
+ dt: 1.0,
+ size: (640, 480),
+ title: None,
+ x_label: None,
+ x_max: None,
+ x_min: None,
+ y_max: None,
+ y_min: None,
+ }
+ }
+}
+
+impl Plot {
+ pub fn plot<const R: usize>(
+ &self,
+ path: impl AsRef<std::path::Path>,
+ data: &[Vect<f64, R>],
+ series: &[(&str, RGBColor); R],
+ ) {
+ let root = BitMapBackend::new(path.as_ref(), self.size).into_drawing_area();
+ root.fill(&WHITE).unwrap();
+ let mut chart = ChartBuilder::on(&root);
+ if let Some(title) = &self.title {
+ chart.caption(
+ title,
+ FontDesc::new(FontFamily::SansSerif, 22.0, FontStyle::Normal),
+ );
+ }
+
+ let x_range =
+ self.x_min.unwrap_or(0.0)..self.x_max.unwrap_or_else(|| data.len() as f64 * self.dt);
+ let y_range = self
+ .y_min
+ .unwrap_or_else(|| find_extremum(data.iter().map(Matrix::min), Ordering::Less).unwrap())
+ ..self.y_max.unwrap_or_else(|| {
+ find_extremum(data.iter().map(Matrix::max), Ordering::Greater).unwrap()
+ });
+
+ let mut chart = chart
+ .margin_right(12)
+ .y_label_area_size(30)
+ .x_label_area_size(30)
+ .build_cartesian_2d(x_range, y_range)
+ .unwrap();
+
+ let mut x_axis = chart.configure_mesh();
+ if let Some(x_label) = &self.x_label {
+ x_axis.x_desc(x_label);
+ }
+ x_axis.draw().unwrap();
+
+ for (i, (label, color)) in series.into_iter().enumerate() {
+ chart
+ .draw_series(LineSeries::new(
+ data.iter()
+ .enumerate()
+ .map(|(x, y)| (x as f64 * self.dt, y[i])),
+ color,
+ ))
+ .unwrap()
+ .label(*label)
+ .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], color.clone()));
+ }
+
+ chart
+ .configure_series_labels()
+ .background_style(WHITE.mix(0.8))
+ .border_style(BLACK)
+ .draw()
+ .unwrap();
+ }
+}
+
+fn find_extremum<T: PartialOrd, I: IntoIterator<Item = T>>(iter: I, ord: Ordering) -> Option<T> {
+ let mut iter = iter.into_iter();
+ let extremum = iter.next()?;
+ Some(iter.fold(extremum, |extremum, i| {
+ if i.partial_cmp(&extremum) == Some(ord) {
+ i
+ } else {
+ extremum
+ }
+ }))
+}
diff --git a/src/solver.rs b/src/solver.rs
new file mode 100644
index 0000000..05b26db
--- /dev/null
+++ b/src/solver.rs
@@ -0,0 +1,7 @@
+pub mod euler;
+
+use crate::{model::*, util::*};
+
+pub trait Solver<T, S: Settings, M: Model<T, S, D>, const D: usize> {
+ fn f(&self, model: &M, x: Vect<T, D>) -> Vect<T, D>;
+}
diff --git a/src/solver/euler.rs b/src/solver/euler.rs
new file mode 100644
index 0000000..65307c8
--- /dev/null
+++ b/src/solver/euler.rs
@@ -0,0 +1,41 @@
+use crate::{
+ model::{Model, Settings},
+ solver::Solver,
+ util::*,
+};
+
+use nalgebra::{base::*, ComplexField};
+
+#[derive(Clone)]
+pub struct ExplicitEuler<T: Copy + ComplexField<RealField = T>> {
+ pub dt: T,
+}
+
+impl<T: Copy + ComplexField<RealField = T>, S: Settings, M: Model<T, S, D>, const D: usize>
+ Solver<T, S, M, D> for ExplicitEuler<T>
+{
+ fn f(&self, model: &M, x: Vect<T, D>) -> Vect<T, D> {
+ model.f(x) * self.dt + x
+ }
+}
+
+#[derive(Clone)]
+pub struct ImplicitEuler<T: Clone + ComplexField<RealField = T>> {
+ pub dt: T,
+ pub tol: T,
+ pub niters: usize,
+}
+
+impl<
+ T: Copy + ComplexField<RealField = T> + PartialOrd,
+ S: Settings,
+ M: Model<T, S, D>,
+ const D: usize,
+ > Solver<T, S, M, D> for ImplicitEuler<T>
+where
+ Const<D>: ToTypenum + DimMin<Const<D>, Output = Const<D>>,
+{
+ fn f(&self, model: &M, x: Vect<T, D>) -> Vect<T, D> {
+ newton(model, x, self.dt, self.tol, self.niters)
+ }
+}
diff --git a/src/util.rs b/src/util.rs
new file mode 100644
index 0000000..75184bd
--- /dev/null
+++ b/src/util.rs
@@ -0,0 +1,39 @@
+use crate::model::{Model, Settings};
+
+use nalgebra::{base::*, ComplexField};
+
+/// Matrix `<Type, Rows, Columns>`
+pub type Mat<T, const R: usize, const C: usize> =
+ Matrix<T, Const<R>, Const<C>, ArrayStorage<T, R, C>>;
+/// Vector `<Type, Dimension>` (column matrix)
+pub type Vect<T, const R: usize> = Matrix<T, Const<R>, U1, ArrayStorage<T, R, 1>>;
+
+pub fn newton<
+ T: Copy + Scalar + ComplexField<RealField = T> + PartialOrd,
+ S: Settings,
+ const D: usize,
+>(
+ model: &impl Model<T, S, D>,
+ x0: Vect<T, D>,
+ dt: T,
+ tol: T,
+ niters: usize,
+) -> Vect<T, D>
+where
+ Const<D>: ToTypenum + DimMin<Const<D>, Output = Const<D>>,
+{
+ let mut x = x0;
+
+ for _ in 0..niters {
+ if let Some(m) = (Mat::<T, D, D>::identity() - model.df(x) * dt).try_inverse() {
+ let dx = m * (x - x0 - model.f(x) * dt);
+ if dx.norm() < tol {
+ break;
+ }
+ x -= dx;
+ } else {
+ break;
+ }
+ }
+ x
+}