changeset
merge dev to default
Mon, 03 Feb 2025 19:22:16 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Mon, 03 Feb 2025 19:22:16 -0500
- changeset 90
- b3c35d16affe
- parent 25
- d14c877e14b7 (current diff)
- parent 89
- 8513a10c5dd9 (diff)
- child 91
- db870f2a2cde
merge dev to default
--- a/Cargo.lock Tue Feb 20 12:33:16 2024 -0500 +++ b/Cargo.lock Mon Feb 03 19:22:16 2025 -0500 @@ -1,11 +1,12 @@ # This file is automatically @generated by Cargo. # It is not intended for manual editing. -version = 3 +version = 4 [[package]] name = "alg_tools" -version = "0.1.0" +version = "0.3.0" dependencies = [ + "anyhow", "colored", "cpu-time", "csv", @@ -17,10 +18,16 @@ "rayon", "serde", "serde_json", - "trait-set", + "simba", ] [[package]] +name = "anyhow" +version = "1.0.95" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34ac096ce696dc2fcabef30516bb13c0a68a11d30131d3df6f04711467681b04" + +[[package]] name = "approx" version = "0.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -31,47 +38,22 @@ [[package]] name = "autocfg" -version = "1.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" - -[[package]] -name = "bitflags" -version = "2.4.0" +version = "1.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b4682ae6287fcf752ecaabbfcc7b6f9b72aa33933dc23a554d853aea8eea8635" - -[[package]] -name = "bstr" -version = "0.2.17" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ba3569f383e8f1598449f1a423e72e99569137b47740b1da11ef19af3d5c3223" -dependencies = [ - "lazy_static", - "memchr", - "regex-automata", - "serde", -] +checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" [[package]] name = "bytemuck" -version = "1.14.0" +version = "1.20.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "374d28ec25809ee0e23827c2ab573d729e293f281dfe393500e7ad618baa61c6" - -[[package]] -name = "cfg-if" -version = "1.0.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" +checksum = "8b37c88a63ffd85d15b406896cc343916d7cf57838a847b3a6f2ca5d39a5695a" [[package]] name = "colored" -version = "2.0.4" +version = "2.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2674ec482fbc38012cf31e6c42ba0177b431a0cb6f15fe40efa5aab1bda516f6" +checksum = "cbf2150cce219b664a8a70df7a1f933836724b503f8a413af9365b4dcc4d90b8" dependencies = [ - "is-terminal", "lazy_static", "windows-sys", ] @@ -88,46 +70,37 @@ [[package]] name = "crossbeam-deque" -version = "0.8.3" +version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ce6fd6f855243022dcecf8702fef0c297d4338e226845fe067f6341ad9fa0cef" +checksum = "613f8cc01fe9cf1a3eb3d7f488fd2fa8388403e97039e2f73692932e291a770d" dependencies = [ - "cfg-if", "crossbeam-epoch", "crossbeam-utils", ] [[package]] name = "crossbeam-epoch" -version = "0.9.15" +version = "0.9.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ae211234986c545741a7dc064309f67ee1e5ad243d0e48335adc0484d960bcc7" +checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" dependencies = [ - "autocfg", - "cfg-if", "crossbeam-utils", - "memoffset", - "scopeguard", ] [[package]] name = "crossbeam-utils" -version = "0.8.16" +version = "0.8.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5a22b2d63d4d1dc0b7f1b6b2747dd0088008a9be28b6ddf0b1e7d335e3037294" -dependencies = [ - "cfg-if", -] +checksum = "22ec99545bb0ed0ea7bb9b8e1e9122ea386ff8a48c0922e43f36d45ab09e0e80" [[package]] name = "csv" -version = "1.1.6" +version = "1.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "22813a6dc45b335f9bade10bf7271dc477e81113e89eb251a0bc2a8a81c536e1" +checksum = "acdc4883a9c96732e4733212c01447ebd805833b7275a73ca3ee080fd77afdaf" dependencies = [ - "bstr", "csv-core", - "itoa 0.4.8", + "itoa", "ryu", "serde", ] @@ -143,81 +116,42 @@ [[package]] name = "either" -version = "1.9.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a26ae43d7bcc3b814de94796a5e736d4029efb0ee900c12e2d54c993ad1a1e07" - -[[package]] -name = "errno" -version = "0.3.5" +version = "1.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac3e13f66a2f95e32a39eaa81f6b95d42878ca0e1db0c7543723dfe12557e860" -dependencies = [ - "libc", - "windows-sys", -] - -[[package]] -name = "hermit-abi" -version = "0.3.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d77f7ec81a6d05a3abb01ab6eb7590f6083d08449fe5a1c8b1e620283546ccb7" - -[[package]] -name = "is-terminal" -version = "0.4.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cb0889898416213fab133e1d33a0e5858a48177452750691bde3666d0fdbaf8b" -dependencies = [ - "hermit-abi", - "rustix", - "windows-sys", -] +checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" [[package]] name = "itertools" -version = "0.10.5" +version = "0.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b0fd2260e829bddf4cb6ea802289de2f86d6a7a690192fbe91b3f46e0f2c8473" +checksum = "413ee7dfc52ee1a4949ceeb7dbc8a33f2d6c088194d9f922fb8318faf1f01186" dependencies = [ "either", ] [[package]] name = "itoa" -version = "0.4.8" +version = "1.0.14" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b71991ff56294aa922b450139ee08b3bfc70982c6b2c7562771375cf73542dd4" - -[[package]] -name = "itoa" -version = "1.0.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "af150ab688ff2122fcef229be89cb50dd66af9e01a4ff320cc137eecc9bacc38" +checksum = "d75a2a4b1b190afb6f5425f10f6a8f959d2ea0b9c2b1d79553551850539e4674" [[package]] name = "lazy_static" -version = "1.4.0" +version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" [[package]] name = "libc" -version = "0.2.149" +version = "0.2.168" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a08173bc88b7955d1b3145aa561539096c421ac8debde8cbc3612ec635fee29b" - -[[package]] -name = "linux-raw-sys" -version = "0.4.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "da2479e8c062e40bf0066ffa0bc823de0a9368974af99c9f6df941d2c231e03f" +checksum = "5aaeb2981e0606ca11d79718f8bb01164f1d6ed75080182d3abf017e6d244b6d" [[package]] name = "matrixmultiply" -version = "0.3.8" +version = "0.3.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2" +checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a" dependencies = [ "autocfg", "rawpointer", @@ -225,24 +159,15 @@ [[package]] name = "memchr" -version = "2.6.4" +version = "2.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f665ee40bc4a3c5590afb1e9677db74a508659dfd71e126420da8274909a0167" - -[[package]] -name = "memoffset" -version = "0.9.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c" -dependencies = [ - "autocfg", -] +checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" [[package]] name = "nalgebra" -version = "0.31.4" +version = "0.33.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "20bd243ab3dbb395b39ee730402d2e5405e448c75133ec49cc977762c4cba3d1" +checksum = "26aecdf64b707efd1310e3544d709c5c0ac61c13756046aaaba41be5c4f66a3b" dependencies = [ "approx", "matrixmultiply", @@ -256,20 +181,20 @@ [[package]] name = "nalgebra-macros" -version = "0.1.0" +version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218" +checksum = "254a5372af8fc138e36684761d3c0cdb758a4410e938babcff1c860ce14ddbfc" dependencies = [ "proc-macro2", "quote", - "syn 1.0.109", + "syn 2.0.90", ] [[package]] name = "num" -version = "0.4.1" +version = "0.4.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b05180d69e3da0e530ba2a1dae5110317e49e3b7f3d41be227dc5f92e49ee7af" +checksum = "35bd024e8b2ff75562e5f34e7f4905839deb4b22955ef5e73d2fea1b9813cb23" dependencies = [ "num-bigint", "num-complex", @@ -281,39 +206,37 @@ [[package]] name = "num-bigint" -version = "0.4.4" +version = "0.4.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "608e7659b5c3d7cba262d894801b9ec9d00de989e8a82bd4bef91d08da45cdc0" +checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9" dependencies = [ - "autocfg", "num-integer", "num-traits", ] [[package]] name = "num-complex" -version = "0.4.4" +version = "0.4.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1ba157ca0885411de85d6ca030ba7e2a83a28636056c7c699b07c8b6f7383214" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" dependencies = [ "num-traits", ] [[package]] name = "num-integer" -version = "0.1.45" +version = "0.1.46" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" dependencies = [ - "autocfg", "num-traits", ] [[package]] name = "num-iter" -version = "0.1.43" +version = "0.1.45" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7d03e6c028c5dc5cac6e2dec0efda81fc887605bb3d884578bb6d6bf7514e252" +checksum = "1429034a0490724d0075ebb2bc9e875d6503c3cf69e235a8941aa757d83ef5bf" dependencies = [ "autocfg", "num-integer", @@ -322,11 +245,10 @@ [[package]] name = "num-rational" -version = "0.4.1" +version = "0.4.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" dependencies = [ - "autocfg", "num-bigint", "num-integer", "num-traits", @@ -334,9 +256,9 @@ [[package]] name = "num-traits" -version = "0.2.17" +version = "0.2.19" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "39e3200413f237f41ab11ad6d161bc7239c84dcb631773ccd7de3dfe4b5c267c" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", ] @@ -353,24 +275,24 @@ [[package]] name = "paste" -version = "1.0.14" +version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" [[package]] name = "proc-macro2" -version = "1.0.69" +version = "1.0.92" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "134c189feb4956b20f6f547d2cf727d4c0fe06722b20a0eec87ed445a97f92da" +checksum = "37d3544b3f2748c54e147655edb5025752e2303145b5aefb3c3ea2c78b973bb0" dependencies = [ "unicode-ident", ] [[package]] name = "quote" -version = "1.0.33" +version = "1.0.37" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae" +checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af" dependencies = [ "proc-macro2", ] @@ -383,9 +305,9 @@ [[package]] name = "rayon" -version = "1.8.0" +version = "1.10.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c27db03db7734835b3f53954b534c91069375ce6ccaa2e065441e07d9b6cdb1" +checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa" dependencies = [ "either", "rayon-core", @@ -393,90 +315,66 @@ [[package]] name = "rayon-core" -version = "1.12.0" +version = "1.12.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5ce3fb6ad83f861aac485e76e1985cd109d9a3713802152be56c3b1f0e0658ed" +checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2" dependencies = [ "crossbeam-deque", "crossbeam-utils", ] [[package]] -name = "regex-automata" -version = "0.1.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6c230d73fb8d8c1b9c0b3135c5142a8acee3a0558fb8db5cf1cb65f8d7862132" - -[[package]] -name = "rustix" -version = "0.38.19" +name = "ryu" +version = "1.0.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "745ecfa778e66b2b63c88a61cb36e0eea109e803b0b86bf9879fbc77c70e86ed" -dependencies = [ - "bitflags", - "errno", - "libc", - "linux-raw-sys", - "windows-sys", -] - -[[package]] -name = "ryu" -version = "1.0.15" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1ad4cc8da4ef723ed60bced201181d83791ad433213d8c24efffda1eec85d741" +checksum = "f3cb5ba0dc43242ce17de99c180e96db90b235b8a9fdc9543c96d2209116bd9f" [[package]] name = "safe_arch" -version = "0.7.1" +version = "0.7.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354" +checksum = "c3460605018fdc9612bce72735cba0d27efbcd9904780d44c7e3a9948f96148a" dependencies = [ "bytemuck", ] [[package]] -name = "scopeguard" -version = "1.2.0" +name = "serde" +version = "1.0.216" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" - -[[package]] -name = "serde" -version = "1.0.189" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e422a44e74ad4001bdc8eede9a4570ab52f71190e9c076d14369f38b9200537" +checksum = "0b9781016e935a97e8beecf0c933758c97a5520d32930e460142b4cd80c6338e" dependencies = [ "serde_derive", ] [[package]] name = "serde_derive" -version = "1.0.189" +version = "1.0.216" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1e48d1f918009ce3145511378cf68d613e3b3d9137d67272562080d68a2b32d5" +checksum = "46f859dbbf73865c6627ed570e78961cd3ac92407a2d117204c49232485da55e" dependencies = [ "proc-macro2", "quote", - "syn 2.0.38", + "syn 2.0.90", ] [[package]] name = "serde_json" -version = "1.0.107" +version = "1.0.133" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6b420ce6e3d8bd882e9b243c6eed35dbc9a6110c9769e74b584e0d68d1f20c65" +checksum = "c7fceb2473b9166b2294ef05efcb65a3db80803f0b03ef86a5fc88a2b85ee377" dependencies = [ - "itoa 1.0.9", + "itoa", + "memchr", "ryu", "serde", ] [[package]] name = "simba" -version = "0.7.3" +version = "0.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2f3fd720c48c53cace224ae62bef1bbff363a70c68c4802a78b5cc6159618176" +checksum = "b3a386a501cd104797982c15ae17aafe8b9261315b5d07e3ec803f2ea26be0fa" dependencies = [ "approx", "num-complex", @@ -498,9 +396,9 @@ [[package]] name = "syn" -version = "2.0.38" +version = "2.0.90" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e96b79aaa137db8f61e26363a0c9b47d8b4ec75da28b7d1d614c2303e232408b" +checksum = "919d3b74a5dd0ccd15aeb8f93e7006bd9e14c295087c9896a110f490752bcf31" dependencies = [ "proc-macro2", "quote", @@ -508,17 +406,6 @@ ] [[package]] -name = "trait-set" -version = "0.2.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "875c4c873cc824e362fa9a9419ffa59807244824275a44ad06fec9684fff08f2" -dependencies = [ - "proc-macro2", - "quote", - "syn 1.0.109", -] - -[[package]] name = "typenum" version = "1.17.0" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -526,15 +413,15 @@ [[package]] name = "unicode-ident" -version = "1.0.12" +version = "1.0.14" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" +checksum = "adb9e6ca4f869e1180728b7950e35922a7fc6397f7b641499e8f3ef06e50dc83" [[package]] name = "wide" -version = "0.7.12" +version = "0.7.30" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ebecebefc38ff1860b4bc47550bbfa63af5746061cf0d29fcd7fa63171602598" +checksum = "58e6db2670d2be78525979e9a5f9c69d296fd7d670549fe9ebf70f8708cb5019" dependencies = [ "bytemuck", "safe_arch",
--- a/Cargo.toml Tue Feb 20 12:33:16 2024 -0500 +++ b/Cargo.toml Mon Feb 03 19:22:16 2025 -0500 @@ -1,8 +1,8 @@ [package] name = "alg_tools" -version = "0.1.0" +version = "0.3.0" edition = "2021" -rust-version = "1.67" +rust-version = "1.85" authors = ["Tuomo Valkonen <tuomov@iki.fi>"] description = "Utility routines for iterative numerical algorithms" homepage = "https://tuomov.iki.fi/software/alg_tools/" @@ -13,17 +13,18 @@ [dependencies] serde = { version = "1.0", features = ["derive"] } -csv = "~1.1.6" -nalgebra = "~0.31.0" +csv = "~1.3.1" +nalgebra = "~0.33.0" num-traits = { version = "~0.2.14", features = ["std"] } -colored = "~2.0.0" -trait-set = "~0.2.0" +colored = "~2.1.0" num = "~0.4.0" -itertools = "~0.10.3" +itertools = "~0.13.0" numeric_literals = "~0.2.0" cpu-time = "~1.0.0" -serde_json = "~1.0.85" +serde_json = { version = "~1.0.85", features = ["std"] } rayon = "1.5.3" +simba = "0.9.0" +anyhow = "1.0.95" [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html" ] @@ -31,3 +32,9 @@ [profile.release] debug = true + +[features] +default = ["nightly"] +use_custom_thread_pool = [] +nightly = [] # enable for higher-performance implementations of some things +
--- a/README.md Tue Feb 20 12:33:16 2024 -0500 +++ b/README.md Mon Feb 03 19:22:16 2025 -0500 @@ -8,6 +8,10 @@ * [Linear operator], [mapping], [Euclidean space], and [norm] abstractions. Matrices and vectors are supported via [nalgebra]. There is also abstraction for [`AXPY`][AXPY] and [`GEMV`][GEMV] operations. + * A facility to create [`Instance`][Instance]s of other types, for easy application + of functions to various concrete and reference types. + * Abstraction of [Fenchel conjugates] and [proximal operators] of convex + functions. * Small (on stack) [vectors] and [cubes] that implement the relevant abstractions and vector space operations. * Multi-dimensional [linear grids], including the familiar-from-Matlab @@ -56,4 +60,7 @@ [associated traits]: crate::types [vectors]: crate::loc::Loc [cubes]: crate::sets::Cube + [Instance]: crate::instance::Instance + [Fenchel conjugates]: crate::convex::Conjugable + [proximal operators]: crate::convex::Prox
--- a/rust-toolchain.toml Tue Feb 20 12:33:16 2024 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,2 +0,0 @@ -[toolchain] -channel = "nightly"
--- a/src/bisection_tree/aggregator.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/bisection_tree/aggregator.rs Mon Feb 03 19:22:16 2025 -0500 @@ -4,6 +4,7 @@ use crate::types::*; use crate::sets::Set; +use crate::instance::Instance; /// Trait for aggregating information about a branch of a [bisection tree][super::BT]. /// @@ -136,10 +137,11 @@ } impl<F : Float> Set<F> for Bounds<F> { - fn contains(&self, item : &F) -> bool { + fn contains<I : Instance<F>>(&self, item : I) -> bool { + let v = item.own(); let &Bounds(l, u) = self; debug_assert!(l <= u); - l <= *item && *item <= u + l <= v && v <= u } }
--- a/src/bisection_tree/btfn.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/bisection_tree/btfn.rs Mon Feb 03 19:22:16 2025 -0500 @@ -4,7 +4,10 @@ use std::marker::PhantomData; use std::sync::Arc; use crate::types::Float; -use crate::mapping::{Apply, Mapping}; +use crate::mapping::{ + Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space, + BasicDecomposition, +}; //use crate::linops::{Apply, Linear}; use crate::sets::Set; use crate::sets::Cube; @@ -40,10 +43,22 @@ } impl<F : Float, G, BT, const N : usize> +Space for BTFN<F, G, BT, N> +where + G : SupportGenerator<F, N, Id=BT::Data>, + G::SupportType : LocalAnalysis<F, BT::Agg, N>, + BT : BTImpl<F, N> +{ + type Decomp = BasicDecomposition; +} + +impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N> -where G : SupportGenerator<F, N, Id=BT::Data>, - G::SupportType : LocalAnalysis<F, BT::Agg, N>, - BT : BTImpl<F, N> { +where + G : SupportGenerator<F, N, Id=BT::Data>, + G::SupportType : LocalAnalysis<F, BT::Agg, N>, + BT : BTImpl<F, N> +{ /// Create a new BTFN from a support generator and a pre-initialised bisection tree. /// @@ -386,42 +401,52 @@ make_btfn_unaryop!(Neg, neg); - - // -// Mapping +// Apply, Mapping, Differentiate // -impl<'a, F : Float, G, BT, V, const N : usize> Apply<&'a Loc<F, N>> +impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F, N>> for BTFN<F, G, BT, N> -where BT : BTImpl<F, N>, - G : SupportGenerator<F, N, Id=BT::Data>, - G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<&'a Loc<F, N>, Output = V>, - V : Sum { +where + BT : BTImpl<F, N>, + G : SupportGenerator<F, N, Id=BT::Data>, + G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>, + V : Sum + Space, +{ - type Output = V; + type Codomain = V; - fn apply(&self, x : &'a Loc<F, N>) -> Self::Output { - self.bt.iter_at(x) - .map(|&d| self.generator.support_for(d).apply(x)).sum() + fn apply<I : Instance<Loc<F,N>>>(&self, x : I) -> Self::Codomain { + let xc = x.cow(); + self.bt.iter_at(&*xc) + .map(|&d| self.generator.support_for(d).apply(&*xc)).sum() } } -impl<F : Float, G, BT, V, const N : usize> Apply<Loc<F, N>> +impl<F : Float, G, BT, V, const N : usize> DifferentiableImpl<Loc<F, N>> for BTFN<F, G, BT, N> -where BT : BTImpl<F, N>, - G : SupportGenerator<F, N, Id=BT::Data>, - G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<Loc<F, N>, Output = V>, - V : Sum { +where + BT : BTImpl<F, N>, + G : SupportGenerator<F, N, Id=BT::Data>, + G::SupportType : LocalAnalysis<F, BT::Agg, N> + + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>, + V : Sum + Space, +{ - type Output = V; + type Derivative = V; - fn apply(&self, x : Loc<F, N>) -> Self::Output { - self.bt.iter_at(&x) - .map(|&d| self.generator.support_for(d).apply(x)).sum() + fn differential_impl<I : Instance<Loc<F, N>>>(&self, x :I) -> Self::Derivative { + let xc = x.cow(); + self.bt.iter_at(&*xc) + .map(|&d| self.generator.support_for(d).differential(&*xc)) + .sum() } } +// +// GlobalAnalysis +// + impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg> for BTFN<F, G, BT, N> where BT : BTImpl<F, N>, @@ -480,7 +505,7 @@ /// /// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers. /// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`. -pub trait P2Minimise<U, F : Float> : Set<U> { +pub trait P2Minimise<U : Space, F : Float> : Set<U> { /// Minimise `g` over the set presented by `Self`. /// /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`. @@ -806,7 +831,7 @@ impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N> where BT : BTSearch<F, N, Agg=Bounds<F>>, G : SupportGenerator<F, N, Id=BT::Data>, - G::SupportType : Mapping<Loc<F, N>,Codomain=F> + G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>, Cube<F, N> : P2Minimise<Loc<F, N>, F> {
--- a/src/bisection_tree/either.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/bisection_tree/either.rs Mon Feb 03 19:22:16 2025 -0500 @@ -3,8 +3,14 @@ use std::sync::Arc; use crate::types::*; -use crate::mapping::Apply; -use crate::iter::{Mappable,MapF,MapZ}; +use crate::mapping::{ + Instance, + Mapping, + DifferentiableImpl, + DifferentiableMapping, + Space, +}; +use crate::iter::{Mappable, MapF, MapZ}; use crate::sets::Cube; use crate::loc::Loc; @@ -177,12 +183,17 @@ } } -impl<F, S1, S2, X> Apply<X> for EitherSupport<S1, S2> -where S1 : Apply<X, Output=F>, - S2 : Apply<X, Output=F> { - type Output = F; +impl<F, S1, S2, X> Mapping<X> for EitherSupport<S1, S2> +where + F : Space, + X : Space, + S1 : Mapping<X, Codomain=F>, + S2 : Mapping<X, Codomain=F>, +{ + type Codomain = F; + #[inline] - fn apply(&self, x : X) -> F { + fn apply<I : Instance<X>>(&self, x : I) -> F { match self { EitherSupport::Left(ref a) => a.apply(x), EitherSupport::Right(ref b) => b.apply(x), @@ -190,6 +201,24 @@ } } +impl<X, S1, S2, O> DifferentiableImpl<X> for EitherSupport<S1, S2> +where + O : Space, + X : Space, + S1 : DifferentiableMapping<X, DerivativeDomain=O>, + S2 : DifferentiableMapping<X, DerivativeDomain=O>, +{ + type Derivative = O; + + #[inline] + fn differential_impl<I : Instance<X>>(&self, x : I) -> O { + match self { + EitherSupport::Left(ref a) => a.differential(x), + EitherSupport::Right(ref b) => b.differential(x), + } + } +} + macro_rules! make_either_scalarop_rhs { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl<F : Float, G1, G2>
--- a/src/bisection_tree/refine.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/bisection_tree/refine.rs Mon Feb 03 19:22:16 2025 -0500 @@ -364,8 +364,13 @@ let mut container = container_arc.lock().unwrap(); // Safe: we just created arg_b and have a mutable exclusive // reference to self containing it. + #[cfg(feature = "nightly")] unsafe { Arc::get_mut_unchecked(arc_b) } .stage_refine(domain, &mut *container); + #[cfg(not(feature = "nightly"))] + Arc::get_mut(arc_b).unwrap() + .stage_refine(domain, &mut *container); + return Err(container) }, _ => unreachable!("This cannot happen"),
--- a/src/bisection_tree/support.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/bisection_tree/support.rs Mon Feb 03 19:22:16 2025 -0500 @@ -1,35 +1,23 @@ /*! -Traits for representing the support of a [`Apply`], and analysing the mapping on a [`Cube`]. +Traits for representing the support of a [`Mapping`], and analysing the mapping on a [`Cube`]. */ use serde::Serialize; use std::ops::{MulAssign,DivAssign,Neg}; use crate::types::{Float, Num}; use crate::maputil::map2; -use crate::mapping::Apply; +use crate::mapping::{ + Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space +}; use crate::sets::Cube; use crate::loc::Loc; use super::aggregator::Bounds; use crate::norms::{Norm, L1, L2, Linfinity}; - -/// A trait for encoding constant [`Float`] values -pub trait Constant : Copy + Sync + Send + 'static + std::fmt::Debug + Into<Self::Type> { - /// The type of the value - type Type : Float; - /// Returns the value of the constant - fn value(&self) -> Self::Type; -} +pub use crate::operator_arithmetic::{Weighted, Constant}; -impl<F : Float> Constant for F { - type Type = F; - #[inline] - fn value(&self) -> F { *self } -} - - -/// A trait for working with the supports of [`Apply`]s. +/// A trait for working with the supports of [`Mapping`]s. /// -/// Apply is not a super-trait to allow more general use. +/// `Mapping` is not a super-trait to allow more general use. pub trait Support<F : Num, const N : usize> : Sized + Sync + Send + 'static { /// Return a cube containing the support of the function represented by `self`. /// @@ -62,15 +50,9 @@ fn shift(self, x : Loc<F, N>) -> Shift<Self, F, N> { Shift { shift : x, base_fn : self } } - - /// Multiply `self` by the scalar `a`. - #[inline] - fn weigh<C : Constant<Type=F>>(self, a : C) -> Weighted<Self, C> { - Weighted { weight : a, base_fn : self } - } } -/// Trait for globally analysing a property `A` of a [`Apply`]. +/// Trait for globally analysing a property `A` of a [`Mapping`]. /// /// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as /// [`Bounds`][super::aggregator::Bounds]. @@ -91,7 +73,7 @@ // } // } -/// Trait for locally analysing a property `A` of a [`Apply`] (implementing [`Support`]) +/// Trait for locally analysing a property `A` of a [`Mapping`] (implementing [`Support`]) /// within a [`Cube`]. /// /// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as @@ -105,10 +87,10 @@ fn local_analysis(&self, cube : &Cube<F, N>) -> A; } -/// Trait for determining the upper and lower bounds of an float-valued [`Apply`]. +/// Trait for determining the upper and lower bounds of an float-valued [`Mapping`]. /// /// This is a blanket-implemented alias for [`GlobalAnalysis`]`<F, Bounds<F>>` -/// [`Apply`] is not a supertrait to allow flexibility in the implementation of either +/// [`Mapping`] is not a supertrait to allow flexibility in the implementation of either /// reference or non-reference arguments. pub trait Bounded<F : Float> : GlobalAnalysis<F, Bounds<F>> { /// Return lower and upper bounds for the values of of `self`. @@ -120,28 +102,30 @@ impl<F : Float, T : GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T { } -/// Shift of [`Support`] and [`Apply`]; output of [`Support::shift`]. +/// Shift of [`Support`] and [`Mapping`]; output of [`Support::shift`]. #[derive(Copy,Clone,Debug,Serialize)] // Serialize! but not implemented by Loc. pub struct Shift<T, F, const N : usize> { shift : Loc<F, N>, base_fn : T, } -impl<'a, T, V, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Shift<T,F,N> -where T : Apply<Loc<F, N>, Output=V> { - type Output = V; +impl<'a, T, V : Space, F : Float, const N : usize> Mapping<Loc<F, N>> for Shift<T,F,N> +where T : Mapping<Loc<F, N>, Codomain=V> { + type Codomain = V; + #[inline] - fn apply(&self, x : &'a Loc<F, N>) -> Self::Output { - self.base_fn.apply(x - &self.shift) + fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain { + self.base_fn.apply(x.own() - &self.shift) } } -impl<'a, T, V, F : Float, const N : usize> Apply<Loc<F, N>> for Shift<T,F,N> -where T : Apply<Loc<F, N>, Output=V> { - type Output = V; +impl<'a, T, V : Space, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>> for Shift<T,F,N> +where T : DifferentiableMapping<Loc<F, N>, DerivativeDomain=V> { + type Derivative = V; + #[inline] - fn apply(&self, x : Loc<F, N>) -> Self::Output { - self.base_fn.apply(x - &self.shift) + fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative { + self.base_fn.differential(x.own() - &self.shift) } } @@ -200,38 +184,6 @@ impl_shift_norm!(L1 L2 Linfinity); -/// Weighting of a [`Support`] and [`Apply`] by scalar multiplication; -/// output of [`Support::weigh`]. -#[derive(Copy,Clone,Debug,Serialize)] -pub struct Weighted<T, C : Constant> { - /// The weight - pub weight : C, - /// The base [`Support`] or [`Apply`] being weighted. - pub base_fn : T, -} - -impl<'a, T, V, F : Float, C, const N : usize> Apply<&'a Loc<F, N>> for Weighted<T, C> -where T : for<'b> Apply<&'b Loc<F, N>, Output=V>, - V : std::ops::Mul<F,Output=V>, - C : Constant<Type=F> { - type Output = V; - #[inline] - fn apply(&self, x : &'a Loc<F, N>) -> Self::Output { - self.base_fn.apply(x) * self.weight.value() - } -} - -impl<'a, T, V, F : Float, C, const N : usize> Apply<Loc<F, N>> for Weighted<T, C> -where T : Apply<Loc<F, N>, Output=V>, - V : std::ops::Mul<F,Output=V>, - C : Constant<Type=F> { - type Output = V; - #[inline] - fn apply(&self, x : Loc<F, N>) -> Self::Output { - self.base_fn.apply(x) * self.weight.value() - } -} - impl<'a, T, F : Float, C, const N : usize> Support<F,N> for Weighted<T, C> where T : Support<F, N>, C : Constant<Type=F> { @@ -331,30 +283,21 @@ impl_weighted_norm!(L1 L2 Linfinity); -/// Normalisation of [`Support`] and [`Apply`] to L¹ norm 1. +/// Normalisation of [`Support`] and [`Mapping`] to L¹ norm 1. /// /// Currently only scalar-valued functions are supported. #[derive(Copy, Clone, Debug, Serialize, PartialEq)] pub struct Normalised<T>( - /// The base [`Support`] or [`Apply`]. + /// The base [`Support`] or [`Mapping`]. pub T ); -impl<'a, T, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Normalised<T> -where T : Norm<F, L1> + for<'b> Apply<&'b Loc<F, N>, Output=F> { - type Output = F; +impl<'a, T, F : Float, const N : usize> Mapping<Loc<F, N>> for Normalised<T> +where T : Norm<F, L1> + Mapping<Loc<F,N>, Codomain=F> { + type Codomain = F; + #[inline] - fn apply(&self, x : &'a Loc<F, N>) -> Self::Output { - let w = self.0.norm(L1); - if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w } - } -} - -impl<'a, T, F : Float, const N : usize> Apply<Loc<F, N>> for Normalised<T> -where T : Norm<F, L1> + Apply<Loc<F,N>, Output=F> { - type Output = F; - #[inline] - fn apply(&self, x : Loc<F, N>) -> Self::Output { + fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain { let w = self.0.norm(L1); if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w } } @@ -449,7 +392,7 @@ : MulAssign<F> + DivAssign<F> + Neg<Output=Self> + Clone + Sync + Send + 'static { /// The identification type type Id : 'static + Copy; - /// The type of the [`Support`] (often also a [`Apply`]). + /// The type of the [`Support`] (often also a [`Mapping`]). type SupportType : 'static + Support<F, N>; /// An iterator over all the [`Support`]s of the generator. type AllDataIter<'a> : Iterator<Item=(Self::Id, Self::SupportType)> where Self : 'a;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/collection.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,54 @@ +/*! +Abstract collections of objects. +*/ + +use crate::loc::Loc; + +/// An abstract collection of elements. +pub trait Collection : IntoIterator<Item = Self::Element> { + /// Type of elements of the collection + type Element; + /// Iterator over references to elements of the collection + type RefsIter<'a> : Iterator<Item=&'a Self::Element> where Self : 'a; + + /// Returns an iterator over references to elements of the collection. + fn iter_refs(&self) -> Self::RefsIter<'_>; +} + +/// An abstract collection of mutable elements. +pub trait CollectionMut : Collection { + /// Iterator over references to elements of the collection + type RefsIterMut<'a> : Iterator<Item=&'a mut Self::Element> where Self : 'a; + + /// Returns an iterator over references to elements of the collection. + fn iter_refs_mut(&mut self) -> Self::RefsIterMut<'_>; +} + +/// Helps implement Collection and CollectionMut for slice-like collections. +#[macro_export] +macro_rules! slice_like_collection { + ($type : ty where $($where:tt)*) => { + impl<$($where)*> Collection for $type { + type Element = E; + type RefsIter<'_a> = std::slice::Iter<'_a, E> where Self : '_a; + + #[inline] + fn iter_refs(&self) -> Self::RefsIter<'_> { + self.iter() + } + } + + impl<$($where)*> CollectionMut for $type { + type RefsIterMut<'_a> = std::slice::IterMut<'_a, E> where Self : '_a; + + #[inline] + fn iter_refs_mut(&mut self) -> Self::RefsIterMut<'_> { + self.iter_mut() + } + } + } +} + +slice_like_collection!(Vec<E> where E); +slice_like_collection!([E; N] where E, const N : usize); +slice_like_collection!(Loc<E, N> where E, const N : usize);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/convex.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,290 @@ +/*! +Some convex analysis basics +*/ + +use std::marker::PhantomData; +use crate::types::*; +use crate::mapping::{Mapping, Space}; +use crate::linops::IdOp; +use crate::instance::{Instance, InstanceMut, DecompositionMut}; +use crate::operator_arithmetic::{Constant, Weighted}; +use crate::norms::*; + +/// Trait for convex mappings. Has no features, just serves as a constraint +/// +/// TODO: should constrain `Mapping::Codomain` to implement a partial order, +/// but this makes everything complicated with little benefit. +pub trait ConvexMapping<Domain : Space, F : Num = f64> : Mapping<Domain, Codomain = F> +{} + +/// Trait for mappings with a Fenchel conjugate +/// +/// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need +/// not be convex. +pub trait Conjugable<Domain : HasDual<F>, F : Num = f64> : Mapping<Domain> { + type Conjugate<'a> : ConvexMapping<Domain::DualSpace, F> where Self : 'a; + + fn conjugate(&self) -> Self::Conjugate<'_>; +} + +/// Trait for mappings with a Fenchel preconjugate +/// +/// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`], +/// but a `Preconjugable` mapping has to be convex. +pub trait Preconjugable<Domain, Predual, F : Num = f64> : ConvexMapping<Domain, F> +where + Domain : Space, + Predual : HasDual<F> +{ + type Preconjugate<'a> : Mapping<Predual> where Self : 'a; + + fn preconjugate(&self) -> Self::Preconjugate<'_>; +} + +/// Trait for mappings with a proximap map +/// +/// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need +/// not be convex. +pub trait Prox<Domain : Space> : Mapping<Domain> { + type Prox<'a> : Mapping<Domain, Codomain=Domain> where Self : 'a; + + /// Returns a proximal mapping with weight τ + fn prox_mapping(&self, τ : Self::Codomain) -> Self::Prox<'_>; + + /// Calculate the proximal mapping with weight τ + fn prox<I : Instance<Domain>>(&self, τ : Self::Codomain, z : I) -> Domain { + self.prox_mapping(τ).apply(z) + } + + /// Calculate the proximal mapping with weight τ in-place + fn prox_mut<'b>(&self, τ : Self::Codomain, y : &'b mut Domain) + where + &'b mut Domain : InstanceMut<Domain>, + Domain:: Decomp : DecompositionMut<Domain>, + for<'a> &'a Domain : Instance<Domain>, + { + *y = self.prox(τ, &*y); + } +} + + +/// Constraint to the unit ball of the norm described by `E`. +pub struct NormConstraint<F : Float, E : NormExponent> { + radius : F, + norm : NormMapping<F, E>, +} + +impl<Domain, E, F> ConvexMapping<Domain, F> for NormMapping<F, E> +where + Domain : Space, + E : NormExponent, + F : Float, + Self : Mapping<Domain, Codomain=F> +{} + + +impl<F, E, Domain> Mapping<Domain> for NormConstraint<F, E> +where + Domain : Space + Norm<F, E>, + F : Float, + E : NormExponent, +{ + type Codomain = F; + + fn apply<I : Instance<Domain>>(&self, d : I) -> F { + if d.eval(|x| x.norm(self.norm.exponent)) <= self.radius { + F::ZERO + } else { + F::INFINITY + } + } +} + +impl<Domain, E, F> ConvexMapping<Domain, F> for NormConstraint<F, E> +where + Domain : Space, + E : NormExponent, + F : Float, + Self : Mapping<Domain, Codomain=F> +{} + + +impl<E, F, Domain> Conjugable<Domain, F> for NormMapping<F, E> +where + E : HasDualExponent, + F : Float, + Domain : HasDual<F> + Norm<F, E> + Space, + <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp> +{ + type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a; + + fn conjugate(&self) -> Self::Conjugate<'_> { + NormConstraint { + radius : F::ONE, + norm : self.exponent.dual_exponent().as_mapping() + } + } +} + +impl<C, E, F, Domain> Conjugable<Domain, F> for Weighted<NormMapping<F, E>, C> +where + C : Constant<Type = F>, + E : HasDualExponent, + F : Float, + Domain : HasDual<F> + Norm<F, E> + Space, + <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp> +{ + type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a; + + fn conjugate(&self) -> Self::Conjugate<'_> { + NormConstraint { + radius : self.weight.value(), + norm : self.base_fn.exponent.dual_exponent().as_mapping() + } + } +} + +impl<Domain, E, F> Prox<Domain> for NormConstraint<F, E> +where + Domain : Space + Norm<F, E>, + E : NormExponent, + F : Float, + NormProjection<F, E> : Mapping<Domain, Codomain=Domain>, +{ + type Prox<'a> = NormProjection<F, E> where Self : 'a; + + #[inline] + fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> { + assert!(self.radius >= F::ZERO); + NormProjection{ radius : self.radius, exponent : self.norm.exponent } + } +} + +/// Projection to the unit ball of the norm described by `E`. +pub struct NormProjection<F : Float, E : NormExponent> { + radius : F, + exponent : E, +} + +/* +impl<F, Domain> Mapping<Domain> for NormProjection<F, L2> +where + Domain : Space + Euclidean<F> + std::ops::MulAssign<F>, + F : Float, +{ + type Codomain = Domain; + + fn apply<I : Instance<Domain>>(&self, d : I) -> Domain { + d.own().proj_ball2(self.radius) + } +} +*/ + +impl<F, E, Domain> Mapping<Domain> for NormProjection<F, E> +where + Domain : Space + Projection<F, E>, + F : Float, + E : NormExponent, +{ + type Codomain = Domain; + + fn apply<I : Instance<Domain>>(&self, d : I) -> Domain { + d.own().proj_ball(self.radius, self.exponent) + } +} + + +/// The zero mapping +pub struct Zero<Domain : Space, F : Num>(PhantomData<(Domain, F)>); + +impl<Domain : Space, F : Num> Zero<Domain, F> { + pub fn new() -> Self { + Zero(PhantomData) + } +} + +impl<Domain : Space, F : Num> Mapping<Domain> for Zero<Domain, F> { + type Codomain = F; + + /// Compute the value of `self` at `x`. + fn apply<I : Instance<Domain>>(&self, _x : I) -> Self::Codomain { + F::ZERO + } +} + +impl<Domain : Space, F : Num> ConvexMapping<Domain, F> for Zero<Domain, F> { } + + +impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for Zero<Domain, F> { + type Conjugate<'a> = ZeroIndicator<Domain::DualSpace, F> where Self : 'a; + + #[inline] + fn conjugate(&self) -> Self::Conjugate<'_> { + ZeroIndicator::new() + } +} + +impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F> +where + Domain : Space, + Predual : HasDual<F> +{ + type Preconjugate<'a> = ZeroIndicator<Predual, F> where Self : 'a; + + #[inline] + fn preconjugate(&self) -> Self::Preconjugate<'_> { + ZeroIndicator::new() + } +} + +impl<Domain : Space + Clone, F : Num> Prox<Domain> for Zero<Domain, F> { + type Prox<'a> = IdOp<Domain> where Self : 'a; + + #[inline] + fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> { + IdOp::new() + } +} + + +/// The zero indicator +pub struct ZeroIndicator<Domain : Space, F : Num>(PhantomData<(Domain, F)>); + +impl<Domain : Space, F : Num> ZeroIndicator<Domain, F> { + pub fn new() -> Self { + ZeroIndicator(PhantomData) + } +} + +impl<Domain : Normed<F>, F : Float> Mapping<Domain> for ZeroIndicator<Domain, F> { + type Codomain = F; + + /// Compute the value of `self` at `x`. + fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain { + x.eval(|x̃| if x̃.is_zero() { F::ZERO } else { F::INFINITY }) + } +} + +impl<Domain : Normed<F>, F : Float> ConvexMapping<Domain, F> for ZeroIndicator<Domain, F> { } + +impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for ZeroIndicator<Domain, F> { + type Conjugate<'a> = Zero<Domain::DualSpace, F> where Self : 'a; + + #[inline] + fn conjugate(&self) -> Self::Conjugate<'_> { + Zero::new() + } +} + +impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F> +where + Domain : Normed<F>, + Predual : HasDual<F> +{ + type Preconjugate<'a> = Zero<Predual, F> where Self : 'a; + + #[inline] + fn preconjugate(&self) -> Self::Preconjugate<'_> { + Zero::new() + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/direct_product.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,506 @@ +/*! +Direct products of the form $A \times B$. + +TODO: This could be easily much more generic if `derive_more` could derive arithmetic +operations on references. +*/ + +use core::ops::{Mul,MulAssign,Div,DivAssign,Add,AddAssign,Sub,SubAssign,Neg}; +use std::clone::Clone; +use serde::{Serialize, Deserialize}; +use crate::types::{Num, Float}; +use crate::{maybe_lifetime, maybe_ref}; +use crate::euclidean::Euclidean; +use crate::instance::{Instance, InstanceMut, Decomposition, DecompositionMut, MyCow}; +use crate::mapping::Space; +use crate::linops::AXPY; +use crate::loc::Loc; +use crate::norms::{Norm, PairNorm, NormExponent, Normed, HasDual, L2}; + +#[derive(Debug,Clone,Copy,PartialEq,Eq,Serialize,Deserialize)] +pub struct Pair<A, B> (pub A, pub B); + +impl<A, B> Pair<A,B> { + pub fn new(a : A, b : B) -> Pair<A,B> { Pair(a, b) } +} + +impl<A, B> From<(A,B)> for Pair<A,B> { + #[inline] + fn from((a, b) : (A, B)) -> Pair<A,B> { Pair(a, b) } +} + +impl<A, B> From<Pair<A,B>> for (A,B) { + #[inline] + fn from(Pair(a, b) : Pair<A, B>) -> (A,B) { (a, b) } +} + +macro_rules! impl_binop { + (($a : ty, $b : ty), $trait : ident, $fn : ident, $refl:ident, $refr:ident) => { + impl_binop!(@doit: $a, $b, $trait, $fn; + maybe_lifetime!($refl, &'l Pair<$a,$b>), + (maybe_lifetime!($refl, &'l $a), + maybe_lifetime!($refl, &'l $b)); + maybe_lifetime!($refr, &'r Pair<Ai,Bi>), + (maybe_lifetime!($refr, &'r Ai), + maybe_lifetime!($refr, &'r Bi)); + $refl, $refr); + }; + + (@doit: $a:ty, $b:ty, + $trait:ident, $fn:ident; + $self:ty, ($aself:ty, $bself:ty); + $in:ty, ($ain:ty, $bin:ty); + $refl:ident, $refr:ident) => { + impl<'l, 'r, Ai, Bi> $trait<$in> + for $self + where $aself: $trait<$ain>, + $bself: $trait<$bin> { + type Output = Pair<<$aself as $trait<$ain>>::Output, + <$bself as $trait<$bin>>::Output>; + + #[inline] + fn $fn(self, y : $in) -> Self::Output { + Pair(maybe_ref!($refl, self.0).$fn(maybe_ref!($refr, y.0)), + maybe_ref!($refl, self.1).$fn(maybe_ref!($refr, y.1))) + } + } + }; +} + +macro_rules! impl_assignop { + (($a : ty, $b : ty), $trait : ident, $fn : ident, $refr:ident) => { + impl_assignop!(@doit: $a, $b, + $trait, $fn; + maybe_lifetime!($refr, &'r Pair<Ai,Bi>), + (maybe_lifetime!($refr, &'r Ai), + maybe_lifetime!($refr, &'r Bi)); + $refr); + }; + (@doit: $a : ty, $b : ty, + $trait:ident, $fn:ident; + $in:ty, ($ain:ty, $bin:ty); + $refr:ident) => { + impl<'r, Ai, Bi> $trait<$in> + for Pair<$a,$b> + where $a: $trait<$ain>, + $b: $trait<$bin> { + #[inline] + fn $fn(&mut self, y : $in) -> () { + self.0.$fn(maybe_ref!($refr, y.0)); + self.1.$fn(maybe_ref!($refr, y.1)); + } + } + } +} + +macro_rules! impl_scalarop { + (($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident, $refl:ident) => { + impl_scalarop!(@doit: $field, + $trait, $fn; + maybe_lifetime!($refl, &'l Pair<$a,$b>), + (maybe_lifetime!($refl, &'l $a), + maybe_lifetime!($refl, &'l $b)); + $refl); + }; + (@doit: $field : ty, + $trait:ident, $fn:ident; + $self:ty, ($aself:ty, $bself:ty); + $refl:ident) => { + // Scalar as Rhs + impl<'l> $trait<$field> + for $self + where $aself: $trait<$field>, + $bself: $trait<$field> { + type Output = Pair<<$aself as $trait<$field>>::Output, + <$bself as $trait<$field>>::Output>; + #[inline] + fn $fn(self, a : $field) -> Self::Output { + Pair(maybe_ref!($refl, self.0).$fn(a), + maybe_ref!($refl, self.1).$fn(a)) + } + } + } +} + +// Not used due to compiler overflow +#[allow(unused_macros)] +macro_rules! impl_scalarlhs_op { + (($a : ty, $b : ty), $field : ty, $trait:ident, $fn:ident, $refr:ident) => { + impl_scalarlhs_op!(@doit: $trait, $fn, + maybe_lifetime!($refr, &'r Pair<$a,$b>), + (maybe_lifetime!($refr, &'r $a), + maybe_lifetime!($refr, &'r $b)); + $refr, $field); + }; + (@doit: $trait:ident, $fn:ident, + $in:ty, ($ain:ty, $bin:ty); + $refr:ident, $field:ty) => { + impl<'r> $trait<$in> + for $field + where $field : $trait<$ain> + + $trait<$bin> { + type Output = Pair<<$field as $trait<$ain>>::Output, + <$field as $trait<$bin>>::Output>; + #[inline] + fn $fn(self, x : $in) -> Self::Output { + Pair(self.$fn(maybe_ref!($refr, x.0)), + self.$fn(maybe_ref!($refr, x.1))) + } + } + }; +} + +macro_rules! impl_scalar_assignop { + (($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => { + impl<'r> $trait<$field> + for Pair<$a, $b> + where $a: $trait<$field>, $b: $trait<$field> { + #[inline] + fn $fn(&mut self, a : $field) -> () { + self.0.$fn(a); + self.1.$fn(a); + } + } + } +} + +macro_rules! impl_unaryop { + (($a : ty, $b : ty), $trait:ident, $fn:ident, $refl:ident) => { + impl_unaryop!(@doit: $trait, $fn; + maybe_lifetime!($refl, &'l Pair<$a,$b>), + (maybe_lifetime!($refl, &'l $a), + maybe_lifetime!($refl, &'l $b)); + $refl); + }; + (@doit: $trait:ident, $fn:ident; + $self:ty, ($aself:ty, $bself:ty); + $refl : ident) => { + impl<'l> $trait + for $self + where $aself: $trait, + $bself: $trait { + type Output = Pair<<$aself as $trait>::Output, + <$bself as $trait>::Output>; + #[inline] + fn $fn(self) -> Self::Output { + Pair(maybe_ref!($refl, self.0).$fn(), + maybe_ref!($refl, self.1).$fn()) + } + } + } +} + +#[macro_export] +macro_rules! impl_pair_vectorspace_ops { + (($a:ty, $b:ty), $field:ty) => { + impl_pair_vectorspace_ops!(@binary, ($a, $b), Add, add); + impl_pair_vectorspace_ops!(@binary, ($a, $b), Sub, sub); + impl_pair_vectorspace_ops!(@assign, ($a, $b), AddAssign, add_assign); + impl_pair_vectorspace_ops!(@assign, ($a, $b), SubAssign, sub_assign); + impl_pair_vectorspace_ops!(@scalar, ($a, $b), $field, Mul, mul); + impl_pair_vectorspace_ops!(@scalar, ($a, $b), $field, Div, div); + // Compiler overflow + // $( + // impl_pair_vectorspace_ops!(@scalar_lhs, ($a, $b), $field, $impl_scalarlhs_op, Mul, mul); + // )* + impl_pair_vectorspace_ops!(@scalar_assign, ($a, $b), $field, MulAssign, mul_assign); + impl_pair_vectorspace_ops!(@scalar_assign, ($a, $b), $field, DivAssign, div_assign); + impl_pair_vectorspace_ops!(@unary, ($a, $b), Neg, neg); + }; + (@binary, ($a : ty, $b : ty), $trait : ident, $fn : ident) => { + impl_binop!(($a, $b), $trait, $fn, ref, ref); + impl_binop!(($a, $b), $trait, $fn, ref, noref); + impl_binop!(($a, $b), $trait, $fn, noref, ref); + impl_binop!(($a, $b), $trait, $fn, noref, noref); + }; + (@assign, ($a : ty, $b : ty), $trait : ident, $fn :ident) => { + impl_assignop!(($a, $b), $trait, $fn, ref); + impl_assignop!(($a, $b), $trait, $fn, noref); + }; + (@scalar, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn :ident) => { + impl_scalarop!(($a, $b), $field, $trait, $fn, ref); + impl_scalarop!(($a, $b), $field, $trait, $fn, noref); + }; + (@scalar_lhs, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => { + impl_scalarlhs_op!(($a, $b), $field, $trait, $fn, ref); + impl_scalarlhs_op!(($a, $b), $field, $trait, $fn, noref); + }; + (@scalar_assign, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => { + impl_scalar_assignop!(($a, $b), $field, $trait, $fn); + }; + (@unary, ($a : ty, $b : ty), $trait : ident, $fn : ident) => { + impl_unaryop!(($a, $b), $trait, $fn, ref); + impl_unaryop!(($a, $b), $trait, $fn, noref); + }; +} + +impl_pair_vectorspace_ops!((f32, f32), f32); +impl_pair_vectorspace_ops!((f64, f64), f64); + +type PairOutput<F, A, B> = Pair<<A as Euclidean<F>>::Output, <B as Euclidean<F>>::Output>; + +impl<A, B, F> Euclidean<F> for Pair<A, B> +where + A : Euclidean<F>, + B : Euclidean<F>, + F : Float, + PairOutput<F, A, B> : Euclidean<F>, + Self : Sized + + Mul<F, Output=PairOutput<F, A, B>> + MulAssign<F> + + Div<F, Output=PairOutput<F, A, B>> + DivAssign<F> + + Add<Self, Output=PairOutput<F, A, B>> + + Sub<Self, Output=PairOutput<F, A, B>> + + for<'b> Add<&'b Self, Output=PairOutput<F, A, B>> + + for<'b> Sub<&'b Self, Output=PairOutput<F, A, B>> + + AddAssign<Self> + for<'b> AddAssign<&'b Self> + + SubAssign<Self> + for<'b> SubAssign<&'b Self> + + Neg<Output=PairOutput<F, A, B>> +{ + type Output = PairOutput<F, A, B>; + + fn dot<I : Instance<Self>>(&self, other : I) -> F { + let Pair(u, v) = other.decompose(); + self.0.dot(u) + self.1.dot(v) + } + + fn norm2_squared(&self) -> F { + self.0.norm2_squared() + self.1.norm2_squared() + } + + fn dist2_squared<I : Instance<Self>>(&self, other : I) -> F { + let Pair(u, v) = other.decompose(); + self.0.dist2_squared(u) + self.1.dist2_squared(v) + } +} + +impl<F, A, B, U, V> AXPY<F, Pair<U, V>> for Pair<A, B> +where + U : Space, + V : Space, + A : AXPY<F, U>, + B : AXPY<F, V>, + F : Num, + Self : MulAssign<F>, + Pair<A, B> : MulAssign<F>, + Pair<A::Owned, B::Owned> : AXPY<F, Pair<U, V>>, +{ + + type Owned = Pair<A::Owned, B::Owned>; + + fn axpy<I : Instance<Pair<U,V>>>(&mut self, α : F, x : I, β : F) { + let Pair(u, v) = x.decompose(); + self.0.axpy(α, u, β); + self.1.axpy(α, v, β); + } + + fn copy_from<I : Instance<Pair<U,V>>>(&mut self, x : I) { + let Pair(u, v) = x.decompose(); + self.0.copy_from(u); + self.1.copy_from(v); + } + + fn scale_from<I : Instance<Pair<U,V>>>(&mut self, α : F, x : I) { + let Pair(u, v) = x.decompose(); + self.0.scale_from(α, u); + self.1.scale_from(α, v); + } + + /// Return a similar zero as `self`. + fn similar_origin(&self) -> Self::Owned { + Pair(self.0.similar_origin(), self.1.similar_origin()) + } + + /// Set self to zero. + fn set_zero(&mut self) { + self.0.set_zero(); + self.1.set_zero(); + } +} + +/// [`Decomposition`] for working with [`Pair`]s. +#[derive(Copy, Clone, Debug)] +pub struct PairDecomposition<D, Q>(D, Q); + +impl<A : Space, B : Space> Space for Pair<A, B> { + type Decomp = PairDecomposition<A::Decomp, B::Decomp>; +} + +impl<A, B, D, Q> Decomposition<Pair<A, B>> for PairDecomposition<D,Q> +where + A : Space, + B : Space, + D : Decomposition<A>, + Q : Decomposition<B>, +{ + type Decomposition<'b> = Pair<D::Decomposition<'b>, Q::Decomposition<'b>> where Pair<A, B> : 'b; + type Reference<'b> = Pair<D::Reference<'b>, Q::Reference<'b>> where Pair<A, B> : 'b; + + #[inline] + fn lift<'b>(Pair(u, v) : Self::Reference<'b>) -> Self::Decomposition<'b> { + Pair(D::lift(u), Q::lift(v)) + } +} + +impl<A, B, U, V, D, Q> Instance<Pair<A, B>, PairDecomposition<D, Q>> for Pair<U, V> +where + A : Space, + B : Space, + D : Decomposition<A>, + Q : Decomposition<B>, + U : Instance<A, D>, + V : Instance<B, Q>, +{ + #[inline] + fn decompose<'b>(self) + -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Decomposition<'b> + where Self : 'b, Pair<A, B> : 'b + { + Pair(self.0.decompose(), self.1.decompose()) + } + + #[inline] + fn ref_instance(&self) + -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Reference<'_> + { + Pair(self.0.ref_instance(), self.1.ref_instance()) + } + + #[inline] + fn cow<'b>(self) -> MyCow<'b, Pair<A, B>> where Self : 'b{ + MyCow::Owned(Pair(self.0.own(), self.1.own())) + } + + #[inline] + fn own(self) -> Pair<A, B> { + Pair(self.0.own(), self.1.own()) + } +} + + +impl<'a, A, B, U, V, D, Q> Instance<Pair<A, B>, PairDecomposition<D, Q>> for &'a Pair<U, V> +where + A : Space, + B : Space, + D : Decomposition<A>, + Q : Decomposition<B>, + U : Instance<A, D>, + V : Instance<B, Q>, + &'a U : Instance<A, D>, + &'a V : Instance<B, Q>, +{ + #[inline] + fn decompose<'b>(self) + -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Decomposition<'b> + where Self : 'b, Pair<A, B> : 'b + { + Pair(D::lift(self.0.ref_instance()), Q::lift(self.1.ref_instance())) + } + + #[inline] + fn ref_instance(&self) + -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Reference<'_> + { + Pair(self.0.ref_instance(), self.1.ref_instance()) + } + + #[inline] + fn cow<'b>(self) -> MyCow<'b, Pair<A, B>> where Self : 'b { + MyCow::Owned(self.own()) + } + + #[inline] + fn own(self) -> Pair<A, B> { + let Pair(ref u, ref v) = self; + Pair(u.own(), v.own()) + } + +} + +impl<A, B, D, Q> DecompositionMut<Pair<A, B>> for PairDecomposition<D,Q> +where + A : Space, + B : Space, + D : DecompositionMut<A>, + Q : DecompositionMut<B>, +{ + type ReferenceMut<'b> = Pair<D::ReferenceMut<'b>, Q::ReferenceMut<'b>> where Pair<A, B> : 'b; +} + +impl<A, B, U, V, D, Q> InstanceMut<Pair<A, B>, PairDecomposition<D, Q>> for Pair<U, V> +where + A : Space, + B : Space, + D : DecompositionMut<A>, + Q : DecompositionMut<B>, + U : InstanceMut<A, D>, + V : InstanceMut<B, Q>, +{ + #[inline] + fn ref_instance_mut(&mut self) + -> <PairDecomposition<D, Q> as DecompositionMut<Pair<A, B>>>::ReferenceMut<'_> + { + Pair(self.0.ref_instance_mut(), self.1.ref_instance_mut()) + } +} + +impl<'a, A, B, U, V, D, Q> InstanceMut<Pair<A, B>, PairDecomposition<D, Q>> for &'a mut Pair<U, V> +where + A : Space, + B : Space, + D : DecompositionMut<A>, + Q : DecompositionMut<B>, + U : InstanceMut<A, D>, + V : InstanceMut<B, Q>, +{ + #[inline] + fn ref_instance_mut(&mut self) + -> <PairDecomposition<D, Q> as DecompositionMut<Pair<A, B>>>::ReferenceMut<'_> + { + Pair(self.0.ref_instance_mut(), self.1.ref_instance_mut()) + } +} + + +impl<F, A, B, ExpA, ExpB, ExpJ> Norm<F, PairNorm<ExpA, ExpB, ExpJ>> +for Pair<A,B> +where + F : Num, + ExpA : NormExponent, + ExpB : NormExponent, + ExpJ : NormExponent, + A : Norm<F, ExpA>, + B : Norm<F, ExpB>, + Loc<F, 2> : Norm<F, ExpJ>, +{ + fn norm(&self, PairNorm(expa, expb, expj) : PairNorm<ExpA, ExpB, ExpJ>) -> F { + Loc([self.0.norm(expa), self.1.norm(expb)]).norm(expj) + } +} + + +impl<F : Float, A, B> Normed<F> for Pair<A,B> +where + A : Normed<F>, + B : Normed<F>, +{ + type NormExp = PairNorm<A::NormExp, B::NormExp, L2>; + + #[inline] + fn norm_exponent(&self) -> Self::NormExp { + PairNorm(self.0.norm_exponent(), self.1.norm_exponent(), L2) + } + + #[inline] + fn is_zero(&self) -> bool { + self.0.is_zero() && self.1.is_zero() + } +} + +impl<F : Float, A, B> HasDual<F> for Pair<A,B> +where + A : HasDual<F>, + B : HasDual<F>, + +{ + type DualSpace = Pair<A::DualSpace, B::DualSpace>; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/discrete_gradient.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,489 @@ +/*! +Simple disrete gradient operators + */ +use numeric_literals::replace_float_literals; +use nalgebra::{ + DVector, Matrix, U1, Storage, StorageMut, Dyn +}; +use crate::types::Float; +use crate::instance::Instance; +use crate::linops::{Mapping, Linear, BoundedLinear, Adjointable, GEMV}; +use crate::norms::{Norm, L2}; + +#[derive(Copy, Clone, Debug)] +/// Forward differences with Neumann boundary conditions +pub struct ForwardNeumann; + +#[derive(Copy, Clone, Debug)] +/// Forward differences with Dirichlet boundary conditions +pub struct ForwardDirichlet; + +#[derive(Copy, Clone, Debug)] +/// Backward differences with Dirichlet boundary conditions +pub struct BackwardDirichlet; + +#[derive(Copy, Clone, Debug)] +/// Backward differences with Neumann boundary conditions +pub struct BackwardNeumann; + +/// Finite differences gradient +pub struct Grad< + F : Float + nalgebra::RealField, + B : Discretisation<F>, + const N : usize +> { + dims : [usize; N], + h : F, // may be negative to implement adjoints! + discretisation : B, +} + + +/// Finite differences divergence +pub struct Div< + F : Float + nalgebra::RealField, + B : Discretisation<F>, + const N : usize +> { + dims : [usize; N], + h : F, // may be negative to implement adjoints! + discretisation : B, +} + +/// Internal: classification of a point in a 1D discretisation +pub enum DiscretisationOrInterior { + /// center, forward + LeftBoundary(usize, usize), + /// center, backward + RightBoundary(usize, usize), + /// center, (backward, forward) + Interior(usize, (usize, usize)), +} + +use DiscretisationOrInterior::*; + +/// Trait for different discretisations +pub trait Discretisation<F : Float + nalgebra::RealField> : Copy { + /// Opposite discretisation, appropriate for adjoints with negated cell width. + type Opposite : Discretisation<F>; + + /// Add to appropiate index of `v` (as determined by `b`) the appropriate difference + /// of `x` with cell width `h`. + fn add_diff_mut<SMut, S>( + &self, + v : &mut Matrix<F, Dyn, U1, SMut>, + x : &Matrix<F, Dyn, U1, S>, + α : F, + b : DiscretisationOrInterior, + ) where + SMut : StorageMut<F, Dyn, U1>, + S : Storage<F, Dyn, U1>; + + /// Give the opposite discretisation, appropriate for adjoints with negated `h`. + fn opposite(&self) -> Self::Opposite; + + /// Bound for the corresponding operator norm. + #[replace_float_literals(F::cast_from(literal))] + fn opnorm_bound(&self, h : F) -> F { + // See: Chambolle, “An Algorithm for Total Variation Minimization and Applications”. + // Ok for forward and backward differences. + // + // Fuck nalgebra for polluting everything with its own shit. + num_traits::Float::sqrt(8.0) / h + } +} + +impl<F : Float + nalgebra::RealField> Discretisation<F> for ForwardNeumann { + type Opposite = BackwardDirichlet; + + #[inline] + fn add_diff_mut<SMut, S>( + &self, + v : &mut Matrix<F, Dyn, U1, SMut>, + x : &Matrix<F, Dyn, U1, S>, + α : F, + b : DiscretisationOrInterior, + ) where + SMut : StorageMut<F, Dyn, U1>, + S : Storage<F, Dyn, U1> + { + match b { + Interior(c, (_, f)) | LeftBoundary(c, f) => { v[c] += (x[f] - x[c]) * α }, + RightBoundary(_c, _b) => { }, + } + } + + #[inline] + fn opposite(&self) -> Self::Opposite { + BackwardDirichlet + } +} + +impl<F : Float + nalgebra::RealField> Discretisation<F> for BackwardNeumann { + type Opposite = ForwardDirichlet; + + #[inline] + fn add_diff_mut<SMut, S>( + &self, + v : &mut Matrix<F, Dyn, U1, SMut>, + x : &Matrix<F, Dyn, U1, S>, + α : F, + b : DiscretisationOrInterior, + ) where + SMut : StorageMut<F, Dyn, U1>, + S : Storage<F, Dyn, U1> + { + match b { + Interior(c, (b, _)) | RightBoundary(c, b) => { v[c] += (x[c] - x[b]) * α }, + LeftBoundary(_c, _f) => { }, + } + } + + #[inline] + fn opposite(&self) -> Self::Opposite { + ForwardDirichlet + } +} + +impl<F : Float + nalgebra::RealField> Discretisation<F> for BackwardDirichlet { + type Opposite = ForwardNeumann; + + #[inline] + fn add_diff_mut<SMut, S>( + &self, + v : &mut Matrix<F, Dyn, U1, SMut>, + x : &Matrix<F, Dyn, U1, S>, + α : F, + b : DiscretisationOrInterior, + ) where + SMut : StorageMut<F, Dyn, U1>, + S : Storage<F, Dyn, U1> + { + match b { + Interior(c, (b, _f)) => { v[c] += (x[c] - x[b]) * α }, + LeftBoundary(c, _f) => { v[c] += x[c] * α }, + RightBoundary(c, b) => { v[c] -= x[b] * α }, + } + } + + #[inline] + fn opposite(&self) -> Self::Opposite { + ForwardNeumann + } +} + +impl<F : Float + nalgebra::RealField> Discretisation<F> for ForwardDirichlet { + type Opposite = BackwardNeumann; + + #[inline] + fn add_diff_mut<SMut, S>( + &self, + v : &mut Matrix<F, Dyn, U1, SMut>, + x : &Matrix<F, Dyn, U1, S>, + α : F, + b : DiscretisationOrInterior, + ) where + SMut : StorageMut<F, Dyn, U1>, + S : Storage<F, Dyn, U1> + { + match b { + Interior(c, (_b, f)) => { v[c] += (x[f] - x[c]) * α }, + LeftBoundary(c, f) => { v[c] += x[f] * α }, + RightBoundary(c, _b) => { v[c] -= x[c] * α }, + } + } + + #[inline] + fn opposite(&self) -> Self::Opposite { + BackwardNeumann + } +} + +struct Iter<'a, const N : usize> { + /// Dimensions + dims : &'a [usize; N], + /// Dimension along which to calculate differences + d : usize, + /// Stride along coordinate d + d_stride : usize, + /// Cartesian indices + i : [usize; N], + /// Linear index + k : usize, + /// Maximal linear index + len : usize +} + +impl<'a, const N : usize> Iter<'a, N> { + fn new(dims : &'a [usize; N], d : usize) -> Self { + let d_stride = dims[0..d].iter().product::<usize>(); + let len = dims.iter().product::<usize>(); + Iter{ dims, d, d_stride, i : [0; N], k : 0, len } + } +} + +impl<'a, const N : usize> Iterator for Iter<'a, N> { + type Item = DiscretisationOrInterior; + fn next(&mut self) -> Option<Self::Item> { + let res = if self.k >= self.len { + None + } else { + let cartesian_idx = self.i[self.d]; + let cartesian_max = self.dims[self.d]; + let k = self.k; + + if cartesian_idx == 0 { + Some(LeftBoundary(k, k + self.d_stride)) + } else if cartesian_idx + 1 >= cartesian_max { + Some(RightBoundary(k, k - self.d_stride)) + } else { + Some(Interior(k, (k - self.d_stride, k + self.d_stride))) + } + }; + self.k += 1; + for j in 0..N { + if self.i[j] + 1 < self.dims[j] { + self.i[j] += 1; + break + } + self.i[j] = 0 + } + res + } +} + +impl<F, B, const N : usize> Mapping<DVector<F>> +for Grad<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ + type Codomain = DVector<F>; + fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> { + let mut y = DVector::zeros(N * self.len()); + self.apply_add(&mut y, i); + y + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F, B, const N : usize> GEMV<F, DVector<F>> +for Grad<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ + fn gemv<I : Instance<DVector<F>>>( + &self, y : &mut DVector<F>, α : F, i : I, β : F + ) { + if β == 0.0 { + y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0); + } else if β != 1.0 { + //*y *= β; + y.as_mut_slice().iter_mut().for_each(|x| *x *= β); + } + let h = self.h; + let m = self.len(); + i.eval(|x| { + assert_eq!(x.len(), m); + for d in 0..N { + let mut v = y.generic_view_mut((d*m, 0), (Dyn(m), U1)); + for b in Iter::new(&self.dims, d) { + self.discretisation.add_diff_mut(&mut v, x, α/h, b) + } + } + }) + } +} + +impl<F, B, const N : usize> Mapping<DVector<F>> +for Div<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ + type Codomain = DVector<F>; + fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> { + let mut y = DVector::zeros(self.len()); + self.apply_add(&mut y, i); + y + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F, B, const N : usize> GEMV<F, DVector<F>> +for Div<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ + fn gemv<I : Instance<DVector<F>>>( + &self, y : &mut DVector<F>, α : F, i : I, β : F + ) { + if β == 0.0 { + y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0); + } else if β != 1.0 { + //*y *= β; + y.as_mut_slice().iter_mut().for_each(|x| *x *= β); + } + let h = self.h; + let m = self.len(); + i.eval(|x| { + assert_eq!(x.len(), N * m); + for d in 0..N { + let v = x.generic_view((d*m, 0), (Dyn(m), U1)); + for b in Iter::new(&self.dims, d) { + self.discretisation.add_diff_mut(y, &v, α/h, b) + } + } + }) + } +} + +impl<F, B, const N : usize> Grad<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField +{ + /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions. + pub fn new_for(u : &DVector<F>, h : F, dims : [usize; N], discretisation : B) + -> Option<Self> + { + if u.len() == dims.iter().product::<usize>() { + Some(Grad { dims, h, discretisation } ) + } else { + None + } + } + + fn len(&self) -> usize { + self.dims.iter().product::<usize>() + } +} + + +impl<F, B, const N : usize> Div<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField +{ + /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions. + pub fn new_for(u : &DVector<F>, h : F, dims : [usize; N], discretisation : B) + -> Option<Self> + { + if u.len() == dims.iter().product::<usize>() * N { + Some(Div { dims, h, discretisation } ) + } else { + None + } + } + + fn len(&self) -> usize { + self.dims.iter().product::<usize>() + } +} + +impl<F, B, const N : usize> Linear<DVector<F>> +for Grad<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ +} + +impl<F, B, const N : usize> Linear<DVector<F>> +for Div<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ +} + +impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F> +for Grad<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, + DVector<F> : Norm<F, L2>, +{ + fn opnorm_bound(&self, _ : L2, _ : L2) -> F { + // Fuck nalgebra. + self.discretisation.opnorm_bound(num_traits::Float::abs(self.h)) + } +} + + +impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F> +for Div<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, + DVector<F> : Norm<F, L2>, +{ + fn opnorm_bound(&self, _ : L2, _ : L2) -> F { + // Fuck nalgebra. + self.discretisation.opnorm_bound(num_traits::Float::abs(self.h)) + } +} + +impl<F, B, const N : usize> +Adjointable<DVector<F>, DVector<F>> +for Grad<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ + type AdjointCodomain = DVector<F>; + type Adjoint<'a> = Div<F, B::Opposite, N> where Self : 'a; + + /// Form the adjoint operator of `self`. + fn adjoint(&self) -> Self::Adjoint<'_> { + Div { + dims : self.dims, + h : -self.h, + discretisation : self.discretisation.opposite(), + } + } +} + + +impl<F, B, const N : usize> +Adjointable<DVector<F>, DVector<F>> +for Div<F, B, N> +where + B : Discretisation<F>, + F : Float + nalgebra::RealField, +{ + type AdjointCodomain = DVector<F>; + type Adjoint<'a> = Grad<F, B::Opposite, N> where Self : 'a; + + /// Form the adjoint operator of `self`. + fn adjoint(&self) -> Self::Adjoint<'_> { + Grad { + dims : self.dims, + h : -self.h, + discretisation : self.discretisation.opposite(), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn grad_adjoint() { + let im = DVector::from( (0..9).map(|t| t as f64).collect::<Vec<_>>()); + let v = DVector::from( (0..18).map(|t| t as f64).collect::<Vec<_>>()); + + let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardNeumann).unwrap(); + let a = grad.apply(&im).dot(&v); + let b = grad.adjoint().apply(&v).dot(&im); + assert_eq!(a, b); + + let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardDirichlet).unwrap(); + let a = grad.apply(&im).dot(&v); + let b = grad.adjoint().apply(&v).dot(&im); + assert_eq!(a, b); + + } +}
--- a/src/error.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/error.rs Mon Feb 03 19:22:16 2025 -0500 @@ -2,10 +2,8 @@ Error passing helper types */ -use std::error::Error; - /// A [`Result`] containing `T` or a dynamic error type -pub type DynResult<T> = Result<T, Box<dyn Error>>; +pub type DynResult<T> = Result<T, anyhow::Error>; /// A [`Result`] containing `()` or a dynamic error type pub type DynError = DynResult<()>;
--- a/src/euclidean.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/euclidean.rs Mon Feb 03 19:22:16 2025 -0500 @@ -2,42 +2,37 @@ Euclidean spaces. */ +use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg}; use crate::types::*; -use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg}; - -/// Space (type) with a defined dot product. -/// -/// `U` is the space of the multiplier, and `F` the space of scalars. -/// Since `U` ≠ `Self`, this trait can also implement dual products. -pub trait Dot<U, F> { - fn dot(&self, other : &U) -> F; -} +use crate::instance::Instance; +use crate::norms::{HasDual, Reflexive}; /// Space (type) with Euclidean and vector space structure /// /// The type should implement vector space operations (addition, subtraction, scalar /// multiplication and scalar division) along with their assignment versions, as well -/// as the [`Dot`] product with respect to `Self`. -pub trait Euclidean<F : Float> : Sized + Dot<Self,F> - + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F> - + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F> - + Add<Self, Output=<Self as Euclidean<F>>::Output> - + Sub<Self, Output=<Self as Euclidean<F>>::Output> - + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output> - + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output> - + AddAssign<Self> + for<'b> AddAssign<&'b Self> - + SubAssign<Self> + for<'b> SubAssign<&'b Self> - + Neg<Output=<Self as Euclidean<F>>::Output> { +/// as an inner product. +pub trait Euclidean<F : Float> : HasDual<F, DualSpace=Self> + Reflexive<F> + + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F> + + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F> + + Add<Self, Output=<Self as Euclidean<F>>::Output> + + Sub<Self, Output=<Self as Euclidean<F>>::Output> + + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output> + + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output> + + AddAssign<Self> + for<'b> AddAssign<&'b Self> + + SubAssign<Self> + for<'b> SubAssign<&'b Self> + + Neg<Output=<Self as Euclidean<F>>::Output> +{ type Output : Euclidean<F>; - /// Returns origin of same dimensions as `self`. - fn similar_origin(&self) -> <Self as Euclidean<F>>::Output; + // Inner product + fn dot<I : Instance<Self>>(&self, other : I) -> F; /// Calculate the square of the 2-norm, $\frac{1}{2}\\|x\\|_2^2$, where `self` is $x$. - #[inline] - fn norm2_squared(&self) -> F { - self.dot(self) - } + /// + /// This is not automatically implemented to avoid imposing + /// `for <'a> &'a Self : Instance<Self>` trait bound bloat. + fn norm2_squared(&self) -> F; /// Calculate the square of the 2-norm divided by 2, $\frac{1}{2}\\|x\\|_2^2$, /// where `self` is $x$. @@ -53,11 +48,11 @@ } /// Calculate the 2-distance squared $\\|x-y\\|_2^2$, where `self` is $x$. - fn dist2_squared(&self, y : &Self) -> F; + fn dist2_squared<I : Instance<Self>>(&self, y : I) -> F; /// Calculate the 2-distance $\\|x-y\\|_2$, where `self` is $x$. #[inline] - fn dist2(&self, y : &Self) -> F { + fn dist2<I : Instance<Self>>(&self, y : I) -> F { self.dist2_squared(y).sqrt() }
--- a/src/fe_model/p2_local_model.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/fe_model/p2_local_model.rs Mon Feb 03 19:22:16 2025 -0500 @@ -6,7 +6,8 @@ use crate::loc::Loc; use crate::sets::{Set,NPolygon,SpannedHalfspace}; use crate::linsolve::*; -use crate::euclidean::Dot; +use crate::euclidean::Euclidean; +use crate::instance::Instance; use super::base::{LocalModel,RealLocalModel}; use crate::sets::Cube; use numeric_literals::replace_float_literals; @@ -30,7 +31,8 @@ impl<'a, F : Float> Set<Loc<F,1>> for RealInterval<F> { #[inline] - fn contains(&self, &Loc([x]) : &Loc<F,1>) -> bool { + fn contains<I : Instance<Loc<F, 1>>>(&self, z : I) -> bool { + let &Loc([x]) = z.ref_instance(); let &[Loc([x0]), Loc([x1])] = &self.0; (x0 < x && x < x1) || (x1 < x && x < x0) } @@ -38,11 +40,11 @@ impl<'a, F : Float> Set<Loc<F,2>> for PlanarSimplex<F> { #[inline] - fn contains(&self, x : &Loc<F,2>) -> bool { + fn contains<I : Instance<Loc<F, 2>>>(&self, z : I) -> bool { let &[x0, x1, x2] = &self.0; NPolygon([[x0, x1].spanned_halfspace(), [x1, x2].spanned_halfspace(), - [x2, x0].spanned_halfspace()]).contains(x) + [x2, x0].spanned_halfspace()]).contains(z) } }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/instance.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,282 @@ +/*! +Helper traits to work with references or owned values of types and their subsets. +*/ + +#[derive(Clone, Copy)] +pub enum EitherDecomp<A, B> { + Owned(A), + Borrowed(B), +} + +/// A very basic implementation of [`std::borrow::Cow`] without a [`Clone`] trait dependency. +pub type MyCow<'b, X> = EitherDecomp<X, &'b X>; + +impl<'b, X> std::ops::Deref for MyCow<'b, X> { + type Target = X; + + #[inline] + fn deref(&self) -> &Self::Target { + match self { + EitherDecomp::Owned(x) => &x, + EitherDecomp::Borrowed(x) => x, + } + } +} + +impl<'b, X> MyCow<'b, X> { + #[inline] + pub fn into_owned(self) -> X where X : Clone { + match self { + EitherDecomp::Owned(x) => x, + EitherDecomp::Borrowed(x) => x.clone(), + } + } +} + +/// Trait for abitrary mathematical spaces. +pub trait Space : Instance<Self, Self::Decomp> { + /// Default decomposition for the space + type Decomp : Decomposition<Self>; +} + +#[macro_export] +macro_rules! impl_basic_space { + ($($type:ty)*) => { $( + impl $crate::instance::Space for $type { + type Decomp = $crate::instance::BasicDecomposition; + } + )* }; + ($type:ty where $($where:tt)*) => { + impl<$($where)*> $crate::instance::Space for $type { + type Decomp = $crate::instance::BasicDecomposition; + } + }; +} + +impl_basic_space!(u8 u16 u32 u64 u128 usize + i8 i16 i32 i64 i128 isize + f32 f64); + +/// Marker type for decompositions to be used with [`Instance`]. +pub trait Decomposition<X : Space> : Sized { + /// Possibly owned form of the decomposition + type Decomposition<'b> : Instance<X, Self> where X : 'b; + /// Unlikely owned form of the decomposition. + /// Type for a lightweight intermediate conversion that does not own the original variable. + /// Usually this is just a reference, but may also be a lightweight structure that + /// contains references; see the implementation for [`crate::direct_product::Pair`]. + type Reference<'b> : Instance<X, Self> + Copy where X : 'b; + + /// Left the lightweight reference type into a full decomposition type. + fn lift<'b>(r : Self::Reference<'b>) -> Self::Decomposition<'b>; +} + +/// Most common [`Decomposition`] (into `Either<X, &'b X>`) that allows working with owned +/// values and all sorts of references. +#[derive(Copy, Clone, Debug)] +pub struct BasicDecomposition; + +impl<X : Space + Clone> Decomposition<X> for BasicDecomposition { + type Decomposition<'b> = MyCow<'b, X> where X : 'b; + type Reference<'b> = &'b X where X : 'b; + + #[inline] + fn lift<'b>(r : Self::Reference<'b>) -> Self::Decomposition<'b> { + MyCow::Borrowed(r) + } +} + +/// Helper trait for functions to work with either owned values or references to either the +/// “principal type” `X` or types some present a subset of `X`. In the latter sense, this +/// generalises [`std::borrow::ToOwned`], [`std::borrow::Borrow`], and [`std::borrow::Cow`]. +/// +/// This is used, for example, by [`crate::mapping::Mapping::apply`]. +pub trait Instance<X : Space, D = <X as Space>::Decomp> : Sized where D : Decomposition<X> { + /// Decomposes self according to `decomposer`. + fn decompose<'b>(self) -> D::Decomposition<'b> + where Self : 'b, X : 'b; + + /// Returns a lightweight instance of `self`. + fn ref_instance(&self) -> D::Reference<'_>; + + /// Returns an owned instance of `X`, cloning or converting non-true instances when necessary. + fn own(self) -> X; + + // ************** automatically implemented methods below from here ************** + + /// Returns an owned instance or reference to `X`, converting non-true instances when necessary. + /// + /// Default implementation uses [`Self::own`]. Consumes the input. + fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b { + MyCow::Owned(self.own()) + } + + #[inline] + /// Evaluates `f` on a reference to self. + /// + /// Default implementation uses [`Self::cow`]. Consumes the input. + fn eval<'b, R>(self, f : impl FnOnce(&X) -> R) -> R + where X : 'b, Self : 'b + { + f(&*self.cow()) + } + + #[inline] + /// Evaluates `f` or `g` depending on whether a reference or owned value is available. + /// + /// Default implementation uses [`Self::cow`]. Consumes the input. + fn either<'b, R>( + self, + f : impl FnOnce(X) -> R, + g : impl FnOnce(&X) -> R + ) -> R + where Self : 'b + { + match self.cow() { + EitherDecomp::Owned(x) => f(x), + EitherDecomp::Borrowed(x) => g(x), + } + } +} + + +impl<X : Space + Clone> Instance<X, BasicDecomposition> for X { + #[inline] + fn decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b> + where Self : 'b, X : 'b + { + MyCow::Owned(self) + } + + #[inline] + fn own(self) -> X { + self + } + + #[inline] + fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b { + MyCow::Owned(self) + } + + #[inline] + fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> { + self + } +} + +impl<'a, X : Space + Clone> Instance<X, BasicDecomposition> for &'a X { + #[inline] + fn decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b> + where Self : 'b, X : 'b + { + MyCow::Borrowed(self) + } + + #[inline] + fn own(self) -> X { + self.clone() + } + + #[inline] + fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b { + MyCow::Borrowed(self) + } + + #[inline] + fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> { + *self + } +} + +impl<'a, X : Space + Clone> Instance<X, BasicDecomposition> for &'a mut X { + #[inline] + fn decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b> + where Self : 'b, X : 'b + { + EitherDecomp::Borrowed(self) + } + + #[inline] + fn own(self) -> X { + self.clone() + } + + #[inline] + fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b, X : Clone { + EitherDecomp::Borrowed(self) + } + + #[inline] + fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> { + *self + } +} + +impl<'a, X : Space + Clone> Instance<X, BasicDecomposition> for MyCow<'a, X> { + + #[inline] + fn decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b> + where Self : 'b, X : 'b + { + self + } + + #[inline] + fn own(self) -> X { + match self { + MyCow::Borrowed(a) => a.own(), + MyCow::Owned(b) => b.own() + } + } + + #[inline] + fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b { + match self { + MyCow::Borrowed(a) => a.cow(), + MyCow::Owned(b) => b.cow() + } + } + + #[inline] + fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> { + match self { + MyCow::Borrowed(a) => a, + MyCow::Owned(b) => &b, + } + } +} + +/// Marker type for mutable decompositions to be used with [`InstanceMut`]. +pub trait DecompositionMut<X : Space> : Sized { + type ReferenceMut<'b> : InstanceMut<X, Self> where X : 'b; +} + + +/// Helper trait for functions to work with mutable references. +pub trait InstanceMut<X : Space , D = <X as Space>::Decomp> : Sized where D : DecompositionMut<X> { + /// Returns a mutable decomposition of self. + fn ref_instance_mut(&mut self) -> D::ReferenceMut<'_>; +} + +impl<X : Space> DecompositionMut<X> for BasicDecomposition { + type ReferenceMut<'b> = &'b mut X where X : 'b; +} + +/// This impl may seem pointless, but allows throwaway mutable scratch variables +impl<'a, X : Space> InstanceMut<X, BasicDecomposition> for X { + #[inline] + fn ref_instance_mut(&mut self) + -> <BasicDecomposition as DecompositionMut<X>>::ReferenceMut<'_> + { + self + } +} + +impl<'a, X : Space> InstanceMut<X, BasicDecomposition> for &'a mut X { + #[inline] + fn ref_instance_mut(&mut self) + -> <BasicDecomposition as DecompositionMut<X>>::ReferenceMut<'_> + { + self + } +}
--- a/src/iterate.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/iterate.rs Mon Feb 03 19:22:16 2025 -0500 @@ -19,6 +19,7 @@ .. Default::default() }; let mut x = 1 as float; +# let mut iter_clone = iter.clone(); iter.iterate(|state|{ // This is our computational step x = x + x.sqrt(); @@ -26,7 +27,17 @@ // return current value when requested return x }) -}) +}); +// or alternatively (avoiding problems with moves) +# iter = iter_clone; +for state in iter.iter() { + // This is our computational step + x = x + x.sqrt(); + state.if_verbose(||{ + // return current value when requested + return x + }) +} ``` There is no colon after `state.if_verbose`, because we need to return its value. If you do something after the step, you need to store the result in a variable and return it from the main computational step afterwards. @@ -52,6 +63,8 @@ use std::marker::PhantomData; use std::time::Duration; use std::error::Error; +use std::cell::RefCell; +use std::rc::Rc; use crate::types::*; use crate::logger::*; @@ -114,7 +127,7 @@ /// /// This is the parameter obtained by the closure passed to [`AlgIterator::iterate`] or /// [`AlgIteratorFactory::iterate`]. -pub trait AlgIteratorState { +pub trait AlgIteratorState : Sized { /// Call `call_objective` if this is a verbose iteration. /// /// Verbosity depends on the [`AlgIterator`] that produced this state. @@ -122,15 +135,18 @@ /// The closure `calc_objective` should return an arbitrary value of type `V`, to be inserted /// into the log, or whatever is deemed by the [`AlgIterator`]. For usage instructions see the /// [module documentation][self]. - fn if_verbose<V, E : Error>(&self, calc_objective : impl FnMut() -> V) -> Step<V, E>; + fn if_verbose<V, E : Error>(self, calc_objective : impl FnOnce() -> V) -> Step<V, Self, E>; /// Returns the current iteration count. fn iteration(&self) -> usize; + + /// Indicates whether the iterator is quiet + fn is_quiet(&self) -> bool; } /// Result of a step of an [`AlgIterator`] #[derive(Debug, Serialize)] -pub enum Step<V, Fail : Error = std::convert::Infallible> { +pub enum Step<V, S, Fail : Error = std::convert::Infallible> { /// Iteration should be terminated Terminated, /// Iteration should be terminated due to failure @@ -138,14 +154,14 @@ /// No result this iteration (due to verbosity settings) Quiet, /// Result of this iteration (due to verbosity settings) - Result(V), + Result(V, S), } -impl<V, E : Error> Step<V, E> { +impl<V, S, E : Error> Step<V, S, E> { /// Maps the value contained within the `Step`, if any, by the closure `f`. - pub fn map<U>(self, mut f : impl FnMut(V) -> U) -> Step<U, E> { + pub fn map<U>(self, mut f : impl FnMut(V) -> U) -> Step<U, S, E> { match self { - Step::Result(v) => Step::Result(f(v)), + Step::Result(v, s) => Step::Result(f(v), s), Step::Failure(e) => Step::Failure(e), Step::Quiet => Step::Quiet, Step::Terminated => Step::Terminated, @@ -153,6 +169,12 @@ } } +impl<V, S, E : Error> Default for Step<V, S, E> { + fn default() -> Self { + Step::Quiet + } +} + /// An iterator for algorithms, produced by [`AlgIteratorFactory::prepare`]. /// /// Typically not accessed directly, but transparently produced by an [`AlgIteratorFactory`]. @@ -160,13 +182,26 @@ pub trait AlgIterator : Sized { /// The state type type State : AlgIteratorState; - /// The output type for [`Self::step`]. + /// The output type for [`Self::poststep`] and [`Self::step`]. type Output; - /// The error type for [`Self::step`] and [`Self::iterate`]. - type Err : Error; + /// The input type for [`Self::poststep`]. + type Input; - /// Advance the iterator. - fn step(&mut self) -> Step<Self::Output, Self::Err>; + /// Advance the iterator, performing `step_fn` with the state + fn step<F, E>(&mut self, step_fn : &mut F) -> Step<Self::Output, Self::State, E> + where F : FnMut(Self::State) -> Step<Self::Input, Self::State, E>, + E : Error { + self.prestep().map_or(Step::Terminated, + |state| self.poststep(step_fn(state))) + } + + /// Initial stage of advancing the iterator, before the actual step + fn prestep(&mut self) -> Option<Self::State>; + + /// Handle step result + fn poststep<E>(&mut self, result : Step<Self::Input, Self::State, E>) + -> Step<Self::Output, Self::State, E> + where E : Error; /// Return current iteration count. fn iteration(&self) -> usize { @@ -176,56 +211,29 @@ /// Return current state. fn state(&self) -> Self::State; - /// Iterate the `AlgIterator` until termination. + /// Iterate the `AlgIterator` until termination, erforming `step_fn` on each step. /// /// Returns either `()` or an error if the step closure terminated in [`Step::Failure´]. #[inline] - fn iterate(&mut self) -> Result<(), Self::Err> { + fn iterate<F, E>(&mut self, mut step_fn : F) -> Result<(), E> + where F : FnMut(Self::State) -> Step<Self::Input, Self::State, E>, + E : Error { loop { - match self.step() { + match self.step(&mut step_fn) { Step::Terminated => return Ok(()), Step::Failure(e) => return Err(e), _ => {}, } } } - - /// Converts the `AlgIterator` into a plain [`Iterator`]. - /// - /// [`Step::Quiet`] results are discarded, and [`Step::Failure`] results **panic**. - fn downcast(self) -> AlgIteratorI<Self> { - AlgIteratorI(self) - } -} - -/// Conversion of an `AlgIterator` into a plain [`Iterator`]. -/// -/// The conversion discards [`Step::Quiet`] and **panics** on [`Step::Failure`]. -pub struct AlgIteratorI<A>(A); - -impl<A> Iterator for AlgIteratorI<A> -where A : AlgIterator { - type Item = A::Output; - - fn next(&mut self) -> Option<A::Output> { - loop { - match self.0.step() { - Step::Result(v) => return Some(v), - Step::Failure(e) => panic!("{e:?}"), - Step::Terminated => return None, - Step::Quiet => continue, - } - } - } } /// A factory for producing an [`AlgIterator`]. /// /// For usage instructions see the [module documentation][self]. pub trait AlgIteratorFactory<V> : Sized { - type Iter<F, E> : AlgIterator<State = Self::State, Output = Self::Output, Err = E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter : AlgIterator<State = Self::State, Input = V, Output = Self::Output>; + /// The state type of the corresponding [`AlgIterator`]. /// A reference to this is passed to the closures passed to methods such as [`Self::iterate`]. type State : AlgIteratorState; @@ -235,12 +243,7 @@ type Output; /// Prepare an [`AlgIterator`], consuming the factory. - /// - /// The function `step_fn` should accept a `state` ([`AlgIteratorState`] parameter, and return - /// a [`Step`]. - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + fn prepare(self) -> Self::Iter; /// Iterate the the closure `step`. /// @@ -252,12 +255,12 @@ /// This method is equivalent to [`Self::prepare`] followed by [`AlgIterator::iterate`]. #[inline] fn iterate_fallible<F, E>(self, step : F) -> Result<(), E> - where F : FnMut(&Self::State) -> Step<V, E>, + where F : FnMut(Self::State) -> Step<V, Self::State, E>, E : Error { - self.prepare(step).iterate() + self.prepare().iterate(step) } - /// Iterate the the closure `step`. + /// Iterate the closure `step`. /// /// The closure should accept a `state` parameter (satisfying the trait [`AlgIteratorState`]), /// It should return the output of @@ -269,7 +272,7 @@ /// with the error type `E=`[`std::convert::Infallible`]. #[inline] fn iterate<F>(self, step : F) - where F : FnMut(&Self::State) -> Step<V> { + where F : FnMut(Self::State) -> Step<V, Self::State> { self.iterate_fallible(step).unwrap_or_default() } @@ -285,13 +288,14 @@ /// /// For usage instructions see the [module documentation][self]. #[inline] - fn iterate_data_fallible<F, D, I, E>(self, mut datasource : I, mut step : F) -> Result<(), E> - where F : FnMut(&Self::State, D) -> Step<V, E>, + fn iterate_data_fallible<F, D, I, E>(self, mut datasource : I, mut step : F) + -> Result<(), E> + where F : FnMut(Self::State, D) -> Step<V, Self::State, E>, I : Iterator<Item = D>, E : Error { - self.prepare(move |state| { + self.prepare().iterate(move |state| { datasource.next().map_or(Step::Terminated, |d| step(state, d)) - }).iterate() + }) } /// Iterate the closure `step` with data produced by `datasource`. @@ -306,7 +310,7 @@ /// For usage instructions see the [module documentation][self]. #[inline] fn iterate_data<F, D, I>(self, datasource : I, step : F) - where F : FnMut(&Self::State, D) -> Step<V>, + where F : FnMut(Self::State, D) -> Step<V, Self::State>, I : Iterator<Item = D> { self.iterate_data_fallible(datasource, step).unwrap_or_default() } @@ -395,6 +399,26 @@ /// Is the iterator quiet, i.e., on-verbose? fn is_quiet(&self) -> bool { false } + + /// Returns an an [`std::iter::Iterator`] that can be used in a `for`-loop. + fn iter(self) -> AlgIteratorIterator<Self::Iter> { + AlgIteratorIterator { + algi : Rc::new(RefCell::new(self.prepare())), + } + } + + /// Returns an an [`std::iter::Iterator`] that can be used in a `for`-loop, + /// also inputting an initial iteration status calculated by `f` if needed. + fn iter_init(self, f : impl FnOnce() -> <Self::Iter as AlgIterator>::Input) + -> AlgIteratorIterator<Self::Iter> { + let mut i = self.prepare(); + let st = i.state(); + let step : Step<<Self::Iter as AlgIterator>::Input, Self::State> = st.if_verbose(f); + i.poststep(step); + AlgIteratorIterator { + algi : Rc::new(RefCell::new(i)), + } + } } /// Options for [`BasicAlgIteratorFactory`]. @@ -427,6 +451,11 @@ /// * every 10 iterations from there on until 100 iterations, /// * every 100 iteartinos frmo there on until 1000 iterations, etc. Logarithmic(usize), + /// Same as `Logarithmic`, but $\log_b(n)$ is replaced by $min\{c, \log_b(n)\}$ where $c$ + /// is the given `cap`. For example, with `base=10` and `cap=2`, the first ten iterations + /// will be output, then every tenth iteration, and after 100 iterations, every 100th iteration, + /// without further logarithmic progression. + LogarithmicCap{ base : usize, cap : u32 }, } impl Verbose { @@ -443,6 +472,10 @@ let every = base.pow((iter as float).log(base as float).floor() as u32); iter % every == 0 } + &Verbose::LogarithmicCap{base, cap} => { + let every = base.pow(((iter as float).log(base as float).floor() as u32).min(cap)); + iter % every == 0 + } } } } @@ -467,6 +500,8 @@ verbose : bool, /// Whether results should be calculated. calc : bool, + /// Indicates whether the iteration is quiet + quiet : bool, } /// [`AlgIteratorFactory`] for [`BasicAlgIterator`] @@ -478,11 +513,10 @@ /// The simplest [`AlgIterator`], created by [`BasicAlgIteratorFactory`] #[derive(Clone,Debug)] -pub struct BasicAlgIterator<F, V, E : Error> { +pub struct BasicAlgIterator<V> { options : AlgIteratorOptions, iter : usize, - step_fn : F, - _phantoms : PhantomData<(V, E)>, + _phantoms : PhantomData<V>, } impl AlgIteratorOptions { @@ -500,18 +534,13 @@ impl<V> AlgIteratorFactory<V> for AlgIteratorOptions where V : LogRepr { type State = BasicState; - type Iter<F, E> = BasicAlgIterator<F, V, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = BasicAlgIterator<V>; type Output = V; - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { BasicAlgIterator{ options : self, iter : 0, - step_fn, _phantoms : PhantomData, } } @@ -525,18 +554,13 @@ impl<V> AlgIteratorFactory<V> for BasicAlgIteratorFactory<V> where V : LogRepr { type State = BasicState; - type Iter<F, E> = BasicAlgIterator<F, V, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = BasicAlgIterator<V>; type Output = V; - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { BasicAlgIterator { options : self.options, iter : 0, - step_fn, _phantoms : PhantomData } } @@ -547,33 +571,33 @@ } } -impl<F, V, E> AlgIterator for BasicAlgIterator<F, V, E> -where V : LogRepr, - E : Error, - F : FnMut(&BasicState) -> Step<V, E> { +impl<V> AlgIterator for BasicAlgIterator<V> +where V : LogRepr { type State = BasicState; type Output = V; - type Err = E; + type Input = V; #[inline] - fn step(&mut self) -> Step<V, E> { + fn prestep(&mut self) -> Option<Self::State> { if self.iter >= self.options.max_iter { - Step::Terminated + None } else { self.iter += 1; - let state = self.state(); - let res = (self.step_fn)(&state); - if let Step::Result(ref val) = res { - if state.verbose && !self.options.quiet { - println!("{}{}/{} {}{}", "".dimmed(), - state.iter, - self.options.max_iter, - val.logrepr(), - "".clear()); - } + Some(self.state()) + } + } + + fn poststep<E : Error>(&mut self, res : Step<V, Self::State, E>) -> Step<V, Self::State, E> { + if let Step::Result(ref val, ref state) = res { + if state.verbose && !self.options.quiet { + println!("{}{}/{} {}{}", "".dimmed(), + state.iter, + self.options.max_iter, + val.logrepr(), + "".clear()); } - res } + res } #[inline] @@ -585,19 +609,20 @@ fn state(&self) -> BasicState { let iter = self.iter; let verbose = self.options.verbose_iter.is_verbose(iter); - BasicState{ + BasicState { iter : iter, verbose : verbose, calc : verbose, + quiet : self.options.quiet } } } impl AlgIteratorState for BasicState { #[inline] - fn if_verbose<V, E : Error>(&self, mut calc_objective : impl FnMut() -> V) -> Step<V, E> { + fn if_verbose<V, E : Error>(self, calc_objective : impl FnOnce() -> V) -> Step<V, Self, E> { if self.calc { - Step::Result(calc_objective()) + Step::Result(calc_objective(), self) } else { Step::Quiet } @@ -605,7 +630,12 @@ #[inline] fn iteration(&self) -> usize { - return self.iter; + self.iter + } + + #[inline] + fn is_quiet(&self) -> bool { + self.quiet } } @@ -636,17 +666,13 @@ for StallIteratorFactory<U, BaseFactory> where BaseFactory : AlgIteratorFactory<V, Output=U>, U : SignedNum + PartialOrd { - type Iter<F, E> = StallIterator<U, BaseFactory::Iter<F, E>> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = StallIterator<U, BaseFactory::Iter>; type State = BaseFactory::State; type Output = BaseFactory::Output; - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { StallIterator { - base_iterator : self.base_options.prepare(step_fn), + base_iterator : self.base_options.prepare(), stall : self.stall, previous_value : None, } @@ -662,18 +688,24 @@ where BaseIterator : AlgIterator<Output=U>, U : SignedNum + PartialOrd { type State = BaseIterator::State; - type Output = BaseIterator::Output; - type Err = BaseIterator::Err; + type Output = U; + type Input = BaseIterator::Input; #[inline] - fn step(&mut self) -> Step<U, Self::Err> { - match self.base_iterator.step() { - Step::Result(nv) => { + fn prestep(&mut self) -> Option<Self::State> { + self.base_iterator.prestep() + } + + #[inline] + fn poststep<E>(&mut self, res : Step<Self::Input, Self::State, E>) -> Step<U, Self::State, E> + where E : Error { + match self.base_iterator.poststep(res) { + Step::Result(nv, state) => { let previous_v = self.previous_value; self.previous_value = Some(nv); match previous_v { Some(pv) if (nv - pv).abs() <= self.stall * pv.abs() => Step::Terminated, - _ => Step::Result(nv), + _ => Step::Result(nv, state), } }, val => val, @@ -711,17 +743,13 @@ for ValueIteratorFactory<U, BaseFactory> where BaseFactory : AlgIteratorFactory<V, Output=U>, U : SignedNum + PartialOrd { - type Iter<F, E> = ValueIterator<U, BaseFactory::Iter<F, E>> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = ValueIterator<U, BaseFactory::Iter>; type State = BaseFactory::State; type Output = BaseFactory::Output; - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { ValueIterator { - base_iterator : self.base_options.prepare(step_fn), + base_iterator : self.base_options.prepare(), target : self.target } } @@ -736,17 +764,22 @@ where BaseIterator : AlgIterator<Output=U>, U : SignedNum + PartialOrd { type State = BaseIterator::State; - type Output = BaseIterator::Output; - type Err = BaseIterator::Err; + type Output = U; + type Input = BaseIterator::Input; #[inline] - fn step(&mut self) -> Step<U, Self::Err> { - match self.base_iterator.step() { - Step::Result(v) => { + fn prestep(&mut self) -> Option<Self::State> { + self.base_iterator.prestep() + } + + #[inline] + fn poststep<E>(&mut self, res : Step<Self::Input, Self::State, E>) -> Step<U, Self::State, E> where E : Error{ + match self.base_iterator.poststep(res) { + Step::Result(v, state) => { if v <= self.target { Step::Terminated } else { - Step::Result(v) + Step::Result(v, state) } }, val => val, @@ -792,16 +825,12 @@ where BaseFactory : AlgIteratorFactory<V>, BaseFactory::Output : 'log { type State = BaseFactory::State; - type Iter<F, E> = LoggingIterator<'log, BaseFactory::Output, BaseFactory::Iter<F, E>> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = LoggingIterator<'log, BaseFactory::Output, BaseFactory::Iter>; type Output = (); - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { LoggingIterator { - base_iterator : self.base_options.prepare(step_fn), + base_iterator : self.base_options.prepare(), logger : self.logger, } } @@ -818,12 +847,17 @@ BaseIterator::Output : 'log { type State = BaseIterator::State; type Output = (); - type Err = BaseIterator::Err; + type Input = BaseIterator::Input; #[inline] - fn step(&mut self) -> Step<Self::Output, Self::Err> { - match self.base_iterator.step() { - Step::Result(v) => { + fn prestep(&mut self) -> Option<Self::State> { + self.base_iterator.prestep() + } + + #[inline] + fn poststep<E>(&mut self, res : Step<Self::Input, Self::State, E>) -> Step<(), Self::State, E> where E : Error { + match self.base_iterator.poststep(res) { + Step::Result(v, _) => { self.logger.log(v); Step::Quiet }, @@ -845,7 +879,7 @@ } /// This [`AlgIteratorFactory`] allows output mapping. -/// +/// /// Typically produced with [`AlgIteratorFactory::mapped`]. #[derive(Debug)] pub struct MappingIteratorFactory<G, BaseFactory> { @@ -868,16 +902,12 @@ where BaseFactory : AlgIteratorFactory<V>, G : Fn(usize, BaseFactory::Output) -> U { type State = BaseFactory::State; - type Iter<F, E> = MappingIterator<G, BaseFactory::Iter<F, E>> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = MappingIterator<G, BaseFactory::Iter>; type Output = U; - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { MappingIterator { - base_iterator : self.base_options.prepare(step_fn), + base_iterator : self.base_options.prepare(), map : self.map } } @@ -894,12 +924,17 @@ G : Fn(usize, BaseIterator::Output) -> U { type State = BaseIterator::State; type Output = U; - type Err = BaseIterator::Err; + type Input = BaseIterator::Input; #[inline] - fn step(&mut self) -> Step<Self::Output, Self::Err> { - match self.base_iterator.step() { - Step::Result(v) => Step::Result((self.map)(self.iteration(), v)), + fn prestep(&mut self) -> Option<Self::State> { + self.base_iterator.prestep() + } + + #[inline] + fn poststep<E>(&mut self, res : Step<Self::Input, Self::State, E>) -> Step<Self::Output, Self::State, E> where E : Error { + match self.base_iterator.poststep(res) { + Step::Result(v, state) => Step::Result((self.map)(self.iteration(), v), state), Step::Quiet => Step::Quiet, Step::Terminated => Step::Terminated, Step::Failure(e) => Step::Failure(e), @@ -935,7 +970,11 @@ /// Data `U` with production time attached #[derive(Copy, Clone, Debug, Serialize)] pub struct Timed<U> { + /// CPU time taken pub cpu_time : Duration, + /// Iteration number + pub iter : usize, + /// User data //#[serde(flatten)] pub data : U } @@ -951,16 +990,12 @@ for TimingIteratorFactory<BaseFactory> where BaseFactory : AlgIteratorFactory<V> { type State = BaseFactory::State; - type Iter<F, E> = TimingIterator<BaseFactory::Iter<F, E>> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error; + type Iter = TimingIterator<BaseFactory::Iter>; type Output = Timed<BaseFactory::Output>; - fn prepare<F, E>(self, step_fn : F) -> Self::Iter<F, E> - where F : FnMut(&Self::State) -> Step<V, E>, - E : Error { + fn prepare(self) -> Self::Iter { TimingIterator { - base_iterator : self.0.prepare(step_fn), + base_iterator : self.0.prepare(), start_time : ProcessTime::now() } } @@ -976,16 +1011,22 @@ where BaseIterator : AlgIterator { type State = BaseIterator::State; type Output = Timed<BaseIterator::Output>; - type Err = BaseIterator::Err; + type Input = BaseIterator::Input; + + #[inline] + fn prestep(&mut self) -> Option<Self::State> { + self.base_iterator.prestep() + } #[inline] - fn step(&mut self) -> Step<Self::Output, Self::Err> { - match self.base_iterator.step() { - Step::Result(data) => { + fn poststep<E>(&mut self, res : Step<Self::Input, Self::State, E>) -> Step<Self::Output, Self::State, E> where E : Error { + match self.base_iterator.poststep(res) { + Step::Result(data, state) => { Step::Result(Timed{ cpu_time : self.start_time.elapsed(), + iter : self.iteration(), data - }) + }, state) }, Step::Quiet => Step::Quiet, Step::Terminated => Step::Terminated, @@ -1005,6 +1046,82 @@ } // +// New for-loop interface +// + +pub struct AlgIteratorIterator<I : AlgIterator> { + algi : Rc<RefCell<I>>, +} + +pub struct AlgIteratorIteration<I : AlgIterator> { + state : I::State, + algi : Rc<RefCell<I>>, +} + +impl<I : AlgIterator> std::iter::Iterator for AlgIteratorIterator<I> { + type Item = AlgIteratorIteration<I>; + + fn next(&mut self) -> Option<Self::Item> { + let algi = self.algi.clone(); + RefCell::borrow_mut(&self.algi).prestep().map(|state| AlgIteratorIteration { + state, + algi, + }) + } +} + +/// Types of errors that may occur +#[derive(Debug,PartialEq,Eq)] +pub enum IterationError { + /// [`AlgIteratorIteration::if_verbose_check`] is not called in iteration order. + ReportingOrderingError +} + +impl<I : AlgIterator> AlgIteratorIteration<I> { + /// Call `call_objective` if this is a verbose iteration. + /// + /// Verbosity depends on the [`AlgIterator`] that produced this state. + /// + /// The closure `calc_objective` should return an arbitrary value of type `V`, to be inserted + /// into the log, or whatever is deemed by the [`AlgIterator`]. For usage instructions see the + /// [module documentation][self]. + /// + /// This function may panic if result reporting is not ordered correctly (an unlikely mistake + /// if using this facility correctly). For a version that propagates errors, see + /// [`Self::if_verbose_check`]. + pub fn if_verbose(self, calc_objective : impl FnOnce() -> I::Input) { + self.if_verbose_check(calc_objective).unwrap() + } + + /// Version of [`Self::if_verbose`] that propagates errors instead of panicking. + pub fn if_verbose_check(self, calc_objective : impl FnOnce() -> I::Input) + -> Result<(), IterationError> { + let mut algi = match RefCell::try_borrow_mut(&self.algi) { + Err(_) => return Err(IterationError::ReportingOrderingError), + Ok(algi) => algi + }; + if self.state.iteration() != algi.iteration() { + Err(IterationError::ReportingOrderingError) + } else { + let res : Step<I::Input, I::State, std::convert::Infallible> + = self.state.if_verbose(calc_objective); + algi.poststep(res); + Ok(()) + } + } + + /// Returns the current iteration count. + pub fn iteration(&self) -> usize { + self.state.iteration() + } + + /// Indicates whether the iterator is quiet + pub fn is_quiet(&self) -> bool { + self.state.is_quiet() + } +} + +// // Tests // @@ -1050,4 +1167,44 @@ .collect::<Vec<int>>()) } } + + #[test] + fn iteration_for_loop() { + let options = AlgIteratorOptions{ + max_iter : 10, + verbose_iter : Verbose::Every(3), + .. Default::default() + }; + + { + let mut start = 1 as int; + for state in options.iter() { + start = start * 2; + state.if_verbose(|| start) + } + assert_eq!(start, (2 as int).pow(10)); + } + + { + let mut start = 1 as int; + let mut log = Logger::new(); + let factory = options.instantiate() + .with_iteration_number() + .into_log(&mut log); + for state in factory.iter() { + start = start * 2; + state.if_verbose(|| start) + } + assert_eq!(start, (2 as int).pow(10)); + assert_eq!(log.data() + .iter() + .map(|LogItem{ data : v, iter : _ }| v.clone()) + .collect::<Vec<int>>(), + (1..10).map(|i| (2 as int).pow(i)) + .skip(2) + .step_by(3) + .collect::<Vec<int>>()) + } + } + }
--- a/src/lib.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/lib.rs Mon Feb 03 19:22:16 2025 -0500 @@ -7,22 +7,16 @@ #![allow(mixed_script_confusables)] #![allow(confusable_idents)] -#![feature(maybe_uninit_uninit_array,maybe_uninit_array_assume_init,maybe_uninit_slice)] -#![feature(try_trait_v2_residual,try_trait_v2)] - -#![feature(array_methods)] - -#![feature(arc_unwrap_or_clone)] - -#![feature(float_minimum_maximum)] - -#![feature(get_mut_unchecked)] - -// They don't work: -//#![feature(negative_impls)] -//#![feature(specialization)] +#![cfg_attr(feature = "nightly", + feature(maybe_uninit_uninit_array,maybe_uninit_array_assume_init,maybe_uninit_slice), + feature(float_minimum_maximum), + feature(get_mut_unchecked), + feature(cow_is_borrowed), +)] pub mod types; +pub mod instance; +pub mod collection; pub mod nanleast; pub mod error; pub mod parallelism; @@ -45,6 +39,10 @@ pub mod fe_model; pub mod bisection_tree; pub mod nalgebra_support; - +pub(crate) mod metaprogramming; +pub mod direct_product; +pub mod convex; +pub mod discrete_gradient; +pub mod operator_arithmetic; pub use types::*;
--- a/src/linops.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/linops.rs Mon Feb 03 19:22:16 2025 -0500 @@ -4,58 +4,79 @@ use numeric_literals::replace_float_literals; use std::marker::PhantomData; +use serde::Serialize; use crate::types::*; -use serde::Serialize; -pub use crate::mapping::Apply; +pub use crate::mapping::{Mapping, Space, Composition}; +use crate::direct_product::Pair; +use crate::instance::Instance; +use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm}; /// Trait for linear operators on `X`. -pub trait Linear<X> : Apply<X, Output=Self::Codomain> - + for<'a> Apply<&'a X, Output=Self::Codomain> { - type Codomain; -} +pub trait Linear<X : Space> : Mapping<X> +{ } /// Efficient in-place summation. #[replace_float_literals(F::cast_from(literal))] -pub trait AXPY<F : Num, X = Self> { +pub trait AXPY<F, X = Self> : Space + std::ops::MulAssign<F> +where + F : Num, + X : Space, +{ + type Owned : AXPY<F, X>; + /// Computes `y = βy + αx`, where `y` is `Self`. - fn axpy(&mut self, α : F, x : &X, β : F); + fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F); /// Copies `x` to `self`. - fn copy_from(&mut self, x : &X) { + fn copy_from<I : Instance<X>>(&mut self, x : I) { self.axpy(1.0, x, 0.0) } /// Computes `y = αx`, where `y` is `Self`. - fn scale_from(&mut self, α : F, x : &X) { + fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) { self.axpy(α, x, 0.0) } + + /// Return a similar zero as `self`. + fn similar_origin(&self) -> Self::Owned; + + /// Set self to zero. + fn set_zero(&mut self); } /// Efficient in-place application for [`Linear`] operators. #[replace_float_literals(F::cast_from(literal))] -pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> { +pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> { /// Computes `y = αAx + βy`, where `A` is `Self`. - fn gemv(&self, y : &mut Y, α : F, x : &X, β : F); + fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F); + #[inline] /// Computes `y = Ax`, where `A` is `Self` - fn apply_mut(&self, y : &mut Y, x : &X){ + fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ self.gemv(y, 1.0, x, 0.0) } + #[inline] /// Computes `y += Ax`, where `A` is `Self` - fn apply_add(&self, y : &mut Y, x : &X){ + fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ self.gemv(y, 1.0, x, 1.0) } } /// Bounded linear operators -pub trait BoundedLinear<X> : Linear<X> { - type FloatType : Float; +pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X> +where + F : Num, + X : Space + Norm<F, XExp>, + XExp : NormExponent, + CodExp : NormExponent +{ /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. /// This is not expected to be the norm, just any bound on it that can be - /// reasonably implemented. - fn opnorm_bound(&self) -> Self::FloatType; + /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm + /// in `X`, and `codexp` in the codomain. + fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F; } // Linear operator application into mutable target. The [`AsRef`] bound @@ -69,16 +90,16 @@ }*/ /// Trait for forming the adjoint operator of `Self`. -pub trait Adjointable<X,Yʹ> : Linear<X> { - type AdjointCodomain; +pub trait Adjointable<X, Yʹ> : Linear<X> +where + X : Space, + Yʹ : Space, +{ + type AdjointCodomain : Space; type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a; /// Form the adjoint operator of `self`. fn adjoint(&self) -> Self::Adjoint<'_>; - - /*fn adjoint_apply(&self, y : &Yʹ) -> Self::AdjointCodomain { - self.adjoint().apply(y) - }*/ } /// Trait for forming a preadjoint of an operator. @@ -88,67 +109,576 @@ /// operator. The space `Ypre` is the predual of its codomain, and should be the /// domain of the adjointed operator. `Self::Preadjoint` should be /// [`Adjointable`]`<'a,Ypre,X>`. -pub trait Preadjointable<X,Ypre> : Linear<X> { - type PreadjointCodomain; - type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a; +/// We do not make additional restrictions on `Self::Preadjoint` (in particular, it +/// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint +/// have the full space as the codomain, etc. +pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> { + type PreadjointCodomain : Space; + type Preadjoint<'a> : Linear< + Ypre, Codomain=Self::PreadjointCodomain + > where Self : 'a; - /// Form the preadjoint operator of `self`. + /// Form the adjoint operator of `self`. fn preadjoint(&self) -> Self::Preadjoint<'_>; } -/// Adjointable operators $A: X → Y$ on between reflexive spaces $X$ and $Y$. -pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {} -impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {} +/// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$. +pub trait SimplyAdjointable<X : Space> : Adjointable<X,<Self as Mapping<X>>::Codomain> {} +impl<'a,X : Space, T> SimplyAdjointable<X> for T +where T : Adjointable<X,<Self as Mapping<X>>::Codomain> {} /// The identity operator #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] pub struct IdOp<X> (PhantomData<X>); impl<X> IdOp<X> { - fn new() -> IdOp<X> { IdOp(PhantomData) } + pub fn new() -> IdOp<X> { IdOp(PhantomData) } } -impl<X> Apply<X> for IdOp<X> { - type Output = X; +impl<X : Clone + Space> Mapping<X> for IdOp<X> { + type Codomain = X; - fn apply(&self, x : X) -> X { - x + fn apply<I : Instance<X>>(&self, x : I) -> X { + x.own() } } -impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone { - type Output = X; - - fn apply(&self, x : &'a X) -> X { - x.clone() - } -} - -impl<X> Linear<X> for IdOp<X> where X : Clone { - type Codomain = X; -} - +impl<X : Clone + Space> Linear<X> for IdOp<X> +{ } #[replace_float_literals(F::cast_from(literal))] -impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone { +impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> +where + Y : AXPY<F, X>, + X : Clone + Space +{ // Computes `y = αAx + βy`, where `A` is `Self`. - fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) { + fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { y.axpy(α, x, β) } - fn apply_mut(&self, y : &mut Y, x : &X){ + fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ y.copy_from(x); } } -impl<X> BoundedLinear<X> for IdOp<X> where X : Clone { - type FloatType = float; - fn opnorm_bound(&self) -> float { 1.0 } +impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> +where + X : Space + Clone + Norm<F, E>, + F : Num, + E : NormExponent +{ + fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE } } -impl<X> Adjointable<X,X> for IdOp<X> where X : Clone { +impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> { type AdjointCodomain=X; type Adjoint<'a> = IdOp<X> where X : 'a; + fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } } +impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> { + type PreadjointCodomain=X; + type Preadjoint<'a> = IdOp<X> where X : 'a; + + fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() } +} + + +/// The zero operator +#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] +pub struct ZeroOp<'a, X, XD, Y, F> { + zero : &'a Y, // TODO: don't pass this in `new`; maybe not even store. + dual_or_predual_zero : XD, + _phantoms : PhantomData<(X, Y, F)>, +} + +// TODO: Need to make Zero in Instance. + +impl<'a, F : Num, X : Space, XD, Y : Space + Clone> ZeroOp<'a, X, XD, Y, F> { + pub fn new(zero : &'a Y, dual_or_predual_zero : XD) -> Self { + ZeroOp{ zero, dual_or_predual_zero, _phantoms : PhantomData } + } +} + +impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> { + type Codomain = Y; + + fn apply<I : Instance<X>>(&self, _x : I) -> Y { + self.zero.clone() + } +} + +impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F> +{ } + +#[replace_float_literals(F::cast_from(literal))] +impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F> +where + F : Num, + Y : AXPY<F, Y> + Clone, + X : Space +{ + // Computes `y = αAx + βy`, where `A` is `Self`. + fn gemv<I : Instance<X>>(&self, y : &mut Y, _α : F, _x : I, β : F) { + *y *= β; + } + + fn apply_mut<I : Instance<X>>(&self, y : &mut Y, _x : I){ + y.set_zero(); + } +} + +impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F> +where + X : Space + Norm<F, E1>, + Y : AXPY<F> + Clone + Norm<F, E2>, + F : Num, + E1 : NormExponent, + E2 : NormExponent, +{ + fn opnorm_bound(&self, _xexp : E1, _codexp : E2) -> F { F::ZERO } +} + +impl<'a, F : Num, X, XD, Y, Yprime : Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F> +where + X : Space, + Y : AXPY<F> + Clone + 'static, + XD : AXPY<F> + Clone + 'static, +{ + type AdjointCodomain = XD; + type Adjoint<'b> = ZeroOp<'b, Yprime, (), XD, F> where Self : 'b; + // () means not (pre)adjointable. + + fn adjoint(&self) -> Self::Adjoint<'_> { + ZeroOp::new(&self.dual_or_predual_zero, ()) + } +} + +impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F> +where + F : Num, + X : Space, + Y : AXPY<F> + Clone, + Ypre : Space, + XD : AXPY<F> + Clone + 'static, +{ + type PreadjointCodomain = XD; + type Preadjoint<'b> = ZeroOp<'b, Ypre, (), XD, F> where Self : 'b; + // () means not (pre)adjointable. + + fn preadjoint(&self) -> Self::Preadjoint<'_> { + ZeroOp::new(&self.dual_or_predual_zero, ()) + } +} + +impl<S, T, E, X> Linear<X> for Composition<S, T, E> +where + X : Space, + T : Linear<X>, + S : Linear<T::Codomain> +{ } + +impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E> +where + F : Num, + X : Space, + T : Linear<X>, + S : GEMV<F, T::Codomain, Y>, +{ + fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { + self.outer.gemv(y, α, self.inner.apply(x), β) + } + + /// Computes `y = Ax`, where `A` is `Self` + fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ + self.outer.apply_mut(y, self.inner.apply(x)) + } + + /// Computes `y += Ax`, where `A` is `Self` + fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ + self.outer.apply_add(y, self.inner.apply(x)) + } +} + +impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp> +where + F : Num, + X : Space + Norm<F, Xexp>, + Z : Space + Norm<F, Zexp>, + Xexp : NormExponent, + Yexp : NormExponent, + Zexp : NormExponent, + T : BoundedLinear<X, Xexp, Zexp, F, Codomain=Z>, + S : BoundedLinear<Z, Zexp, Yexp, F>, +{ + fn opnorm_bound(&self, xexp : Xexp, yexp : Yexp) -> F { + let zexp = self.intermediate_norm_exponent; + self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp) + } +} + +/// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$. +pub struct RowOp<S, T>(pub S, pub T); + +use std::ops::Add; + +impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T> +where + A : Space, + B : Space, + S : Mapping<A>, + T : Mapping<B>, + S::Codomain : Add<T::Codomain>, + <S::Codomain as Add<T::Codomain>>::Output : Space, + +{ + type Codomain = <S::Codomain as Add<T::Codomain>>::Output; + + fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { + let Pair(a, b) = x.decompose(); + self.0.apply(a) + self.1.apply(b) + } +} + +impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T> +where + A : Space, + B : Space, + S : Linear<A>, + T : Linear<B>, + S::Codomain : Add<T::Codomain>, + <S::Codomain as Add<T::Codomain>>::Output : Space, +{ } + + +impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T> +where + U : Space, + V : Space, + S : GEMV<F, U, Y>, + T : GEMV<F, V, Y>, + F : Num, + Self : Linear<Pair<U, V>, Codomain=Y> +{ + fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) { + let Pair(u, v) = x.decompose(); + self.0.gemv(y, α, u, β); + self.1.gemv(y, α, v, F::ONE); + } + + fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) { + let Pair(u, v) = x.decompose(); + self.0.apply_mut(y, u); + self.1.apply_add(y, v); + } + + /// Computes `y += Ax`, where `A` is `Self` + fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) { + let Pair(u, v) = x.decompose(); + self.0.apply_add(y, u); + self.1.apply_add(y, v); + } +} + +/// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$. +pub struct ColOp<S, T>(pub S, pub T); + +impl<A, S, T> Mapping<A> for ColOp<S, T> +where + A : Space, + S : Mapping<A>, + T : Mapping<A>, +{ + type Codomain = Pair<S::Codomain, T::Codomain>; + + fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain { + Pair(self.0.apply(a.ref_instance()), self.1.apply(a)) + } +} + +impl<A, S, T> Linear<A> for ColOp<S, T> +where + A : Space, + S : Mapping<A>, + T : Mapping<A>, +{ } + +impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T> +where + X : Space, + S : GEMV<F, X, A>, + T : GEMV<F, X, B>, + F : Num, + Self : Linear<X, Codomain=Pair<A, B>> +{ + fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { + self.0.gemv(&mut y.0, α, x.ref_instance(), β); + self.1.gemv(&mut y.1, α, x, β); + } + + fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){ + self.0.apply_mut(&mut y.0, x.ref_instance()); + self.1.apply_mut(&mut y.1, x); + } + + /// Computes `y += Ax`, where `A` is `Self` + fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){ + self.0.apply_add(&mut y.0, x.ref_instance()); + self.1.apply_add(&mut y.1, x); + } +} + + +impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T> +where + A : Space, + B : Space, + Yʹ : Space, + S : Adjointable<A, Yʹ>, + T : Adjointable<B, Yʹ>, + Self : Linear<Pair<A, B>>, + // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< + // Yʹ, + // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain> + // >, +{ + type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>; + type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; + + fn adjoint(&self) -> Self::Adjoint<'_> { + ColOp(self.0.adjoint(), self.1.adjoint()) + } +} + +impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T> +where + A : Space, + B : Space, + Yʹ : Space, + S : Preadjointable<A, Yʹ>, + T : Preadjointable<B, Yʹ>, + Self : Linear<Pair<A, B>>, + for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< + Yʹ, Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>, + >, +{ + type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>; + type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; + + fn preadjoint(&self) -> Self::Preadjoint<'_> { + ColOp(self.0.preadjoint(), self.1.preadjoint()) + } +} + + +impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> +where + A : Space, + Xʹ : Space, + Yʹ : Space, + R : Space + ClosedAdd, + S : Adjointable<A, Xʹ, AdjointCodomain = R>, + T : Adjointable<A, Yʹ, AdjointCodomain = R>, + Self : Linear<A>, + // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< + // Pair<Xʹ,Yʹ>, + // Codomain=R, + // >, +{ + type AdjointCodomain = R; + type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; + + fn adjoint(&self) -> Self::Adjoint<'_> { + RowOp(self.0.adjoint(), self.1.adjoint()) + } +} + +impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> +where + A : Space, + Xʹ : Space, + Yʹ : Space, + R : Space + ClosedAdd, + S : Preadjointable<A, Xʹ, PreadjointCodomain = R>, + T : Preadjointable<A, Yʹ, PreadjointCodomain = R>, + Self : Linear<A>, + for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< + Pair<Xʹ,Yʹ>, Codomain = R, + >, +{ + type PreadjointCodomain = R; + type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; + + fn preadjoint(&self) -> Self::Preadjoint<'_> { + RowOp(self.0.preadjoint(), self.1.preadjoint()) + } +} + +/// Diagonal operator +pub struct DiagOp<S, T>(pub S, pub T); + +impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T> +where + A : Space, + B : Space, + S : Mapping<A>, + T : Mapping<B>, +{ + type Codomain = Pair<S::Codomain, T::Codomain>; + + fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { + let Pair(a, b) = x.decompose(); + Pair(self.0.apply(a), self.1.apply(b)) + } +} + +impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T> +where + A : Space, + B : Space, + S : Linear<A>, + T : Linear<B>, +{ } + +impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T> +where + A : Space, + B : Space, + U : Space, + V : Space, + S : GEMV<F, U, A>, + T : GEMV<F, V, B>, + F : Num, + Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>, +{ + fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { + let Pair(u, v) = x.decompose(); + self.0.gemv(&mut y.0, α, u, β); + self.1.gemv(&mut y.1, α, v, β); + } + + fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){ + let Pair(u, v) = x.decompose(); + self.0.apply_mut(&mut y.0, u); + self.1.apply_mut(&mut y.1, v); + } + + /// Computes `y += Ax`, where `A` is `Self` + fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){ + let Pair(u, v) = x.decompose(); + self.0.apply_add(&mut y.0, u); + self.1.apply_add(&mut y.1, v); + } +} + +impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> +where + A : Space, + B : Space, + Xʹ: Space, + Yʹ : Space, + R : Space, + S : Adjointable<A, Xʹ>, + T : Adjointable<B, Yʹ>, + Self : Linear<Pair<A, B>>, + for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< + Pair<Xʹ,Yʹ>, Codomain=R, + >, +{ + type AdjointCodomain = R; + type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; + + fn adjoint(&self) -> Self::Adjoint<'_> { + DiagOp(self.0.adjoint(), self.1.adjoint()) + } +} + +impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> +where + A : Space, + B : Space, + Xʹ: Space, + Yʹ : Space, + R : Space, + S : Preadjointable<A, Xʹ>, + T : Preadjointable<B, Yʹ>, + Self : Linear<Pair<A, B>>, + for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< + Pair<Xʹ,Yʹ>, Codomain=R, + >, +{ + type PreadjointCodomain = R; + type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; + + fn preadjoint(&self) -> Self::Preadjoint<'_> { + DiagOp(self.0.preadjoint(), self.1.preadjoint()) + } +} + +/// Block operator +pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>; + + +macro_rules! pairnorm { + ($expj:ty) => { + impl<F, A, B, S, T, ExpA, ExpB, ExpR> + BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> + for RowOp<S, T> + where + F : Float, + A : Space + Norm<F, ExpA>, + B : Space + Norm<F, ExpB>, + S : BoundedLinear<A, ExpA, ExpR, F>, + T : BoundedLinear<B, ExpB, ExpR, F>, + S::Codomain : Add<T::Codomain>, + <S::Codomain as Add<T::Codomain>>::Output : Space, + ExpA : NormExponent, + ExpB : NormExponent, + ExpR : NormExponent, + { + fn opnorm_bound( + &self, + PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>, + expr : ExpR + ) -> F { + // An application of the triangle inequality bounds the norm by the maximum + // of the individual norms. A simple observation shows this to be exact. + let na = self.0.opnorm_bound(expa, expr); + let nb = self.1.opnorm_bound(expb, expr); + na.max(nb) + } + } + + impl<F, A, S, T, ExpA, ExpS, ExpT> + BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> + for ColOp<S, T> + where + F : Float, + A : Space + Norm<F, ExpA>, + S : BoundedLinear<A, ExpA, ExpS, F>, + T : BoundedLinear<A, ExpA, ExpT, F>, + ExpA : NormExponent, + ExpS : NormExponent, + ExpT : NormExponent, + { + fn opnorm_bound( + &self, + expa : ExpA, + PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj> + ) -> F { + // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, + // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} + let ns = self.0.opnorm_bound(expa, exps); + let nt = self.1.opnorm_bound(expa, expt); + ns.max(nt) + } + } + } +} + +pairnorm!(L1); +pairnorm!(L2); +pairnorm!(Linfinity); +
--- a/src/linsolve.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/linsolve.rs Mon Feb 03 19:22:16 2025 -0500 @@ -3,6 +3,7 @@ */ use crate::types::Float; +#[cfg(feature = "nightly")] use std::mem::MaybeUninit; /// Gaussian elimination for $AX=B$, where $A$ and $B$ are both stored in `ab`, @@ -45,27 +46,44 @@ // Solve UAX=UB for X where UA with U presenting the transformations above an // upper triangular matrix. + // + // If the "nightly" feature is enabled, we will use an uninitialised array for a + // little bit of efficiency. // This use of MaybeUninit assumes F : Copy. Otherwise undefined behaviour may occur. - let mut x : [[MaybeUninit<F>; K]; M] = core::array::from_fn(|_| MaybeUninit::uninit_array::<K>() ); - //unsafe { std::mem::MaybeUninit::uninit().assume_init() }; - for i in (0..M).rev() { - for 𝓁 in 0..K { - let mut tmp = ab[i][M+𝓁]; - for j in (i+1)..M { - tmp -= ab[i][j] * unsafe { *(x[j][𝓁].assume_init_ref()) }; + #[cfg(feature = "nightly")] + { + let mut x : [[MaybeUninit<F>; K]; M] = core::array::from_fn(|_| MaybeUninit::uninit_array::<K>() ); + //unsafe { std::mem::MaybeUninit::uninit().assume_init() }; + for i in (0..M).rev() { + for 𝓁 in 0..K { + let mut tmp = ab[i][M+𝓁]; + for j in (i+1)..M { + tmp -= ab[i][j] * unsafe { *(x[j][𝓁].assume_init_ref()) }; + } + tmp /= ab[i][i]; + x[i][𝓁].write(tmp); } - tmp /= ab[i][i]; - x[i][𝓁].write(tmp); + } + unsafe { + //core::intrinsics::assert_inhabited::<[[F; K]; M]>(); + (&x as *const _ as *const [[F; K]; M]).read() } } - //unsafe { MaybeUninit::array_assume_init(x) }; - let xinit = unsafe { - //core::intrinsics::assert_inhabited::<[[F; K]; M]>(); - (&x as *const _ as *const [[F; K]; M]).read() - }; - - //std::mem::forget(x); - xinit + #[cfg(not(feature = "nightly"))] + { + let mut x : [[F; K]; M] = [[F::ZERO; K]; M]; + for i in (0..M).rev() { + for 𝓁 in 0..K { + let mut tmp = ab[i][M+𝓁]; + for j in (i+1)..M { + tmp -= ab[i][j] * x[j][𝓁]; + } + tmp /= ab[i][i]; + x[i][𝓁] = tmp; + } + } + x + } } /// Gaussian elimination for $Ax=b$, where $A$ and $b$ are both stored in `ab`,
--- a/src/loc.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/loc.rs Mon Feb 03 19:22:16 2025 -0500 @@ -5,13 +5,17 @@ use std::ops::{Add,Sub,AddAssign,SubAssign,Mul,Div,MulAssign,DivAssign,Neg,Index,IndexMut}; use std::slice::{Iter,IterMut}; +use std::fmt::{Display, Formatter}; use crate::types::{Float,Num,SignedNum}; use crate::maputil::{FixedLength,FixedLengthMut,map1,map2,map1_mut,map2_mut}; use crate::euclidean::*; use crate::norms::*; -use crate::linops::AXPY; +use crate::linops::{AXPY, Mapping, Linear}; +use crate::instance::{Instance, BasicDecomposition}; +use crate::mapping::Space; use serde::ser::{Serialize, Serializer, SerializeSeq}; + /// A container type for (short) `N`-dimensional vectors of element type `F`. /// /// Supports basic operations of an [`Euclidean`] space, several [`Norm`]s, and @@ -22,6 +26,19 @@ pub [F; N] ); +impl<F : Display, const N : usize> Display for Loc<F, N>{ + // Required method + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + write!(f, "[")?; + let mut comma = ""; + for e in self.iter() { + write!(f, "{comma}{e}")?; + comma = ", "; + } + write!(f, "]") + } +} + // Need to manually implement as [F; N] serialisation is provided only for some N. impl<F, const N : usize> Serialize for Loc<F, N> where @@ -135,6 +152,14 @@ } } +impl<F> Loc<F, 1> { + #[inline] + pub fn flatten1d(self) -> F { + let Loc([v]) = self; + v + } +} + impl<F, const N : usize> From<Loc<F, N>> for [F; N] { #[inline] fn from(other : Loc<F, N>) -> [F; N] { @@ -237,6 +262,20 @@ make_binop!(Add, add, AddAssign, add_assign); make_binop!(Sub, sub, SubAssign, sub_assign); +impl<F : Float, const N : usize> std::iter::Sum for Loc<F, N> { + fn sum<I: Iterator<Item = Loc<F, N>>>(mut iter: I) -> Self { + match iter.next() { + None => Self::ORIGIN, + Some(mut v) => { + for w in iter { + v += w + } + v + } + } + } +} + macro_rules! make_scalarop_rhs { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl<F : Num, const N : usize> $trait<F> for Loc<F, N> { @@ -390,24 +429,17 @@ domination!(Linfinity, L2); domination!(L2, L1); -impl<F : Num,const N : usize> Dot<Loc<F, N>,F> for Loc<F, N> { +impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> { + type Output = Self; + /// This implementation is not stabilised as it's meant to be used for very small vectors. /// Use [`nalgebra`] for larger vectors. #[inline] - fn dot(&self, other : &Loc<F, N>) -> F { + fn dot<I : Instance<Self>>(&self, other : I) -> F { self.0.iter() - .zip(other.0.iter()) + .zip(other.ref_instance().0.iter()) .fold(F::ZERO, |m, (&v, &w)| m + v * w) } -} - -impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> { - type Output = Self; - - #[inline] - fn similar_origin(&self) -> Self { - Self::ORIGIN - } /// This implementation is not stabilised as it's meant to be used for very small vectors. /// Use [`nalgebra`] for larger vectors. @@ -416,9 +448,9 @@ self.iter().fold(F::ZERO, |m, &v| m + v * v) } - fn dist2_squared(&self, other : &Self) -> F { + fn dist2_squared<I : Instance<Self>>(&self, other : I) -> F { self.iter() - .zip(other.iter()) + .zip(other.ref_instance().iter()) .fold(F::ZERO, |m, (&v, &w)| { let d = v - w; m + d * d }) } @@ -433,10 +465,11 @@ } #[inline] - fn dist2(&self, other : &Self) -> F { + fn dist2<I : Instance<Self>>(&self, other : I) -> F { // Optimisation for N==1 that avoids squaring and square rooting. + let otherr = other.ref_instance(); if N==1 { - unsafe { *self.0.get_unchecked(0) - *other.0.get_unchecked(0) }.abs() + unsafe { *self.0.get_unchecked(0) - *otherr.0.get_unchecked(0) }.abs() } else { self.dist2_squared(other).sqrt() } @@ -444,9 +477,94 @@ } impl<F : Num, const N : usize> Loc<F, N> { + /// Origin point pub const ORIGIN : Self = Loc([F::ZERO; N]); } +impl<F : Num + std::ops::Neg<Output=F>, const N : usize> Loc<F, N> { + /// Reflects along the given coordinate + pub fn reflect(mut self, i : usize) -> Self { + self[i] = -self[i]; + self + } + + /// Reflects along the given coordinate sequences + pub fn reflect_many<I : IntoIterator<Item=usize>>(mut self, idxs : I) -> Self { + for i in idxs { + self[i]=-self[i] + } + self + } +} + +impl<F : std::ops::Neg<Output=F>> Loc<F, 2> { + /// Reflect x coordinate + pub fn reflect_x(self) -> Self { + let Loc([x, y]) = self; + [-x, y].into() + } + + /// Reflect y coordinate + pub fn reflect_y(self) -> Self { + let Loc([x, y]) = self; + [x, -y].into() + } +} + +impl<F : Float> Loc<F, 2> { + /// Rotate by angle φ + pub fn rotate(self, φ : F) -> Self { + let Loc([x, y]) = self; + let sin_φ = φ.sin(); + let cos_φ = φ.cos(); + [cos_φ * x - sin_φ * y, + sin_φ * x + cos_φ * y].into() + } +} + +impl<F : std::ops::Neg<Output=F>> Loc<F, 3> { + /// Reflect x coordinate + pub fn reflect_x(self) -> Self { + let Loc([x, y, z]) = self; + [-x, y, z].into() + } + + /// Reflect y coordinate + pub fn reflect_y(self) -> Self { + let Loc([x, y, z]) = self; + [x, -y, z].into() + } + + /// Reflect y coordinate + pub fn reflect_z(self) -> Self { + let Loc([x, y, z]) = self; + [x, y, -z].into() + } +} + +impl<F : Float> Loc<F, 3> { + /// Rotate by angles (yaw, pitch, roll) + pub fn rotate(self, Loc([φ, ψ, θ]) : Self) -> Self { + let Loc([mut x, mut y, mut z]) = self; + let sin_φ = φ.sin(); + let cos_φ = φ.cos(); + [x, y, z] = [cos_φ * x - sin_φ *y, + sin_φ * x + cos_φ * y, + z]; + let sin_ψ = ψ.sin(); + let cos_ψ = ψ.cos(); + [x, y, z] = [cos_ψ * x + sin_ψ * z, + y, + -sin_ψ * x + cos_ψ * z]; + let sin_θ = θ.sin(); + let cos_θ = θ.cos(); + [x, y, z] = [x, + cos_θ * y - sin_θ * z, + sin_θ * y + cos_θ * z]; + [x, y, z].into() + } +} + impl<F : Float,const N : usize> StaticEuclidean<F> for Loc<F, N> { #[inline] fn origin() -> Self { @@ -454,6 +572,31 @@ } } + +/// The default norm for `Loc` is [`L2`]. +impl<F : Float, const N : usize> Normed<F> for Loc<F, N> { + type NormExp = L2; + + #[inline] + fn norm_exponent(&self) -> Self::NormExp { + L2 + } + + // #[inline] + // fn similar_origin(&self) -> Self { + // [F::ZERO; N].into() + // } + + #[inline] + fn is_zero(&self) -> bool { + self.norm2_squared() == F::ZERO + } +} + +impl<F : Float, const N : usize> HasDual<F> for Loc<F, N> { + type DualSpace = Self; +} + impl<F : Float, const N : usize> Norm<F, L2> for Loc<F, N> { #[inline] fn norm(&self, _ : L2) -> F { self.norm2() } @@ -461,9 +604,20 @@ impl<F : Float, const N : usize> Dist<F, L2> for Loc<F, N> { #[inline] - fn dist(&self, other : &Self, _ : L2) -> F { self.dist2(other) } + fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> F { self.dist2(other) } } +/* Implemented automatically as Euclidean. +impl<F : Float, const N : usize> Projection<F, L2> for Loc<F, N> { + #[inline] + fn proj_ball_mut(&mut self, ρ : F, nrm : L2) { + let n = self.norm(nrm); + if n > ρ { + *self *= ρ/n; + } + } +}*/ + impl<F : Float, const N : usize> Norm<F, L1> for Loc<F, N> { /// This implementation is not stabilised as it's meant to be used for very small vectors. /// Use [`nalgebra`] for larger vectors. @@ -475,9 +629,9 @@ impl<F : Float, const N : usize> Dist<F, L1> for Loc<F, N> { #[inline] - fn dist(&self, other : &Self, _ : L1) -> F { + fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> F { self.iter() - .zip(other.iter()) + .zip(other.ref_instance().iter()) .fold(F::ZERO, |m, (&v, &w)| m + (v-w).abs() ) } } @@ -500,9 +654,9 @@ impl<F : Float, const N : usize> Dist<F, Linfinity> for Loc<F, N> { #[inline] - fn dist(&self, other : &Self, _ : Linfinity) -> F { + fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> F { self.iter() - .zip(other.iter()) + .zip(other.ref_instance().iter()) .fold(F::ZERO, |m, (&v, &w)| m.max((v-w).abs())) } } @@ -536,19 +690,48 @@ } } -impl<F : Num, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> { + +impl<F : Num, const N : usize> Space for Loc<F, N> { + type Decomp = BasicDecomposition; +} + +impl<F : Float, const N : usize> Mapping<Loc<F, N>> for Loc<F, N> { + type Codomain = F; + + fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain { + x.eval(|x̃| self.dot(x̃)) + } +} + +impl<F : Float, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { } + +impl<F : Float, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> { + type Owned = Self; #[inline] - fn axpy(&mut self, α : F, x : &Loc<F, N>, β : F) { - if β == F::ZERO { - map2_mut(self, x, |yi, xi| { *yi = α * (*xi) }) - } else { - map2_mut(self, x, |yi, xi| { *yi = β * (*yi) + α * (*xi) }) - } + fn axpy<I : Instance<Loc<F, N>>>(&mut self, α : F, x : I, β : F) { + x.eval(|x̃| { + if β == F::ZERO { + map2_mut(self, x̃, |yi, xi| { *yi = α * (*xi) }) + } else { + map2_mut(self, x̃, |yi, xi| { *yi = β * (*yi) + α * (*xi) }) + } + }) } #[inline] - fn copy_from(&mut self, x : &Loc<F, N>) { - map2_mut(self, x, |yi, xi| *yi = *xi ) + fn copy_from<I : Instance<Loc<F, N>>>(&mut self, x : I) { + x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi )) + } + + #[inline] + fn similar_origin(&self) -> Self::Owned { + Self::ORIGIN + } + + #[inline] + fn set_zero(&mut self) { + *self = Self::ORIGIN; } } +
--- a/src/logger.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/logger.rs Mon Feb 03 19:22:16 2025 -0500 @@ -28,6 +28,11 @@ pub fn data(&self) -> &Vec<V> { &self.data } + + /// Map the log with `g`. + pub fn map<W>(self, g : impl FnMut(V) -> W) -> Logger<W> { + Logger { data : self.data.into_iter().map(g).collect() } + } } impl<'a, V : Serialize + 'a> TableDump<'a> for Logger<V> {
--- a/src/mapping.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/mapping.rs Mon Feb 03 19:22:16 2025 -0500 @@ -3,115 +3,248 @@ */ use std::marker::PhantomData; -use crate::types::{Float}; -use serde::Serialize; +use std::borrow::Cow; +use crate::types::{Num, Float, ClosedMul}; use crate::loc::Loc; +pub use crate::instance::{Instance, Decomposition, BasicDecomposition, Space}; +use crate::norms::{Norm, NormExponent}; +use crate::operator_arithmetic::{Weighted, Constant}; -/// Trait for application of `Self` as a mathematical function or operator on `X`. -pub trait Apply<X> { - type Output; +/// A mapping from `Domain` to `Self::Codomain`. +pub trait Mapping<Domain : Space> { + type Codomain : Space; /// Compute the value of `self` at `x`. - fn apply(&self, x : X) -> Self::Output; -} + fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain; -/// This blanket implementation is a workaround helper to Rust trait system limitations. -/// -/// It is introduced because the trait system does not allow blanket implementations of both -/// [`Apply<X>`] and [`Apply<&'a X>`]. With this, the latter is implemented automatically for -/// the reference, which can be sufficient to apply the operation in another blanket implementation. -impl<'a, T, X> Apply<X> for &'a T where T : Apply<X> { - type Output = <T as Apply<X>>::Output; + #[inline] + /// Form the composition `self ∘ other` + fn compose<X : Space, T : Mapping<X, Codomain=Domain>>(self, other : T) + -> Composition<Self, T> + where + Self : Sized + { + Composition{ outer : self, inner : other, intermediate_norm_exponent : () } + } + #[inline] - fn apply(&self, x : X) -> Self::Output { - (*self).apply(x) + /// Form the composition `self ∘ other`, assigning a norm to the inermediate space + fn compose_with_norm<F, X, T, E>( + self, other : T, norm : E + ) -> Composition<Self, T, E> + where + Self : Sized, + X : Space, + T : Mapping<X, Codomain=Domain>, + E : NormExponent, + Domain : Norm<F, E>, + F : Num + { + Composition{ outer : self, inner : other, intermediate_norm_exponent : norm } + } + + /// Multiply `self` by the scalar `a`. + #[inline] + fn weigh<C>(self, a : C) -> Weighted<Self, C> + where + Self : Sized, + C : Constant, + Self::Codomain : ClosedMul<C::Type>, + { + Weighted { weight : a, base_fn : self } } } -/// A mapping from `Domain` to `Codomain`. -/// -/// This is automatically implemented when the relevant [`Apply`] are implemented. -pub trait Mapping<Domain> : Apply<Domain, Output=Self::Codomain> - + for<'a> Apply<&'a Domain, Output=Self::Codomain> { - type Codomain; -} - -impl<Domain, Codomain, T> Mapping<Domain> for T -where T : Apply<Domain, Output=Codomain> + for<'a> Apply<&'a Domain, Output=Codomain> { - type Codomain = Codomain; -} - - -/// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to `F` a [`Float`]. -pub trait RealMapping<F : Float, const N : usize> : Mapping<Loc<F, N>, Codomain = F> {} +/// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<F, N>`] to `F`. +pub trait RealMapping<F : Float, const N : usize> +: Mapping<Loc<F, N>, Codomain = F> {} impl<F : Float, T, const N : usize> RealMapping<F, N> for T where T : Mapping<Loc<F, N>, Codomain = F> {} - -/// Trait for calculation the differential of `Self` as a mathematical function on `X`. -pub trait Differentiate<X> { - type Output; +/// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to [`Loc<F, M>`]. +pub trait RealVectorField<F : Float, const N : usize, const M : usize> +: Mapping<Loc<F, N>, Codomain = Loc<F, M>> {} - /// Compute the differential of `self` at `x`. - fn differential(&self, x : X) -> Self::Output; -} - +impl<F : Float, T, const N : usize, const M : usize> RealVectorField<F, N, M> for T +where T : Mapping<Loc<F, N>, Codomain = Loc<F, M>> {} /// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials /// `Differential`. /// -/// This is automatically implemented when the relevant [`Differentiate`] are implemented. -pub trait DifferentiableMapping<Domain> -: Mapping<Domain> - + Differentiate<Domain, Output=Self::Differential> - + for<'a> Differentiate<&'a Domain, Output=Self::Differential>{ - type Differential; +/// This is automatically implemented when [`DifferentiableImpl`] is. +pub trait DifferentiableMapping<Domain : Space> : Mapping<Domain> { + type DerivativeDomain : Space; + type Differential<'b> : Mapping<Domain, Codomain=Self::DerivativeDomain> where Self : 'b; + + /// Calculate differential at `x` + fn differential<I : Instance<Domain>>(&self, x : I) -> Self::DerivativeDomain; + + /// Form the differential mapping of `self`. + fn diff(self) -> Self::Differential<'static>; + + /// Form the differential mapping of `self`. + fn diff_ref(&self) -> Self::Differential<'_>; } +/// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from +/// [`Loc<F, N>`] to `F`. +pub trait DifferentiableRealMapping<F : Float, const N : usize> +: DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain = Loc<F, N>> {} -impl<Domain, Differential, T> DifferentiableMapping<Domain> for T -where T : Mapping<Domain> - + Differentiate<Domain, Output=Differential> - + for<'a> Differentiate<&'a Domain, Output=Differential> { - type Differential = Differential; +impl<F : Float, T, const N : usize> DifferentiableRealMapping<F, N> for T +where T : DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain = Loc<F, N>> {} + +/// Helper trait for implementing [`DifferentiableMapping`] +pub trait DifferentiableImpl<X : Space> : Sized { + type Derivative : Space; + + /// Compute the differential of `self` at `x`, consuming the input. + fn differential_impl<I : Instance<X>>(&self, x : I) -> Self::Derivative; } -/// A sum of [`Mapping`]s. -#[derive(Serialize, Debug, Clone)] -pub struct Sum<Domain, M : Mapping<Domain>> { - components : Vec<M>, - _domain : PhantomData<Domain>, -} +impl<T, Domain> DifferentiableMapping<Domain> for T +where + Domain : Space, + T : Clone + Mapping<Domain> + DifferentiableImpl<Domain> +{ + type DerivativeDomain = T::Derivative; + type Differential<'b> = Differential<'b, Domain, Self> where Self : 'b; + + #[inline] + fn differential<I : Instance<Domain>>(&self, x : I) -> Self::DerivativeDomain { + self.differential_impl(x) + } -impl<Domain, M : Mapping<Domain>> Sum<Domain, M> { - /// Construct from an iterator. - pub fn new<I : Iterator<Item = M>>(iter : I) -> Self { - Sum { components : iter.collect(), _domain : PhantomData } + fn diff(self) -> Differential<'static, Domain, Self> { + Differential{ g : Cow::Owned(self), _space : PhantomData } + } + + fn diff_ref(&self) -> Differential<'_, Domain, Self> { + Differential{ g : Cow::Borrowed(self), _space : PhantomData } } } -impl<Domain : Copy, M> Apply<Domain> for Sum<Domain, M> -where M : Mapping<Domain>, - M::Codomain : std::iter::Sum { - type Output = M::Codomain; +/// Container for the differential [`Mapping`] of a [`DifferentiableMapping`]. +pub struct Differential<'a, X, G : Clone> { + g : Cow<'a, G>, + _space : PhantomData<X> +} + +impl<'a, X, G : Clone> Differential<'a, X, G> { + pub fn base_fn(&self) -> &G { + &self.g + } +} + +impl<'a, X, G> Mapping<X> for Differential<'a, X, G> +where + X : Space, + G : Clone + DifferentiableMapping<X> +{ + type Codomain = G::DerivativeDomain; + + #[inline] + fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain { + (*self.g).differential(x) + } +} - fn apply(&self, x : Domain) -> Self::Output { - self.components.iter().map(|c| c.apply(x)).sum() +/// Container for flattening [`Loc`]`<F, 1>` codomain of a [`Mapping`] to `F`. +pub struct FlattenedCodomain<X, F, G> { + g : G, + _phantoms : PhantomData<(X, F)> +} + +impl<F : Space, X, G> Mapping<X> for FlattenedCodomain<X, F, G> +where + X : Space, + G: Mapping<X, Codomain=Loc<F, 1>> +{ + type Codomain = F; + + #[inline] + fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain { + self.g.apply(x).flatten1d() + } +} + +/// An auto-trait for constructing a [`FlattenCodomain`] structure for +/// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`. +pub trait FlattenCodomain<X : Space, F> : Mapping<X, Codomain=Loc<F, 1>> + Sized { + /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`. + fn flatten_codomain(self) -> FlattenedCodomain<X, F, Self> { + FlattenedCodomain{ g : self, _phantoms : PhantomData } } } -impl<Domain, M> Differentiate<Domain> for Sum<Domain, M> -where M : DifferentiableMapping<Domain>, - M :: Codomain : std::iter::Sum, - M :: Differential : std::iter::Sum, - Domain : Copy { +impl<X : Space, F, G : Sized + Mapping<X, Codomain=Loc<F, 1>>> FlattenCodomain<X, F> for G {} + +/// Container for dimensional slicing [`Loc`]`<F, N>` codomain of a [`Mapping`] to `F`. +pub struct SlicedCodomain<'a, X, F, G : Clone, const N : usize> { + g : Cow<'a, G>, + slice : usize, + _phantoms : PhantomData<(X, F)> +} - type Output = M::Differential; +impl<'a, X, F, G, const N : usize> Mapping<X> for SlicedCodomain<'a, X, F, G, N> +where + X : Space, + F : Copy + Space, + G : Mapping<X, Codomain=Loc<F, N>> + Clone, +{ + type Codomain = F; - fn differential(&self, x : Domain) -> Self::Output { - self.components.iter().map(|c| c.differential(x)).sum() + #[inline] + fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain { + let tmp : [F; N] = (*self.g).apply(x).into(); + // Safety: `slice_codomain` below checks the range. + unsafe { *tmp.get_unchecked(self.slice) } } } + +/// An auto-trait for constructing a [`FlattenCodomain`] structure for +/// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`. +pub trait SliceCodomain<X : Space, F : Copy, const N : usize> + : Mapping<X, Codomain=Loc<F, N>> + Clone + Sized +{ + /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`. + fn slice_codomain(self, slice : usize) -> SlicedCodomain<'static, X, F, Self, N> { + assert!(slice < N); + SlicedCodomain{ g : Cow::Owned(self), slice, _phantoms : PhantomData } + } + + /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`. + fn slice_codomain_ref(&self, slice : usize) -> SlicedCodomain<'_, X, F, Self, N> { + assert!(slice < N); + SlicedCodomain{ g : Cow::Borrowed(self), slice, _phantoms : PhantomData } + } +} + +impl<X : Space, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>> + Clone, const N : usize> +SliceCodomain<X, F, N> +for G {} + + +/// The composition S ∘ T. `E` is for storing a `NormExponent` for the intermediate space. +pub struct Composition<S, T, E = ()> { + pub outer : S, + pub inner : T, + pub intermediate_norm_exponent : E +} + +impl<S, T, X, E> Mapping<X> for Composition<S, T, E> +where + X : Space, + T : Mapping<X>, + S : Mapping<T::Codomain> +{ + type Codomain = S::Codomain; + + #[inline] + fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain { + self.outer.apply(self.inner.apply(x)) + } +}
--- a/src/maputil.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/maputil.rs Mon Feb 03 19:22:16 2025 -0500 @@ -2,6 +2,7 @@ Utilities for mapping over various container types. */ +#[cfg(feature = "nightly")] use std::mem::MaybeUninit; use itertools::izip; @@ -322,8 +323,13 @@ /// /// If `iter.next()` panicks, all items already yielded by the iterator are /// dropped. +#[cfg(feature = "nightly")] #[inline] -pub(crate) fn collect_into_array_unchecked<T, I : Iterator<Item=T>, const N: usize>(mut iter: I) -> [T; N] +pub(crate) fn collect_into_array_unchecked< + T, + I : Iterator<Item=T>, + const N: usize +>(mut iter: I) -> [T; N] { if N == 0 { // SAFETY: An empty array is always inhabited and has no validity invariants. @@ -375,6 +381,19 @@ unreachable!("Something went wrong with iterator length") } +#[cfg(not(feature = "nightly"))] +#[inline] +pub(crate) fn collect_into_array_unchecked< + T, + I : Iterator<Item=T>, + const N: usize +>(iter: I) -> [T; N] +{ + match iter.collect::<Vec<T>>().try_into() { + Ok(a) => a, + Err(_) => panic!("collect_into_array failure: should not happen"), + } +} #[cfg(test)] mod tests {
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/metaprogramming.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,41 @@ +/*! +Metaprogramming tools +*/ + +/// Reference `x` if so indicated by the first parameter. +/// Typically to be used from another macro. +/// +/// ```ignore +/// maybe_ref!(ref, V) // ➡ &V +/// maybe_ref!(noref, V) // ➡ V +/// ``` +#[macro_export] +macro_rules! maybe_ref { + (ref, $x:expr) => { &$x }; + (noref, $x:expr) => { $x }; + (ref, $x:ty) => { &$x }; + (noref, $x:ty) => { $x }; +} + +/// Choose `a` if first argument is the literal `ref`, otherwise `b`. +#[macro_export] +macro_rules! ifref { + (noref, $a:expr, $b:expr) => { $b }; + (ref, $a:expr, $b:expr) => { $a }; +} + + +/// Annotate `x` with a lifetime if the first parameter +/// Typically to be used from another macro. +/// +/// ```ignore +/// maybe_ref!(ref, &'a V) // ➡ &'a V +/// maybe_ref!(noref, &'a V) // ➡ V +/// ``` +#[macro_export] +macro_rules! maybe_lifetime { + (ref, $x:ty) => { $x }; + (noref, &$lt:lifetime $x:ty) => { $x }; + (noref, &$x:ty) => { $x }; +} +
--- a/src/nalgebra_support.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/nalgebra_support.rs Mon Feb 03 19:22:16 2025 -0500 @@ -1,7 +1,7 @@ /*! Integration with nalgebra. -This module mainly implements [`Euclidean`], [`Norm`], [`Dot`], [`Linear`], etc. for [`nalgebra`] +This module mainly implements [`Euclidean`], [`Norm`], [`Linear`], etc. for [`nalgebra`] matrices and vectors. It also provides [`ToNalgebraRealField`] as a vomit-inducingly ugly workaround to nalgebra force-feeding its own versions of the same basic mathematical methods on `f32` and `f64` as @@ -10,10 +10,9 @@ use nalgebra::{ Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar, - ClosedMul, ClosedAdd, SimdComplexField, Vector, OVector, RealField, + ClosedAddAssign, ClosedMulAssign, SimdComplexField, Vector, OVector, RealField, LpNorm, UniformNorm }; -use nalgebra::Norm as NalgebraNorm; use nalgebra::base::constraint::{ ShapeConstraint, SameNumberOfRows, SameNumberOfColumns }; @@ -23,102 +22,127 @@ use num_traits::identities::{Zero, One}; use crate::linops::*; use crate::euclidean::*; +use crate::mapping::{Space, BasicDecomposition}; use crate::types::Float; use crate::norms::*; +use crate::instance::Instance; -impl<SM,SV,N,M,K,E> Apply<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV: Storage<E,M,K>, - N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, - DefaultAllocator : Allocator<E,N,K>, - DefaultAllocator : Allocator<E,M,K>, - DefaultAllocator : Allocator<E,N,M>, - DefaultAllocator : Allocator<E,M,N> { - type Output = OMatrix<E,N,K>; +impl<SM,N,M,E> Space for Matrix<E,N,M,SM> +where + SM: Storage<E,N,M> + Clone, + N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator : Allocator<N,M>, +{ + type Decomp = BasicDecomposition; +} + +impl<SM,SV,N,M,K,E> Mapping<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> +where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone, + N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator : Allocator<N,K>, + DefaultAllocator : Allocator<M,K>, + DefaultAllocator : Allocator<N,M>, + DefaultAllocator : Allocator<M,N> { + type Codomain = OMatrix<E,N,K>; #[inline] - fn apply(&self, x : Matrix<E,M,K,SV>) -> Self::Output { - self.mul(x) + fn apply<I : Instance<Matrix<E,M,K,SV>>>( + &self, x : I + ) -> Self::Codomain { + x.either(|owned| self.mul(owned), |refr| self.mul(refr)) } } -impl<'a, SM,SV,N,M,K,E> Apply<&'a Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV: Storage<E,M,K>, - N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, - DefaultAllocator : Allocator<E,N,K>, - DefaultAllocator : Allocator<E,M,K>, - DefaultAllocator : Allocator<E,N,M>, - DefaultAllocator : Allocator<E,M,N> { - type Output = OMatrix<E,N,K>; - - #[inline] - fn apply(&self, x : &'a Matrix<E,M,K,SV>) -> Self::Output { - self.mul(x) - } -} impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV: Storage<E,M,K>, - N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, - DefaultAllocator : Allocator<E,N,K>, - DefaultAllocator : Allocator<E,M,K>, - DefaultAllocator : Allocator<E,N,M>, - DefaultAllocator : Allocator<E,M,N> { - type Codomain = OMatrix<E,N,K>; +where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone, + N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator : Allocator<N,K>, + DefaultAllocator : Allocator<M,K>, + DefaultAllocator : Allocator<N,M>, + DefaultAllocator : Allocator<M,N> { } impl<SM,SV1,SV2,N,M,K,E> GEMV<E, Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: StorageMut<E,N,K>, - N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, - DefaultAllocator : Allocator<E,N,K>, - DefaultAllocator : Allocator<E,M,K>, - DefaultAllocator : Allocator<E,N,M>, - DefaultAllocator : Allocator<E,M,N> { +where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: StorageMut<E,N,K>, + N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float, + DefaultAllocator : Allocator<N,K>, + DefaultAllocator : Allocator<M,K>, + DefaultAllocator : Allocator<N,M>, + DefaultAllocator : Allocator<M,N> { #[inline] - fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &Matrix<E,M,K,SV1>, β : E) { - Matrix::gemm(y, α, self, x, β) + fn gemv<I : Instance<Matrix<E,M,K,SV1>>>( + &self, y : &mut Matrix<E,N,K,SV2>, α : E, x : I, β : E + ) { + x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β)) } #[inline] - fn apply_mut<'a>(&self, y : &mut Matrix<E,N,K,SV2>, x : &Matrix<E,M,K,SV1>) { - self.mul_to(x, y) + fn apply_mut<'a, I : Instance<Matrix<E,M,K,SV1>>>(&self, y : &mut Matrix<E,N,K,SV2>, x : I) { + x.eval(|x̃| self.mul_to(x̃, y)) } } impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> -where SM: StorageMut<E,M>, SV1: Storage<E,M>, - M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, - DefaultAllocator : Allocator<E,M> { +where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone, + M : Dim, E : Scalar + Zero + One + Float, + DefaultAllocator : Allocator<M> { + type Owned = OVector<E, M>; #[inline] - fn axpy(&mut self, α : E, x : &Vector<E,M,SV1>, β : E) { - Matrix::axpy(self, α, x, β) + fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) { + x.eval(|x̃| Matrix::axpy(self, α, x̃, β)) + } + + #[inline] + fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) { + y.eval(|ỹ| Matrix::copy_from(self, ỹ)) + } + + #[inline] + fn set_zero(&mut self) { + self.iter_mut().for_each(|e| *e = E::ZERO); } #[inline] - fn copy_from(&mut self, y : &Vector<E,M,SV1>) { - Matrix::copy_from(self, y) + fn similar_origin(&self) -> Self::Owned { + OVector::zeros_generic(M::from_usize(self.len()), Const) } } +/* Implemented automatically as Euclidean. +impl<SM,M,E> Projection<E, L2> for Vector<E,M,SM> +where SM: StorageMut<E,M> + Clone, + M : Dim, E : Scalar + Zero + One + Float + RealField, + DefaultAllocator : Allocator<M> { + #[inline] + fn proj_ball_mut(&mut self, ρ : E, _ : L2) { + let n = self.norm(L2); + if n > ρ { + self.iter_mut().for_each(|v| *v *= ρ/n) + } + } +}*/ + impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM> -where SM: StorageMut<E,M>, - M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float + RealField, - DefaultAllocator : Allocator<E,M> { +where SM: StorageMut<E,M> + Clone, + M : Dim, E : Scalar + Zero + One + Float + RealField, + DefaultAllocator : Allocator<M> { #[inline] fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) { self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) } } -impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>,Matrix<E,N,K,SV2>> +impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: Storage<E,N,K>, - N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField, - DefaultAllocator : Allocator<E,N,K>, - DefaultAllocator : Allocator<E,M,K>, - DefaultAllocator : Allocator<E,N,M>, - DefaultAllocator : Allocator<E,M,N> { +where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: Storage<E,N,K> + Clone, + N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField, + DefaultAllocator : Allocator<N,K>, + DefaultAllocator : Allocator<M,K>, + DefaultAllocator : Allocator<N,M>, + DefaultAllocator : Allocator<M,N> { type AdjointCodomain = OMatrix<E,M,K>; type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a; @@ -128,20 +152,6 @@ } } -impl<E,M,S,Si> Dot<Vector<E,M,Si>,E> -for Vector<E,M,S> -where M : Dim, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One, - S : Storage<E,M>, - Si : Storage<E,M>, - DefaultAllocator : Allocator<E,M> { - - #[inline] - fn dot(&self, other : &Vector<E,M,Si>) -> E { - Vector::<E,M,S>::dot(self, other) - } -} - /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`. #[inline] fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>( @@ -170,15 +180,15 @@ impl<E,M,S> Euclidean<E> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : StorageMut<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { type Output = OVector<E, M>; - + #[inline] - fn similar_origin(&self) -> OVector<E, M> { - OVector::zeros_generic(M::from_usize(self.len()), Const) + fn dot<I : Instance<Self>>(&self, other : I) -> E { + Vector::<E,M,S>::dot(self, other.ref_instance()) } #[inline] @@ -187,17 +197,17 @@ } #[inline] - fn dist2_squared(&self, other : &Self) -> E { - metric_distance_squared(self, other) + fn dist2_squared<I : Instance<Self>>(&self, other : I) -> E { + metric_distance_squared(self, other.ref_instance()) } } impl<E,M,S> StaticEuclidean<E> for Vector<E,M,S> where M : DimName, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : StorageMut<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] fn origin() -> OVector<E, M> { @@ -205,78 +215,109 @@ } } +/// The default norm for `Vector` is [`L2`]. +impl<E,M,S> Normed<E> +for Vector<E,M,S> +where M : Dim, + S : Storage<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { + + type NormExp = L2; + + #[inline] + fn norm_exponent(&self) -> Self::NormExp { + L2 + } + + #[inline] + fn is_zero(&self) -> bool { + Vector::<E,M,S>::norm_squared(self) == E::ZERO + } +} + +impl<E,M,S> HasDual<E> +for Vector<E,M,S> +where M : Dim, + S : Storage<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { + // TODO: Doesn't work with different storage formats. + type DualSpace = Self; +} + impl<E,M,S> Norm<E, L1> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : Storage<E,M>, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] fn norm(&self, _ : L1) -> E { - LpNorm(1).norm(self) + nalgebra::Norm::norm(&LpNorm(1), self) } } impl<E,M,S> Dist<E, L1> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : Storage<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] - fn dist(&self, other : &Self, _ : L1) -> E { - LpNorm(1).metric_distance(self, other) + fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> E { + nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance()) } } impl<E,M,S> Norm<E, L2> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : Storage<E,M>, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] fn norm(&self, _ : L2) -> E { - LpNorm(2).norm(self) + nalgebra::Norm::norm(&LpNorm(2), self) } } impl<E,M,S> Dist<E, L2> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : Storage<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] - fn dist(&self, other : &Self, _ : L2) -> E { - LpNorm(2).metric_distance(self, other) + fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> E { + nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance()) } } impl<E,M,S> Norm<E, Linfinity> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : Storage<E,M>, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] fn norm(&self, _ : Linfinity) -> E { - UniformNorm.norm(self) + nalgebra::Norm::norm(&UniformNorm, self) } } impl<E,M,S> Dist<E, Linfinity> for Vector<E,M,S> where M : Dim, - S : StorageMut<E,M>, - E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, - DefaultAllocator : Allocator<E,M> { + S : Storage<E,M> + Clone, + E : Float + Scalar + Zero + One + RealField, + DefaultAllocator : Allocator<M> { #[inline] - fn dist(&self, other : &Self, _ : Linfinity) -> E { - UniformNorm.metric_distance(self, other) + fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> E { + nalgebra::Norm::metric_distance(&UniformNorm, self, other.ref_instance()) } }
--- a/src/norms.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/norms.rs Mon Feb 03 19:22:16 2025 -0500 @@ -3,18 +3,31 @@ */ use serde::Serialize; +use std::marker::PhantomData; use crate::types::*; use crate::euclidean::*; +use crate::mapping::{Mapping, Space, Instance}; // // Abstract norms // +#[derive(Copy,Clone,Debug)] +/// Helper structure to convert a [`NormExponent`] into a [`Mapping`] +pub struct NormMapping<F : Float, E : NormExponent>{ + pub(crate) exponent : E, + _phantoms : PhantomData<F> +} + /// An exponent for norms. /// -// Just a collection of desirabl attributes for a marker type -pub trait NormExponent : Copy + Send + Sync + 'static {} - +// Just a collection of desirable attributes for a marker type +pub trait NormExponent : Copy + Send + Sync + 'static { + /// Return the norm as a mappin + fn as_mapping<F : Float>(self) -> NormMapping<F, Self> { + NormMapping{ exponent : self, _phantoms : PhantomData } + } +} /// Exponent type for the 1-[`Norm`]. #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] @@ -37,6 +50,15 @@ pub struct L21; impl NormExponent for L21 {} +/// Norms for pairs (a, b). ‖(a,b)‖ = ‖(‖a‖_A, ‖b‖_B)‖_J +/// For use with [`crate::direct_product::Pair`] +#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] +pub struct PairNorm<A, B, J>(pub A, pub B, pub J); + +impl<A, B, J> NormExponent for PairNorm<A, B, J> +where A : NormExponent, B : NormExponent, J : NormExponent {} + + /// A Huber/Moreau–Yosida smoothed [`L1`] norm. (Not a norm itself.) /// /// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher @@ -82,9 +104,9 @@ } /// Trait for distances with respect to a norm. -pub trait Dist<F : Num, Exponent : NormExponent> : Norm<F, Exponent> { +pub trait Dist<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Space { /// Calculate the distance - fn dist(&self, other : &Self, _p : Exponent) -> F; + fn dist<I : Instance<Self>>(&self, other : I, _p : Exponent) -> F; } /// Trait for Euclidean projections to the `Exponent`-[`Norm`]-ball. @@ -97,7 +119,7 @@ /// /// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity)); /// ``` -pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Euclidean<F> +pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Sized where F : Float { /// Projection of `self` to the `q`-norm-ball of radius ρ. fn proj_ball(mut self, ρ : F, q : Exponent) -> Self { @@ -106,7 +128,7 @@ } /// In-place projection of `self` to the `q`-norm-ball of radius ρ. - fn proj_ball_mut(&mut self, ρ : F, _q : Exponent); + fn proj_ball_mut(&mut self, ρ : F, q : Exponent); } /*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E { @@ -149,8 +171,149 @@ } impl<F : Float, E : Euclidean<F>> Dist<F, HuberL1<F>> for E { - fn dist(&self, other : &Self, huber : HuberL1<F>) -> F { + fn dist<I : Instance<Self>>(&self, other : I, huber : HuberL1<F>) -> F { huber.apply(self.dist2_squared(other)) } } +// impl<F : Float, E : Norm<F, L2>> Norm<F, L21> for Vec<E> { +// fn norm(&self, _l21 : L21) -> F { +// self.iter().map(|e| e.norm(L2)).sum() +// } +// } + +// impl<F : Float, E : Dist<F, L2>> Dist<F, L21> for Vec<E> { +// fn dist<I : Instance<Self>>(&self, other : I, _l21 : L21) -> F { +// self.iter().zip(other.iter()).map(|(e, g)| e.dist(g, L2)).sum() +// } +// } + +impl<E, F, Domain> Mapping<Domain> for NormMapping<F, E> +where + F : Float, + E : NormExponent, + Domain : Space + Norm<F, E>, +{ + type Codomain = F; + + #[inline] + fn apply<I : Instance<Domain>>(&self, x : I) -> F { + x.eval(|r| r.norm(self.exponent)) + } +} + +pub trait Normed<F : Num = f64> : Space + Norm<F, Self::NormExp> { + type NormExp : NormExponent; + + fn norm_exponent(&self) -> Self::NormExp; + + #[inline] + fn norm_(&self) -> F { + self.norm(self.norm_exponent()) + } + + // fn similar_origin(&self) -> Self; + + fn is_zero(&self) -> bool; +} + +pub trait HasDual<F : Num = f64> : Normed<F> { + type DualSpace : Normed<F>; +} + +/// Automatically implemented trait for reflexive spaces +pub trait Reflexive<F : Num = f64> : HasDual<F> +where + Self::DualSpace : HasDual<F, DualSpace = Self> +{ } + +impl<F : Num, X : HasDual<F>> Reflexive<F> for X +where + X::DualSpace : HasDual<F, DualSpace = X> +{ } + +pub trait HasDualExponent : NormExponent { + type DualExp : NormExponent; + + fn dual_exponent(&self) -> Self::DualExp; +} + +impl HasDualExponent for L2 { + type DualExp = L2; + + #[inline] + fn dual_exponent(&self) -> Self::DualExp { + L2 + } +} + +impl HasDualExponent for L1 { + type DualExp = Linfinity; + + #[inline] + fn dual_exponent(&self) -> Self::DualExp { + Linfinity + } +} + + +impl HasDualExponent for Linfinity { + type DualExp = L1; + + #[inline] + fn dual_exponent(&self) -> Self::DualExp { + L1 + } +} + +#[macro_export] +macro_rules! impl_weighted_norm { + ($exponent : ty) => { + impl<C, F, D> Norm<F, Weighted<$exponent, C>> for D + where + F : Float, + D : Norm<F, $exponent>, + C : Constant<Type = F>, + { + fn norm(&self, e : Weighted<$exponent, C>) -> F { + let v = e.weight.value(); + assert!(v > F::ZERO); + v * self.norm(e.base_fn) + } + } + + impl<C : Constant> NormExponent for Weighted<$exponent, C> {} + + impl<C : Constant> HasDualExponent for Weighted<$exponent, C> + where $exponent : HasDualExponent { + type DualExp = Weighted<<$exponent as HasDualExponent>::DualExp, C::Type>; + + fn dual_exponent(&self) -> Self::DualExp { + Weighted { + weight : C::Type::ONE / self.weight.value(), + base_fn : self.base_fn.dual_exponent() + } + } + } + + impl<C, F, T> Projection<F, Weighted<$exponent , C>> for T + where + T : Projection<F, $exponent >, + F : Float, + C : Constant<Type = F>, + { + fn proj_ball(self, ρ : F, q : Weighted<$exponent , C>) -> Self { + self.proj_ball(ρ / q.weight.value(), q.base_fn) + } + + fn proj_ball_mut(&mut self, ρ : F, q : Weighted<$exponent , C>) { + self.proj_ball_mut(ρ / q.weight.value(), q.base_fn) + } + } + } +} + +//impl_weighted_norm!(L1); +//impl_weighted_norm!(L2); +//impl_weighted_norm!(Linfinity); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/operator_arithmetic.rs Mon Feb 03 19:22:16 2025 -0500 @@ -0,0 +1,117 @@ +/*! +Arithmetic of [`Mapping`]s. + */ + +use serde::Serialize; +use crate::types::*; +use crate::instance::{Space, Instance}; +use crate::mapping::{Mapping, DifferentiableImpl, DifferentiableMapping}; + +/// A trait for encoding constant [`Float`] values +pub trait Constant : Copy + Sync + Send + 'static + std::fmt::Debug + Into<Self::Type> { + /// The type of the value + type Type : Float; + /// Returns the value of the constant + fn value(&self) -> Self::Type; +} + +impl<F : Float> Constant for F { + type Type = F; + #[inline] + fn value(&self) -> F { *self } +} + +/// Weighting of a [`Mapping`] by scalar multiplication. +#[derive(Copy,Clone,Debug,Serialize)] +pub struct Weighted<T, C : Constant> { + /// The weight + pub weight : C, + /// The base [`Mapping`] being weighted. + pub base_fn : T, +} + +impl<T, C> Weighted<T, C> +where + C : Constant, +{ + /// Construct from an iterator. + pub fn new(weight : C, base_fn : T) -> Self { + Weighted{ weight, base_fn } + } +} + +impl<'a, T, V, D, F, C> Mapping<D> for Weighted<T, C> +where + F : Float, + D : Space, + T : Mapping<D, Codomain=V>, + V : Space + ClosedMul<F>, + C : Constant<Type=F> +{ + type Codomain = V; + + #[inline] + fn apply<I : Instance<D>>(&self, x : I) -> Self::Codomain { + self.base_fn.apply(x) * self.weight.value() + } +} + +impl<'a, T, V, D, F, C> DifferentiableImpl<D> for Weighted<T, C> +where + F : Float, + D : Space, + T : DifferentiableMapping<D, DerivativeDomain=V>, + V : Space + std::ops::Mul<F, Output=V>, + C : Constant<Type=F> +{ + type Derivative = V; + + #[inline] + fn differential_impl<I : Instance<D>>(&self, x : I) -> Self::Derivative { + self.base_fn.differential(x) * self.weight.value() + } +} + +/// A sum of [`Mapping`]s. +#[derive(Serialize, Debug, Clone)] +pub struct MappingSum<M>(Vec<M>); + +impl< M> MappingSum<M> { + /// Construct from an iterator. + pub fn new<I : IntoIterator<Item = M>>(iter : I) -> Self { + MappingSum(iter.into_iter().collect()) + } + + /// Iterate over the component functions of the sum + pub fn iter(&self) -> std::slice::Iter<'_, M> { + self.0.iter() + } +} + +impl<Domain, M> Mapping<Domain> for MappingSum<M> +where + Domain : Space + Clone, + M : Mapping<Domain>, + M::Codomain : std::iter::Sum + Clone +{ + type Codomain = M::Codomain; + + fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain { + let xr = x.ref_instance(); + self.0.iter().map(|c| c.apply(xr)).sum() + } +} + +impl<Domain, M> DifferentiableImpl<Domain> for MappingSum< M> +where + Domain : Space + Clone, + M : DifferentiableMapping<Domain>, + M :: DerivativeDomain : std::iter::Sum +{ + type Derivative = M::DerivativeDomain; + + fn differential_impl<I : Instance<Domain>>(&self, x : I) -> Self::Derivative { + let xr = x.ref_instance(); + self.0.iter().map(|c| c.differential(xr)).sum() + } +}
--- a/src/parallelism.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/parallelism.rs Mon Feb 03 19:22:16 2025 -0500 @@ -16,9 +16,9 @@ Ordering::{Release, Relaxed}, }; -#[cfg(use_custom_thread_pool)] +#[cfg(feature = "use_custom_thread_pool")] type Pool = ThreadPool; -#[cfg(not(use_custom_thread_pool))] +#[cfg(not(feature = "use_custom_thread_pool"))] type Pool = GlobalPool; const ONE : NonZeroUsize = unsafe { NonZeroUsize::new_unchecked(1) }; @@ -27,12 +27,12 @@ static mut POOL : Option<Pool> = None; static INIT: Once = Once::new(); -#[cfg(not(use_custom_thread_pool))] +#[cfg(not(feature = "use_custom_thread_pool"))] mod global_pool { /// This is a nicer way to use the global pool of [`rayon`]. pub struct GlobalPool; - #[cfg(not(use_custom_thread_pool))] + #[cfg(not(feature = "use_custom_thread_pool"))] impl GlobalPool { #[inline] pub fn scope<'scope, OP, R>(&self, op: OP) -> R @@ -44,7 +44,7 @@ } } -#[cfg(not(use_custom_thread_pool))] +#[cfg(not(feature = "use_custom_thread_pool"))] pub use global_pool::GlobalPool; /// Set the number of threads. @@ -57,10 +57,10 @@ let n = n.get(); set_task_overbudgeting((n + 1) / 2); POOL = if n > 1 { - #[cfg(use_custom_thread_pool)] { + #[cfg(feature = "use_custom_thread_pool")] { Some(ThreadPoolBuilder::new().num_threads(n).build().unwrap()) } - #[cfg(not(use_custom_thread_pool))] { + #[cfg(not(feature = "use_custom_thread_pool"))] { ThreadPoolBuilder::new().num_threads(n).build_global().unwrap(); Some(GlobalPool) } @@ -75,6 +75,7 @@ /// The initial value is 1. Calling [`set_num_threads`] sets this to $m = (n + 1) / 2$, where /// $n$ is the number of threads. pub fn set_task_overbudgeting(m : usize) { + #[allow(static_mut_refs)] unsafe { TASK_OVERBUDGETING.store(m, Relaxed) } } @@ -97,6 +98,7 @@ /// If the number of configured threads is less than 2, this is None. /// The pool has [`num_threads`]` - 1` threads. pub fn thread_pool() -> Option<&'static Pool> { + #[allow(static_mut_refs)] unsafe { POOL.as_ref() } } @@ -146,6 +148,7 @@ /// is `None`, the [global setting][set_task_overbudgeting] is used.§ pub fn init(overbudget : Option<usize>) -> Self { let n = num_threads().get(); + #[allow(static_mut_refs)] let m = overbudget.unwrap_or_else(|| unsafe { TASK_OVERBUDGETING.load(Relaxed) }); if n <= 1 { Self::SingleThreaded
--- a/src/sets.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/sets.rs Mon Feb 03 19:22:16 2025 -0500 @@ -3,19 +3,19 @@ */ use std::ops::{RangeFull,RangeFrom,Range,RangeInclusive,RangeTo,RangeToInclusive}; -use std::marker::PhantomData; use crate::types::*; use crate::loc::Loc; -use crate::euclidean::Dot; +use crate::euclidean::Euclidean; +use crate::instance::{Space, Instance}; use serde::Serialize; -mod cube; +pub mod cube; pub use cube::Cube; /// Trait for arbitrary sets. The parameter `U` is the element type. -pub trait Set<U> where U : ?Sized { +pub trait Set<U> where U : Space { /// Check for element containment - fn contains(&self, item : &U) -> bool; + fn contains<I : Instance<U>>(&self, item : I) -> bool; } /// Additional ordering (besides [`PartialOrd`]) of a subfamily of sets: @@ -31,22 +31,27 @@ impl<U, const N : usize> Set<Loc<U, N>> for Cube<U,N> where U : Num + PartialOrd + Sized { - fn contains(&self, item : &Loc<U, N>) -> bool { - self.0.iter().zip(item.iter()).all(|(s, x)| s.contains(x)) + fn contains<I : Instance<Loc<U, N>>>(&self, item : I) -> bool { + self.0.iter().zip(item.ref_instance().iter()).all(|(s, x)| s.contains(x)) } } -impl<U> Set<U> for RangeFull { - fn contains(&self, _item : &U) -> bool { true } +impl<U : Space> Set<U> for RangeFull { + fn contains<I : Instance<U>>(&self, _item : I) -> bool { true } } macro_rules! impl_ranges { ($($range:ident),*) => { $( - impl<U,Idx> Set<U> - for $range<Idx> - where Idx : PartialOrd<U>, U : PartialOrd<Idx> + ?Sized, Idx : PartialOrd { + impl<U,Idx> Set<U> for $range<Idx> + where + Idx : PartialOrd<U>, + U : PartialOrd<Idx> + Space, + Idx : PartialOrd + { #[inline] - fn contains(&self, item : &U) -> bool { $range::contains(&self, item) } + fn contains<I : Instance<U>>(&self, item : I) -> bool { + item.eval(|x| $range::contains(&self, x)) + } } )* } } @@ -57,44 +62,40 @@ /// /// The halfspace is $H = \\{ t v + a \mid a^⊤ v = 0 \\}$, where $v$ is the orthogonal /// vector and $t$ the offset. -/// -/// `U` is the element type, `F` the floating point number type, and `A` the type of the -/// orthogonal (dual) vectors. They need implement [`Dot<U, F>`]. #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] -pub struct Halfspace<A, F, U> where A : Dot<U, F>, F : Float { +pub struct Halfspace<A, F> where A : Euclidean<F>, F : Float { pub orthogonal : A, pub offset : F, - _phantom : PhantomData<U>, } -impl<A,F,U> Halfspace<A,F,U> where A : Dot<U,F>, F : Float { +impl<A,F> Halfspace<A,F> where A : Euclidean<F>, F : Float { #[inline] pub fn new(orthogonal : A, offset : F) -> Self { - Halfspace{ orthogonal : orthogonal, offset : offset, _phantom : PhantomData } + Halfspace{ orthogonal : orthogonal, offset : offset } } } /// Trait for generating a halfspace spanned by another set `Self` of elements of type `U`. -pub trait SpannedHalfspace<F, U> where F : Float { +pub trait SpannedHalfspace<F> where F : Float { /// Type of the orthogonal vector describing the halfspace. - type A : Dot<U, F>; + type A : Euclidean<F>; /// Returns the halfspace spanned by this set. - fn spanned_halfspace(&self) -> Halfspace<Self::A, F, U>; + fn spanned_halfspace(&self) -> Halfspace<Self::A, F>; } // TODO: Gram-Schmidt for higher N. -impl<F : Float> SpannedHalfspace<F,Loc<F, 1>> for [Loc<F, 1>; 2] { +impl<F : Float> SpannedHalfspace<F> for [Loc<F, 1>; 2] { type A = Loc<F, 1>; - fn spanned_halfspace(&self) -> Halfspace<Self::A, F, Loc<F, 1>> { + fn spanned_halfspace(&self) -> Halfspace<Self::A, F> { let (x0, x1) = (self[0], self[1]); Halfspace::new(x1-x0, x0[0]) } } // TODO: Gram-Schmidt for higher N. -impl<F : Float> SpannedHalfspace<F,Loc<F, 2>> for [Loc<F, 2>; 2] { +impl<F : Float> SpannedHalfspace<F> for [Loc<F, 2>; 2] { type A = Loc<F, 2>; - fn spanned_halfspace(&self) -> Halfspace<Self::A, F, Loc<F, 2>> { + fn spanned_halfspace(&self) -> Halfspace<Self::A, F> { let (x0, x1) = (&self[0], &self[1]); let d = x1 - x0; let orthog = loc![d[1], -d[0]]; // We don't normalise for efficiency @@ -103,19 +104,29 @@ } } -impl<A,U,F> Set<U> for Halfspace<A,F,U> where A : Dot<U,F>, F : Float { +impl<A,F> Set<A> for Halfspace<A,F> +where + A : Euclidean<F>, + F : Float, +{ #[inline] - fn contains(&self, item : &U) -> bool { + fn contains<I : Instance<A>>(&self, item : I) -> bool { self.orthogonal.dot(item) >= self.offset } } /// Polygons defined by `N` `Halfspace`s. #[derive(Clone,Copy,Debug,Eq,PartialEq)] -pub struct NPolygon<A, F, U, const N : usize>(pub [Halfspace<A,F,U>; N]) where A : Dot<U,F>, F : Float; +pub struct NPolygon<A, F, const N : usize>(pub [Halfspace<A,F>; N]) +where A : Euclidean<F>, F : Float; -impl<A,U,F,const N : usize> Set<U> for NPolygon<A,F,U,N> where A : Dot<U,F>, F : Float { - fn contains(&self, item : &U) -> bool { - self.0.iter().all(|halfspace| halfspace.contains(item)) +impl<A,F,const N : usize> Set<A> for NPolygon<A,F,N> +where + A : Euclidean<F>, + F : Float, +{ + fn contains<I : Instance<A>>(&self, item : I) -> bool { + let r = item.ref_instance(); + self.0.iter().all(|halfspace| halfspace.contains(r)) } }
--- a/src/sets/cube.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/sets/cube.rs Mon Feb 03 19:22:16 2025 -0500 @@ -333,7 +333,9 @@ impl_common!(u8 u16 u32 u64 u128 usize i8 i16 i32 i64 i128 isize, min, max); + // Any NaN yields NaN +#[cfg(feature = "nightly")] impl_common!(f32 f64, minimum, maximum); impl<U : Num, const N : usize> std::ops::Index<usize> for Cube<U, N> {
--- a/src/types.rs Tue Feb 20 12:33:16 2024 -0500 +++ b/src/types.rs Mon Feb 03 19:22:16 2025 -0500 @@ -13,6 +13,14 @@ pub use num_traits::Float as NumTraitsFloat; // needed to re-export functions. pub use num_traits::cast::AsPrimitive; +pub use simba::scalar::{ + ClosedAdd, ClosedAddAssign, + ClosedSub, ClosedSubAssign, + ClosedMul, ClosedMulAssign, + ClosedDiv, ClosedDivAssign, + ClosedNeg +}; + /// Typical integer type #[allow(non_camel_case_types)] pub type int = i64; @@ -57,7 +65,8 @@ + CastFrom<u128> + CastFrom<usize> + CastFrom<i8> + CastFrom<i16> + CastFrom<i32> + CastFrom<i64> + CastFrom<i128> + CastFrom<isize> - + CastFrom<f32> + CastFrom<f64> { + + CastFrom<f32> + CastFrom<f64> + + crate::instance::Space { const ZERO : Self; const ONE : Self; @@ -84,6 +93,7 @@ const SQRT_2 : Self; const INFINITY : Self; const NEG_INFINITY : Self; + const NAN : Self; const FRAC_2_SQRT_PI : Self; } @@ -133,6 +143,7 @@ const SQRT_2 : Self = std::f64::consts::SQRT_2; const INFINITY : Self = std::f64::INFINITY; const NEG_INFINITY : Self = std::f64::NEG_INFINITY; + const NAN : Self = std::f64::NAN; const FRAC_2_SQRT_PI : Self = std::f64::consts::FRAC_2_SQRT_PI; } @@ -150,6 +161,7 @@ const SQRT_2 : Self = std::f32::consts::SQRT_2; const INFINITY : Self = std::f32::INFINITY; const NEG_INFINITY : Self = std::f32::NEG_INFINITY; + const NAN : Self = std::f32::NAN; const FRAC_2_SQRT_PI : Self = std::f32::consts::FRAC_2_SQRT_PI; }
