# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1738628536 18000
# Node ID b3c35d16affeb46d27d4b70e76199db83cb96eee
# Parent  d14c877e14b70dd4c1ddc9436e4095e84964ef20# Parent  8513a10c5dd937c3fed7159cf9ca04ef46c2a981
merge dev to default

diff -r d14c877e14b7 -r b3c35d16affe Cargo.lock
--- 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",
diff -r d14c877e14b7 -r b3c35d16affe Cargo.toml
--- 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
+
diff -r d14c877e14b7 -r b3c35d16affe README.md
--- 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
 
diff -r d14c877e14b7 -r b3c35d16affe rust-toolchain.toml
--- 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"
diff -r d14c877e14b7 -r b3c35d16affe src/bisection_tree/aggregator.rs
--- 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
     }
 }
 
diff -r d14c877e14b7 -r b3c35d16affe src/bisection_tree/btfn.rs
--- 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> {
 
diff -r d14c877e14b7 -r b3c35d16affe src/bisection_tree/either.rs
--- 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>
diff -r d14c877e14b7 -r b3c35d16affe src/bisection_tree/refine.rs
--- 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"),
diff -r d14c877e14b7 -r b3c35d16affe src/bisection_tree/support.rs
--- 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;
diff -r d14c877e14b7 -r b3c35d16affe src/collection.rs
--- /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);
diff -r d14c877e14b7 -r b3c35d16affe src/convex.rs
--- /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()
+    }
+}
diff -r d14c877e14b7 -r b3c35d16affe src/direct_product.rs
--- /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>;
+}
diff -r d14c877e14b7 -r b3c35d16affe src/discrete_gradient.rs
--- /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);
+
+    }
+}
diff -r d14c877e14b7 -r b3c35d16affe src/error.rs
--- 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<()>;
diff -r d14c877e14b7 -r b3c35d16affe src/euclidean.rs
--- 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()
     }
 
diff -r d14c877e14b7 -r b3c35d16affe src/fe_model/p2_local_model.rs
--- 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)
     }
 }
 
diff -r d14c877e14b7 -r b3c35d16affe src/instance.rs
--- /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
+    }
+}
diff -r d14c877e14b7 -r b3c35d16affe src/iterate.rs
--- 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>>())
+        }
+    }
+
 }
diff -r d14c877e14b7 -r b3c35d16affe src/lib.rs
--- 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::*;
diff -r d14c877e14b7 -r b3c35d16affe src/linops.rs
--- 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);
+
diff -r d14c877e14b7 -r b3c35d16affe src/linsolve.rs
--- 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`,
diff -r d14c877e14b7 -r b3c35d16affe src/loc.rs
--- 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;
     }
 }
+
diff -r d14c877e14b7 -r b3c35d16affe src/logger.rs
--- 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> {
diff -r d14c877e14b7 -r b3c35d16affe src/mapping.rs
--- 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))
+    }
+}
diff -r d14c877e14b7 -r b3c35d16affe src/maputil.rs
--- 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 {
diff -r d14c877e14b7 -r b3c35d16affe src/metaprogramming.rs
--- /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 };
+}
+
diff -r d14c877e14b7 -r b3c35d16affe src/nalgebra_support.rs
--- 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())
     }
 }
 
diff -r d14c877e14b7 -r b3c35d16affe src/norms.rs
--- 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);
+
diff -r d14c877e14b7 -r b3c35d16affe src/operator_arithmetic.rs
--- /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()
+    }
+}
diff -r d14c877e14b7 -r b3c35d16affe src/parallelism.rs
--- 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
diff -r d14c877e14b7 -r b3c35d16affe src/sets.rs
--- 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))
     }
 }
diff -r d14c877e14b7 -r b3c35d16affe src/sets/cube.rs
--- 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> {
diff -r d14c877e14b7 -r b3c35d16affe src/types.rs
--- 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;
 }