diff --git a/Cargo.lock b/Cargo.lock
index 2bb30dbc..43de9609 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1,12 +1,12 @@
 # This file is automatically @generated by Cargo.
 # It is not intended for manual editing.
-version = 3
+version = 4
 
 [[package]]
 name = "atomic_float"
-version = "0.1.0"
+version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "62af46d040ba9df09edc6528dae9d8e49f5f3e82f55b7d2ec31a733c38dbc49d"
+checksum = "628d228f918ac3b82fe590352cc719d30664a0c13ca3a60266fe02c7132d480a"
 
 [[package]]
 name = "autocfg"
@@ -16,9 +16,9 @@ checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26"
 
 [[package]]
 name = "bitflags"
-version = "2.6.0"
+version = "2.8.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de"
+checksum = "8f68f53c83ab957f72c32642f3868eec03eb974d1fb82e453128456482613d36"
 
 [[package]]
 name = "byteorder"
@@ -48,9 +48,9 @@ dependencies = [
 
 [[package]]
 name = "crossbeam-deque"
-version = "0.8.5"
+version = "0.8.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "613f8cc01fe9cf1a3eb3d7f488fd2fa8388403e97039e2f73692932e291a770d"
+checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51"
 dependencies = [
  "crossbeam-epoch",
  "crossbeam-utils",
@@ -67,9 +67,9 @@ dependencies = [
 
 [[package]]
 name = "crossbeam-utils"
-version = "0.8.20"
+version = "0.8.21"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "22ec99545bb0ed0ea7bb9b8e1e9122ea386ff8a48c0922e43f36d45ab09e0e80"
+checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28"
 
 [[package]]
 name = "either"
@@ -85,32 +85,39 @@ checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5"
 
 [[package]]
 name = "fixedbitset"
-version = "0.4.2"
+version = "0.5.7"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0ce7134b9999ecaf8bcd65542e436736ef32ddca1b3e06094cb6ec5755203b80"
+checksum = "1d674e81391d1e1ab681a28d99df07927c6d4aa5b027d7da16ba32d1d21ecd99"
 
 [[package]]
 name = "getrandom"
-version = "0.2.15"
+version = "0.3.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c4567c8db10ae91089c99af84c68c38da3ec2f087c3f82960bcdbf3656b6f4d7"
+checksum = "43a49c392881ce6d5c3b8cb70f98717b7c07aabbdff06687b9030dbfbe2725f8"
 dependencies = [
  "cfg-if",
  "libc",
  "wasi",
+ "windows-targets",
 ]
 
 [[package]]
 name = "hashbrown"
-version = "0.15.0"
+version = "0.15.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1e087f84d4f86bf4b218b927129862374b72199ae7d8657835f1e89000eea4fb"
+checksum = "bf151400ff0baff5465007dd2f3e717f3fe502074ca563069ce3a6629d07b289"
+
+[[package]]
+name = "heck"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
 
 [[package]]
 name = "indexmap"
-version = "2.6.0"
+version = "2.7.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "707907fe3c25f5424cce2cb7e1cbcafee6bdbe735ca90ef77c29e84591e5b9da"
+checksum = "8c9c992b02b5b4c94ea26e32fe5bccb7aa7d9f390ab5c1221ff895bc7ea8b652"
 dependencies = [
  "equivalent",
  "hashbrown",
@@ -118,37 +125,30 @@ dependencies = [
 
 [[package]]
 name = "indoc"
-version = "1.0.9"
+version = "2.0.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "bfa799dd5ed20a7e349f3b4639aa80d74549c81716d9ec4f994c9b5815598306"
+checksum = "b248f5224d1d606005e02c97f5aa4e88eeb230488bcc03bc9ca4d7991399f2b5"
 
 [[package]]
 name = "inventory"
-version = "0.3.15"
+version = "0.3.17"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f958d3d68f4167080a18141e10381e7634563984a537f2a49a30fd8e53ac5767"
+checksum = "3b31349d02fe60f80bbbab1a9402364cad7460626d6030494b08ac4a2075bf81"
+dependencies = [
+ "rustversion",
+]
 
 [[package]]
 name = "libc"
-version = "0.2.161"
+version = "0.2.169"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8e9489c2807c139ffd9c1794f4af0ebe86a828db53ecdc7fea2111d0fed085d1"
+checksum = "b5aba8db14291edd000dfcc4d620c7ebfb122c613afb886ca8803fa4e128a20a"
 
 [[package]]
 name = "libm"
-version = "0.2.8"
+version = "0.2.11"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058"
-
-[[package]]
-name = "lock_api"
-version = "0.4.12"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "07af8b9cdd281b7915f413fa73f29ebd5d55d0d3f0155584dade1ff18cea1b17"
-dependencies = [
- "autocfg",
- "scopeguard",
-]
+checksum = "8355be11b20d696c8f18f6cc018c4e372165b1fa8126cef092399c9951984ffa"
 
 [[package]]
 name = "matrixmultiply"
@@ -162,23 +162,25 @@ dependencies = [
 
 [[package]]
 name = "memoffset"
-version = "0.8.0"
+version = "0.9.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d61c719bcfbcf5d62b3a09efa6088de8c54bc0bfcd3ea7ae39fcc186108b8de1"
+checksum = "488016bfae457b036d996092f6cb448677611ce4449e970ceaf42695203f218a"
 dependencies = [
  "autocfg",
 ]
 
 [[package]]
 name = "ndarray"
-version = "0.15.6"
+version = "0.16.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32"
+checksum = "882ed72dce9365842bf196bdeedf5055305f11fc8c03dee7bb0194a6cad34841"
 dependencies = [
  "matrixmultiply",
  "num-complex",
  "num-integer",
  "num-traits",
+ "portable-atomic",
+ "portable-atomic-util",
  "rawpointer",
 ]
 
@@ -212,9 +214,9 @@ dependencies = [
 
 [[package]]
 name = "numpy"
-version = "0.18.0"
+version = "0.23.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "96b0fee4571867d318651c24f4a570c3f18408cf95f16ccb576b3ce85496a46e"
+checksum = "b94caae805f998a07d33af06e6a3891e38556051b8045c615470a71590e13e78"
 dependencies = [
  "libc",
  "ndarray",
@@ -232,36 +234,28 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "1261fe7e33c73b354eab43b1273a57c8f967d0391e80353e51f764ac02cf6775"
 
 [[package]]
-name = "parking_lot"
-version = "0.12.3"
+name = "petgraph"
+version = "0.7.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f1bf18183cf54e8d6059647fc3063646a1801cf30896933ec2311622cc4b9a27"
+checksum = "3672b37090dbd86368a4145bc067582552b29c27377cad4e0a306c97f9bd7772"
 dependencies = [
- "lock_api",
- "parking_lot_core",
+ "fixedbitset",
+ "indexmap",
 ]
 
 [[package]]
-name = "parking_lot_core"
-version = "0.9.10"
+name = "portable-atomic"
+version = "1.10.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1e401f977ab385c9e4e3ab30627d6f26d00e2c73eef317493c4ec6d468726cf8"
-dependencies = [
- "cfg-if",
- "libc",
- "redox_syscall",
- "smallvec",
- "windows-targets",
-]
+checksum = "280dc24453071f1b63954171985a0b0d30058d287960968b9b2aca264c8d4ee6"
 
 [[package]]
-name = "petgraph"
-version = "0.6.5"
+name = "portable-atomic-util"
+version = "0.2.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b4c5cc86750666a3ed20bdaf5ca2a0344f9c67674cae0515bec2da16fbaa47db"
+checksum = "d8a2f0d8d040d7848a709caf78912debcc3f33ee4b3cac47d73d1e1069e83507"
 dependencies = [
- "fixedbitset",
- "indexmap",
+ "portable-atomic",
 ]
 
 [[package]]
@@ -270,30 +264,31 @@ version = "0.2.20"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "77957b295656769bb8ad2b6a6b09d897d94f05c41b069aede1fcdaa675eaea04"
 dependencies = [
- "zerocopy",
+ "zerocopy 0.7.35",
 ]
 
 [[package]]
 name = "proc-macro2"
-version = "1.0.89"
+version = "1.0.93"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f139b0662de085916d1fb67d2b4169d1addddda1919e696f3252b740b629986e"
+checksum = "60946a68e5f9d28b0dc1c21bb8a97ee7d018a8b322fa57838ba31cc878e22d99"
 dependencies = [
  "unicode-ident",
 ]
 
 [[package]]
 name = "pyo3"
-version = "0.18.3"
+version = "0.23.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e3b1ac5b3731ba34fdaa9785f8d74d17448cd18f30cf19e0c7e7b1fdb5272109"
+checksum = "57fe09249128b3173d092de9523eaa75136bf7ba85e0d69eca241c7939c933cc"
 dependencies = [
  "cfg-if",
  "indoc",
  "inventory",
  "libc",
  "memoffset",
- "parking_lot",
+ "once_cell",
+ "portable-atomic",
  "pyo3-build-config",
  "pyo3-ffi",
  "pyo3-macros",
@@ -302,9 +297,9 @@ dependencies = [
 
 [[package]]
 name = "pyo3-build-config"
-version = "0.18.3"
+version = "0.23.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9cb946f5ac61bb61a5014924910d936ebd2b23b705f7a4a3c40b05c720b079a3"
+checksum = "1cd3927b5a78757a0d71aa9dff669f903b1eb64b54142a9bd9f757f8fde65fd7"
 dependencies = [
  "once_cell",
  "target-lexicon",
@@ -312,9 +307,9 @@ dependencies = [
 
 [[package]]
 name = "pyo3-ffi"
-version = "0.18.3"
+version = "0.23.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fd4d7c5337821916ea2a1d21d1092e8443cf34879e53a0ac653fbb98f44ff65c"
+checksum = "dab6bb2102bd8f991e7749f130a70d05dd557613e39ed2deeee8e9ca0c4d548d"
 dependencies = [
  "libc",
  "pyo3-build-config",
@@ -322,52 +317,54 @@ dependencies = [
 
 [[package]]
 name = "pyo3-macros"
-version = "0.18.3"
+version = "0.23.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a9d39c55dab3fc5a4b25bbd1ac10a2da452c4aca13bb450f22818a002e29648d"
+checksum = "91871864b353fd5ffcb3f91f2f703a22a9797c91b9ab497b1acac7b07ae509c7"
 dependencies = [
  "proc-macro2",
  "pyo3-macros-backend",
  "quote",
- "syn 1.0.109",
+ "syn",
 ]
 
 [[package]]
 name = "pyo3-macros-backend"
-version = "0.18.3"
+version = "0.23.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "97daff08a4c48320587b5224cc98d609e3c27b6d437315bd40b605c98eeb5918"
+checksum = "43abc3b80bc20f3facd86cd3c60beed58c3e2aa26213f3cda368de39c60a27e4"
 dependencies = [
+ "heck",
  "proc-macro2",
+ "pyo3-build-config",
  "quote",
- "syn 1.0.109",
+ "syn",
 ]
 
 [[package]]
 name = "quote"
-version = "1.0.37"
+version = "1.0.38"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af"
+checksum = "0e4dccaaaf89514f546c693ddc140f729f958c247918a13380cccc6078391acc"
 dependencies = [
  "proc-macro2",
 ]
 
 [[package]]
 name = "rand"
-version = "0.8.5"
+version = "0.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404"
+checksum = "3779b94aeb87e8bd4e834cee3650289ee9e0d5677f976ecdb6d219e5f4f6cd94"
 dependencies = [
- "libc",
  "rand_chacha",
  "rand_core",
+ "zerocopy 0.8.14",
 ]
 
 [[package]]
 name = "rand_chacha"
-version = "0.3.1"
+version = "0.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88"
+checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb"
 dependencies = [
  "ppv-lite86",
  "rand_core",
@@ -375,18 +372,19 @@ dependencies = [
 
 [[package]]
 name = "rand_core"
-version = "0.6.4"
+version = "0.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
+checksum = "b08f3c9802962f7e1b25113931d94f43ed9725bebc59db9d0c3e9a23b67e15ff"
 dependencies = [
  "getrandom",
+ "zerocopy 0.8.14",
 ]
 
 [[package]]
 name = "rand_distr"
-version = "0.4.3"
+version = "0.5.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31"
+checksum = "ddc3b5afe4c995c44540865b8ca5c52e6a59fa362da96c5d30886930ddc8da1c"
 dependencies = [
  "num-traits",
  "rand",
@@ -418,49 +416,23 @@ dependencies = [
  "crossbeam-utils",
 ]
 
-[[package]]
-name = "redox_syscall"
-version = "0.5.7"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9b6dfecf2c74bce2466cabf93f6664d6998a69eb21e39f4207930065b27b771f"
-dependencies = [
- "bitflags",
-]
-
 [[package]]
 name = "rustc-hash"
-version = "1.1.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2"
-
-[[package]]
-name = "scopeguard"
-version = "1.2.0"
+version = "2.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49"
+checksum = "c7fb8039b3032c191086b10f11f319a6e99e1e82889c5cc6046f515c9db1d497"
 
 [[package]]
-name = "smallvec"
-version = "1.13.2"
+name = "rustversion"
+version = "1.0.19"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67"
+checksum = "f7c45b9784283f1b2e7fb61b42047c2fd678ef0960d4f6f1eba131594cc369d4"
 
 [[package]]
 name = "syn"
-version = "1.0.109"
+version = "2.0.96"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237"
-dependencies = [
- "proc-macro2",
- "quote",
- "unicode-ident",
-]
-
-[[package]]
-name = "syn"
-version = "2.0.85"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5023162dfcd14ef8f32034d8bcd4cc5ddc61ef7a247c024a33e24e1f24d21b56"
+checksum = "d5d0adab1ae378d7f53bdebc67a39f1f151407ef230f0ce2883572f5d8985c80"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -475,21 +447,24 @@ checksum = "61c41af27dd6d1e27b1b16b489db798443478cef1f06a660c96db617ba5de3b1"
 
 [[package]]
 name = "unicode-ident"
-version = "1.0.13"
+version = "1.0.16"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e91b56cd4cadaeb79bbf1a5645f6b4f8dc5bde8834ad5894a8db35fda9efa1fe"
+checksum = "a210d160f08b701c8721ba1c726c11662f877ea6b7094007e1ca9a1041945034"
 
 [[package]]
 name = "unindent"
-version = "0.1.11"
+version = "0.2.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e1766d682d402817b5ac4490b3c3002d91dfa0d22812f341609f97b08757359c"
+checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce"
 
 [[package]]
 name = "wasi"
-version = "0.11.0+wasi-snapshot-preview1"
+version = "0.13.3+wasi-0.2.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423"
+checksum = "26816d2e1a4a36a2940b96c5296ce403917633dff8f3440e9b236ed6f6bacad2"
+dependencies = [
+ "wit-bindgen-rt",
+]
 
 [[package]]
 name = "windows-targets"
@@ -555,6 +530,15 @@ version = "0.52.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec"
 
+[[package]]
+name = "wit-bindgen-rt"
+version = "0.33.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3268f3d866458b787f390cf61f4bbb563b922d091359f9608842999eaee3943c"
+dependencies = [
+ "bitflags",
+]
+
 [[package]]
 name = "zerocopy"
 version = "0.7.35"
@@ -562,7 +546,16 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "1b9b4fd18abc82b8136838da5d50bae7bdea537c574d8dc1a34ed098d6c166f0"
 dependencies = [
  "byteorder",
- "zerocopy-derive",
+ "zerocopy-derive 0.7.35",
+]
+
+[[package]]
+name = "zerocopy"
+version = "0.8.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a367f292d93d4eab890745e75a778da40909cab4d6ff8173693812f79c4a2468"
+dependencies = [
+ "zerocopy-derive 0.8.14",
 ]
 
 [[package]]
@@ -573,5 +566,16 @@ checksum = "fa4f8080344d4671fb4e831a13ad1e68092748387dfc4f55e356242fae12ce3e"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.85",
+ "syn",
+]
+
+[[package]]
+name = "zerocopy-derive"
+version = "0.8.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d3931cb58c62c13adec22e38686b559c86a30565e16ad6e8510a337cedc611e1"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn",
 ]
diff --git a/Cargo.toml b/Cargo.toml
index 24be24fc..2ac49160 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -16,11 +16,11 @@ name = "cityseer"
 crate-type = ["cdylib"]
 
 [dependencies]
-atomic_float = "0.1.0"
-ndarray = "0.15.6"
-numpy = "0.18.0"
-petgraph = "0.6.3"
-pyo3 = { version="0.18.3", features=["multiple-pymethods"]}
-rand = "0.8.5"
-rand_distr = "0.4.3"
+atomic_float = "1.1.0"
+ndarray = "0.16.1"
+numpy = "0.23.0"
+petgraph = "0.7.1"
+pyo3 = { version = "0.23.4", features = ["multiple-pymethods"] }
+rand = "0.9.0"
+rand_distr = "0.5.0"
 rayon = "1.7.0"
diff --git a/README.md b/README.md
index ccbbd1f7..c101624e 100644
--- a/README.md
+++ b/README.md
@@ -36,5 +36,5 @@ The `cityseer-api` `Python` package addresses a range of issues specific to comp
 
 ## Development
 
-`brew install pdm rust rust-analyzer rustfmt`
-`pdm install`
+`brew install uv rust rust-analyzer rustfmt`
+`uv sync`
diff --git a/docs/public/images/assignment.png b/docs/public/images/assignment.png
index b94e1935..cd025675 100644
Binary files a/docs/public/images/assignment.png and b/docs/public/images/assignment.png differ
diff --git a/docs/public/images/assignment_decomposed.png b/docs/public/images/assignment_decomposed.png
index 3f78d58f..bc4df6c0 100644
Binary files a/docs/public/images/assignment_decomposed.png and b/docs/public/images/assignment_decomposed.png differ
diff --git a/docs/public/images/assignment_plot.png b/docs/public/images/assignment_plot.png
index 7fee2c2a..af8fc2c7 100644
Binary files a/docs/public/images/assignment_plot.png and b/docs/public/images/assignment_plot.png differ
diff --git a/docs/public/images/betas.png b/docs/public/images/betas.png
index 8c80efe6..337eb34f 100644
Binary files a/docs/public/images/betas.png and b/docs/public/images/betas.png differ
diff --git a/docs/public/images/graph.png b/docs/public/images/graph.png
index 862d8963..29de452f 100644
Binary files a/docs/public/images/graph.png and b/docs/public/images/graph.png differ
diff --git a/docs/public/images/graph_clean.png b/docs/public/images/graph_clean.png
index d78df366..5dce7e6e 100644
Binary files a/docs/public/images/graph_clean.png and b/docs/public/images/graph_clean.png differ
diff --git a/docs/public/images/graph_colour.png b/docs/public/images/graph_colour.png
index 25e803e0..68ae3265 100644
Binary files a/docs/public/images/graph_colour.png and b/docs/public/images/graph_colour.png differ
diff --git a/docs/public/images/graph_decomposed.png b/docs/public/images/graph_decomposed.png
index e45f7c9b..b8f87d9f 100644
Binary files a/docs/public/images/graph_decomposed.png and b/docs/public/images/graph_decomposed.png differ
diff --git a/docs/public/images/graph_dual.png b/docs/public/images/graph_dual.png
index 169ad379..650cc9de 100644
Binary files a/docs/public/images/graph_dual.png and b/docs/public/images/graph_dual.png differ
diff --git a/docs/public/images/graph_example.png b/docs/public/images/graph_example.png
index fc145bf4..e461305c 100644
Binary files a/docs/public/images/graph_example.png and b/docs/public/images/graph_example.png differ
diff --git a/docs/public/images/graph_raw.png b/docs/public/images/graph_raw.png
index a88292e1..e9364dc8 100644
Binary files a/docs/public/images/graph_raw.png and b/docs/public/images/graph_raw.png differ
diff --git a/docs/public/images/graph_simple.png b/docs/public/images/graph_simple.png
index 9f5a4bf8..4afc27e2 100644
Binary files a/docs/public/images/graph_simple.png and b/docs/public/images/graph_simple.png differ
diff --git a/docs/public/images/intro_mixed_uses.png b/docs/public/images/intro_mixed_uses.png
index c158d3bb..517f2177 100644
Binary files a/docs/public/images/intro_mixed_uses.png and b/docs/public/images/intro_mixed_uses.png differ
diff --git a/docs/public/images/intro_segment_harmonic.png b/docs/public/images/intro_segment_harmonic.png
index faafb029..220f4641 100644
Binary files a/docs/public/images/intro_segment_harmonic.png and b/docs/public/images/intro_segment_harmonic.png differ
diff --git a/pyproject.toml b/pyproject.toml
index 3613cc6a..3d726509 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "cityseer"
-version = '4.17.5'
+version = '4.18.0'
 description = "Computational tools for network-based pedestrian-scale urban analysis"
 readme = "README.md"
 requires-python = ">=3.10, <3.14"
@@ -43,12 +43,12 @@ dependencies = [
   "requests>=2.27.1",
   "scikit-learn>=1.0.2",
   "tqdm>=4.63.1",
-  "shapely>=2.0.0",
-  "numpy>=2.0.0",
-  "geopandas>=1.0.0",
+  "shapely>=2.0.6",
+  "numpy>=2.2.2",
+  "geopandas>=1.0.1",
   "rasterio>=1.3.9",
   "fiona>=1.9.6",
-  "osmnx>=2.0.0",
+  "osmnx>=2.0.1",
 ]
 
 [project.urls]
diff --git a/pysrc/cityseer/metrics/layers.py b/pysrc/cityseer/metrics/layers.py
index 0c2e039c..449db20b 100644
--- a/pysrc/cityseer/metrics/layers.py
+++ b/pysrc/cityseer/metrics/layers.py
@@ -480,7 +480,7 @@ def compute_mixed_uses(
 
 def compute_stats(
     data_gdf: gpd.GeoDataFrame,
-    stats_column_label: str,
+    stats_column_labels: list[str],
     nodes_gdf: gpd.GeoDataFrame,
     network_structure: rustalgos.NetworkStructure,
     max_netw_assign_dist: int = 400,
@@ -507,8 +507,8 @@ def compute_stats(
         representing data points. The coordinates of data points should correspond as precisely as possible to the
         location of the feature in space; or, in the case of buildings, should ideally correspond to the location of the
         building entrance.
-    stats_column_label: str
-        The column label corresponding to the column in `data_gdf` from which to take numerical information.
+    stats_column_labels: list[str]
+        The column labels corresponding to the columns in `data_gdf` from which to take numerical information.
     nodes_gdf
         A [`GeoDataFrame`](https://geopandas.org/en/stable/docs/user_guide/data_structures.html#geodataframe)
         representing nodes. Best generated with the
@@ -598,18 +598,20 @@ def compute_stats(
     :::
 
     """
-    if stats_column_label not in data_gdf.columns:
-        raise ValueError("The specified numerical stats column name can't be found in the GeoDataFrame.")
     data_map, data_gdf = assign_gdf_to_network(data_gdf, network_structure, max_netw_assign_dist, data_id_col)
     if not config.QUIET_MODE:
         logger.info("Computing statistics.")
     # extract landuses
-    stats_map: dict[str, float] = data_gdf[stats_column_label].to_dict()  # type: ignore
+    stats_maps = []
+    for stats_column_label in stats_column_labels:
+        if stats_column_label not in data_gdf.columns:
+            raise ValueError("The specified numerical stats column name can't be found in the GeoDataFrame.")
+        stats_maps.append(data_gdf[stats_column_label].to_dict())  # type: ignore)
     # stats
     partial_func = partial(
         data_map.stats,
         network_structure=network_structure,
-        numerical_map=stats_map,
+        numerical_maps=stats_maps,
         distances=distances,
         betas=betas,
         angular=angular,
@@ -621,26 +623,27 @@ def compute_stats(
     result = config.wrap_progress(total=network_structure.node_count(), rust_struct=data_map, partial_func=partial_func)
     # unpack the numerical arrays
     distances, betas = rustalgos.pair_distances_and_betas(distances, betas)
-    for dist_key in distances:
-        k = config.prep_gdf_key(f"{stats_column_label}_sum", dist_key, angular=angular, weighted=False)
-        nodes_gdf[k] = result.sum[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_sum", dist_key, angular=angular, weighted=True)
-        nodes_gdf[k] = result.sum_wt[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_mean", dist_key, angular=angular, weighted=False)
-        nodes_gdf[k] = result.mean[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_mean", dist_key, angular=angular, weighted=True)
-        nodes_gdf[k] = result.mean_wt[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_count", dist_key, angular=angular, weighted=False)
-        nodes_gdf[k] = result.count[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_count", dist_key, angular=angular, weighted=True)
-        nodes_gdf[k] = result.count_wt[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_var", dist_key, angular=angular, weighted=False)
-        nodes_gdf[k] = result.variance[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_var", dist_key, angular=angular, weighted=True)
-        nodes_gdf[k] = result.variance_wt[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_max", dist_key, angular=angular)
-        nodes_gdf[k] = result.max[dist_key]  # type: ignore
-        k = config.prep_gdf_key(f"{stats_column_label}_min", dist_key, angular=angular)
-        nodes_gdf[k] = result.min[dist_key]  # type: ignore
+    for idx, stats_column_label in enumerate(stats_column_labels):
+        for dist_key in distances:
+            k = config.prep_gdf_key(f"{stats_column_label}_sum", dist_key, angular=angular, weighted=False)
+            nodes_gdf[k] = result[idx].sum[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_sum", dist_key, angular=angular, weighted=True)
+            nodes_gdf[k] = result[idx].sum_wt[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_mean", dist_key, angular=angular, weighted=False)
+            nodes_gdf[k] = result[idx].mean[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_mean", dist_key, angular=angular, weighted=True)
+            nodes_gdf[k] = result[idx].mean_wt[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_count", dist_key, angular=angular, weighted=False)
+            nodes_gdf[k] = result[idx].count[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_count", dist_key, angular=angular, weighted=True)
+            nodes_gdf[k] = result[idx].count_wt[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_var", dist_key, angular=angular, weighted=False)
+            nodes_gdf[k] = result[idx].variance[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_var", dist_key, angular=angular, weighted=True)
+            nodes_gdf[k] = result[idx].variance_wt[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_max", dist_key, angular=angular)
+            nodes_gdf[k] = result[idx].max[dist_key]  # type: ignore
+            k = config.prep_gdf_key(f"{stats_column_label}_min", dist_key, angular=angular)
+            nodes_gdf[k] = result[idx].min[dist_key]  # type: ignore
 
     return nodes_gdf, data_gdf
diff --git a/pysrc/cityseer/rustalgos.pyi b/pysrc/cityseer/rustalgos.pyi
index 52d6d7fe..a37c459c 100644
--- a/pysrc/cityseer/rustalgos.pyi
+++ b/pysrc/cityseer/rustalgos.pyi
@@ -690,7 +690,7 @@ class DataMap:
     def stats(
         self,
         network_structure: NetworkStructure,
-        numerical_map: dict[str, float],
+        numerical_maps: list[dict[str, float]],
         distances: list[int] | None = None,
         betas: list[float] | None = None,
         angular: bool | None = None,
@@ -698,7 +698,7 @@ class DataMap:
         min_threshold_wt: float | None = None,
         jitter_scale: float | None = None,
         pbar_disabled: bool | None = None,
-    ) -> StatsResult: ...
+    ) -> list[StatsResult]: ...
 
 class Viewshed:
     @classmethod
diff --git a/src/centrality.rs b/src/centrality.rs
index bd36f1dd..1a7f799f 100644
--- a/src/centrality.rs
+++ b/src/centrality.rs
@@ -91,6 +91,7 @@ impl NetworkStructure {
     impedance heuristic - which can be different from metres. Distance map in metres is used for defining max
     distances and computing equivalent distance measures.
     */
+    #[pyo3(signature = (src_idx, max_dist, jitter_scale=None))]
     pub fn dijkstra_tree_shortest(
         &self,
         src_idx: usize,
@@ -114,7 +115,7 @@ impl NetworkStructure {
             distance: 0.0,
         });
         // random number generator
-        let mut rng = rand::thread_rng();
+        let mut rng = rand::rng();
         while let Some(NodeDistance { node_idx, .. }) = active.pop() {
             tree_map[node_idx].visited = true;
             visited_nodes.push(node_idx);
@@ -173,7 +174,7 @@ impl NetworkStructure {
                 // inject jitter
                 let mut jitter: f32 = 0.0;
                 if jitter_scale > 0.0 {
-                    jitter = rng.gen::<f32>() * jitter_scale;
+                    jitter = rng.random::<f32>() * jitter_scale;
                 }
                 /*
                 if impedance less than prior distances for this node then update shortest path
@@ -186,6 +187,7 @@ impl NetworkStructure {
         }
         (visited_nodes, tree_map)
     }
+    #[pyo3(signature = (src_idx, max_dist, jitter_scale=None))]
     pub fn dijkstra_tree_simplest(
         &self,
         src_idx: usize,
@@ -210,7 +212,7 @@ impl NetworkStructure {
             distance: 0.0,
         });
         // random number generator
-        let mut rng = rand::thread_rng();
+        let mut rng = rand::rng();
         while let Some(NodeDistance { node_idx, .. }) = active.pop() {
             tree_map[node_idx].visited = true;
             visited_nodes.push(node_idx);
@@ -269,7 +271,7 @@ impl NetworkStructure {
                 // inject jitter
                 let mut jitter: f32 = 0.0;
                 if jitter_scale > 0.0 {
-                    jitter = rng.gen::<f32>() * jitter_scale;
+                    jitter = rng.random::<f32>() * jitter_scale;
                 }
                 /*
                 if impedance less than prior distances for this node then update shortest path
@@ -284,6 +286,7 @@ impl NetworkStructure {
         }
         (visited_nodes, tree_map)
     }
+    #[pyo3(signature = (src_idx, max_dist, jitter_scale=None))]
     pub fn dijkstra_tree_segment(
         &self,
         src_idx: usize,
@@ -309,7 +312,7 @@ impl NetworkStructure {
             distance: 0.0,
         });
         // random number generator
-        let mut rng = rand::thread_rng();
+        let mut rng = rand::rng();
         while let Some(NodeDistance { node_idx, .. }) = active.pop() {
             tree_map[node_idx].visited = true;
             visited_nodes.push(node_idx);
@@ -363,7 +366,7 @@ impl NetworkStructure {
                 // inject jitter
                 let mut jitter: f32 = 0.0;
                 if jitter_scale > 0.0 {
-                    jitter = rng.gen::<f32>() * jitter_scale;
+                    jitter = rng.random::<f32>() * jitter_scale;
                 }
                 /*
                 if impedance less than prior distances for this node then update shortest path
@@ -384,6 +387,15 @@ impl NetworkStructure {
         }
         (visited_nodes, visited_edges, tree_map, edge_map)
     }
+    #[pyo3(signature = (
+        distances=None,
+        betas=None,
+        compute_closeness=None,
+        compute_betweenness=None,
+        min_threshold_wt=None,
+        jitter_scale=None,
+        pbar_disabled=None
+    ))]
     pub fn local_node_centrality_shortest(
         &self,
         distances: Option<Vec<u32>>,
@@ -532,7 +544,17 @@ impl NetworkStructure {
         });
         Ok(result)
     }
-
+    #[pyo3(signature = (
+        distances=None,
+        betas=None,
+        compute_closeness=None,
+        compute_betweenness=None,
+        min_threshold_wt=None,
+        angular_scaling_unit=None,
+        farness_scaling_offset=None,
+        jitter_scale=None,
+        pbar_disabled=None
+    ))]
     pub fn local_node_centrality_simplest(
         &self,
         distances: Option<Vec<u32>>,
@@ -658,7 +680,15 @@ impl NetworkStructure {
         });
         Ok(result)
     }
-
+    #[pyo3(signature = (
+        distances=None,
+        betas=None,
+        compute_closeness=None,
+        compute_betweenness=None,
+        min_threshold_wt=None,
+        jitter_scale=None,
+        pbar_disabled=None
+    ))]
     pub fn local_segment_centrality(
         &self,
         distances: Option<Vec<u32>>,
diff --git a/src/common.rs b/src/common.rs
index 4e5ba6c6..8bdb2ff9 100644
--- a/src/common.rs
+++ b/src/common.rs
@@ -7,8 +7,10 @@ use std::collections::HashMap;
 use std::f32::consts::PI;
 use std::sync::atomic::Ordering;
 
+/// Minimum threshold weight for distance and beta calculations.
 static MIN_THRESH_WT: f32 = 0.01831563888873418;
 
+/// Represents a 2D coordinate with `x` and `y` values.
 #[pyclass]
 #[derive(Clone, Copy)]
 pub struct Coord {
@@ -17,59 +19,99 @@ pub struct Coord {
     #[pyo3(get)]
     pub y: f32,
 }
+
 #[pymethods]
 impl Coord {
+    /// Creates a new `Coord` instance.
     #[new]
     pub fn new(x: f32, y: f32) -> Self {
         Self { x, y }
     }
+
+    /// Returns the coordinates as a tuple `(x, y)`.
     pub fn xy(&self) -> (f32, f32) {
         (self.x, self.y)
     }
+
+    /// Validates that the coordinates are finite.
     pub fn validate(&self) -> bool {
         self.x.is_finite() && self.y.is_finite()
     }
+
+    /// Calculates the Euclidean distance between this coordinate and another.
     pub fn hypot(&self, other_coord: Coord) -> f32 {
         ((self.x - other_coord.x).powi(2) + (self.y - other_coord.y).powi(2)).sqrt()
     }
+
+    /// Computes the difference between this coordinate and another as a vector.
     pub fn difference(&self, other_coord: Coord) -> Coord {
-        // using Coord struct as a vector
-        let x_diff = self.x - other_coord.x;
-        let y_diff = self.y - other_coord.y;
-        Coord::new(x_diff, y_diff)
+        Coord::new(self.x - other_coord.x, self.y - other_coord.y)
     }
 }
+
+/// Holds metric results, including distances and a 2D matrix of atomic floats.
 pub struct MetricResult {
     pub distances: Vec<u32>,
     pub metric: Vec<Vec<AtomicF32>>,
 }
+
 impl MetricResult {
+    /// Initializes a new `MetricResult` with given distances, size, and initial value.
     pub fn new(distances: Vec<u32>, size: usize, init_val: f32) -> Self {
-        let mut metric = Vec::new();
-        for _d in 0..distances.len() {
-            metric.push(
-                // tricky to initialise for given size
+        let metric = distances
+            .iter()
+            .map(|_| {
                 std::iter::repeat_with(|| AtomicF32::new(init_val))
                     .take(size)
-                    .collect::<Vec<AtomicF32>>(),
-            );
-        }
+                    .collect::<Vec<AtomicF32>>()
+            })
+            .collect();
         Self { distances, metric }
     }
+
+    /// Converts the atomic floats into a Python-compatible format (`PyArray1<f32>`).
     pub fn load(&self) -> HashMap<u32, Py<PyArray1<f32>>> {
-        let mut loaded: HashMap<u32, Py<PyArray1<f32>>> = HashMap::new();
-        for i in 0..self.distances.len() {
-            let dist = self.distances[i];
-            let vec_f32: Vec<f32> = self.metric[i]
+        self.distances
+            .iter()
+            .enumerate()
+            .map(|(i, &dist)| {
+                let vec_f32: Vec<f32> = self.metric[i]
+                    .iter()
+                    .map(|a| a.load(Ordering::SeqCst))
+                    .collect();
+                let array = Python::with_gil(|py| {
+                    vec_f32
+                        .into_pyarray(py)
+                        .to_owned() // This gives us a PyArray, but wrapped in pyo3::Bound
+                        .into() // Convert to the required Py<PyArray> type
+                });
+                (dist, array)
+            })
+            .collect()
+    }
+}
+// Manually implement Clone for MetricResult
+impl Clone for MetricResult {
+    fn clone(&self) -> Self {
+        MetricResult {
+            distances: self.distances.clone(), // Clone the distances (Vec<u32>)
+            metric: self
+                .metric
                 .iter()
-                .map(|a| a.load(Ordering::SeqCst))
-                .collect();
-            let array = Python::with_gil(|py| vec_f32.into_pyarray(py).to_owned());
-            loaded.insert(dist, array);
+                .map(|row| {
+                    row.iter()
+                        .map(|atomic_f32| {
+                            // Here we clone by reading the value and then creating a new AtomicF32
+                            AtomicF32::new(atomic_f32.load(Ordering::SeqCst))
+                        })
+                        .collect()
+                })
+                .collect(),
         }
-        loaded
     }
 }
+
+/// Calculates the rotation angle between two points relative to the origin.
 #[pyfunction]
 pub fn calculate_rotation(point_a: Coord, point_b: Coord) -> f32 {
     let ang_a = point_a.y.atan2(point_a.x);
@@ -77,25 +119,23 @@ pub fn calculate_rotation(point_a: Coord, point_b: Coord) -> f32 {
     let rotation = (ang_a - ang_b) % (2.0 * PI);
     rotation.to_degrees()
 }
-// https://stackoverflow.com/questions/37459121/calculating-angle-between-three-points-but-only-anticlockwise-in-python
-// these two points / angles are relative to the origin
-// pass in difference between the points and origin as vectors
+
+/// Calculates the smallest difference angle between two vectors.
 #[pyfunction]
 pub fn calculate_rotation_smallest(vec_a: Coord, vec_b: Coord) -> f32 {
-    // Convert angles from radians to degrees and calculate the smallest difference angle
-    let ang_a = (vec_a.y.atan2(vec_a.x)).to_degrees();
-    let ang_b = (vec_b.y.atan2(vec_b.x)).to_degrees();
+    let ang_a = vec_a.y.atan2(vec_a.x).to_degrees();
+    let ang_b = vec_b.y.atan2(vec_b.x).to_degrees();
     let diff_angle = (ang_b - ang_a + 180.0) % 360.0 - 180.0;
     diff_angle.abs()
 }
+
+/// Validates that all elements in a 2D NumPy array are finite.
 #[pyfunction]
 pub fn check_numerical_data(data_arr: PyReadonlyArray2<f32>) -> PyResult<()> {
-    // Check the integrity of numeric data arrays.
     let data_slice = data_arr.as_array();
     for inner_arr in data_slice.rows() {
-        for num in inner_arr.iter() {
-            let num_val = *num;
-            if !num_val.is_finite() {
+        for &num in inner_arr.iter() {
+            if !num.is_finite() {
                 return Err(exceptions::PyValueError::new_err(
                     "The numeric data values must be finite.",
                 ));
@@ -105,9 +145,11 @@ pub fn check_numerical_data(data_arr: PyReadonlyArray2<f32>) -> PyResult<()> {
     Ok(())
 }
 
+/// Converts beta values to distances using a logarithmic transformation.
 #[pyfunction]
+#[pyo3(signature = (betas, min_threshold_wt=None))]
 pub fn distances_from_betas(betas: Vec<f32>, min_threshold_wt: Option<f32>) -> PyResult<Vec<u32>> {
-    if betas.len() == 0 {
+    if betas.is_empty() {
         return Err(exceptions::PyValueError::new_err(
             "Empty iterable of betas.",
         ));
@@ -115,34 +157,37 @@ pub fn distances_from_betas(betas: Vec<f32>, min_threshold_wt: Option<f32>) -> P
     let min_threshold_wt = min_threshold_wt.unwrap_or(MIN_THRESH_WT);
     let mut clean: Vec<f32> = Vec::new();
     let mut distances: Vec<u32> = Vec::new();
-    for beta in betas.iter() {
-        if *beta < 0.0 {
+
+    for &beta in &betas {
+        if beta < 0.0 {
             return Err(exceptions::PyValueError::new_err(
                 "Provide the beta value without the leading negative.",
             ));
         }
-        if *beta == 0.0 {
+        if beta == 0.0 {
             return Err(exceptions::PyValueError::new_err(
                 "Provide a beta value greater than zero.",
             ));
         }
-        if clean.contains(beta) || clean.iter().any(|&x| x < *beta) {
+        if clean.contains(&beta) || clean.iter().any(|&x| x < beta) {
             return Err(exceptions::PyValueError::new_err(
                 "Betas must be free of duplicates and sorted in decreasing order.",
             ));
         }
-        clean.push(*beta);
+        clean.push(beta);
         distances.push((min_threshold_wt.ln() / -beta).round() as u32);
     }
     Ok(distances)
 }
 
+/// Converts distances back to beta values.
 #[pyfunction]
+#[pyo3(signature = (distances, min_threshold_wt=None))]
 pub fn betas_from_distances(
     distances: Vec<u32>,
     min_threshold_wt: Option<f32>,
 ) -> PyResult<Vec<f32>> {
-    if distances.len() == 0 {
+    if distances.is_empty() {
         return Err(exceptions::PyValueError::new_err(
             "Empty iterable of distances.",
         ));
@@ -150,112 +195,111 @@ pub fn betas_from_distances(
     let min_threshold_wt = min_threshold_wt.unwrap_or(MIN_THRESH_WT);
     let mut clean: Vec<u32> = Vec::new();
     let mut betas: Vec<f32> = Vec::new();
-    for distance in distances.iter() {
-        if *distance <= 0 {
+
+    for &distance in &distances {
+        if distance == 0 {
             return Err(exceptions::PyValueError::new_err(
                 "Distances must be positive integers.",
             ));
         }
-        if clean.contains(distance) || clean.iter().any(|&x| x > *distance) {
+        if clean.contains(&distance) || clean.iter().any(|&x| x > distance) {
             return Err(exceptions::PyValueError::new_err(
                 "Distances must be free of duplicates and sorted in increasing order.",
             ));
         }
-        clean.push(*distance);
-        betas.push(-min_threshold_wt.ln() / *distance as f32);
+        clean.push(distance);
+        betas.push(-min_threshold_wt.ln() / distance as f32);
     }
     Ok(betas)
 }
 
+/// Pairs distances and betas, ensuring consistency.
 #[pyfunction]
+#[pyo3(signature = (distances=None, betas=None, min_threshold_wt=None))]
 pub fn pair_distances_and_betas(
     distances: Option<Vec<u32>>,
     betas: Option<Vec<f32>>,
     min_threshold_wt: Option<f32>,
 ) -> PyResult<(Vec<u32>, Vec<f32>)> {
-    if distances.is_some() && betas.is_some() {
-        return Err(exceptions::PyValueError::new_err(
-            "Please provide either a distances or betas, not both.",
-        ));
-    }
-    if distances.is_none() && betas.is_none() {
-        return Err(exceptions::PyValueError::new_err(
-            "Please provide either a distances or betas. Neither has been provided",
-        ));
+    match (distances, betas) {
+        (Some(_), Some(_)) => Err(exceptions::PyValueError::new_err(
+            "Please provide either distances or betas, not both.",
+        )),
+        (None, None) => Err(exceptions::PyValueError::new_err(
+            "Please provide either distances or betas. Neither has been provided.",
+        )),
+        (Some(distances), None) => {
+            let betas = betas_from_distances(distances.clone(), min_threshold_wt)?;
+            Ok((distances, betas))
+        }
+        (None, Some(betas)) => {
+            let distances = distances_from_betas(betas.clone(), min_threshold_wt)?;
+            Ok((distances, betas))
+        }
     }
-    let betas = if betas.is_some() {
-        betas.unwrap()
-    } else {
-        betas_from_distances(distances.clone().unwrap(), min_threshold_wt)?
-    };
-    let distances = if distances.is_some() {
-        distances.unwrap()
-    } else {
-        distances_from_betas(betas.clone(), min_threshold_wt)?
-    };
-    Ok((distances, betas))
 }
 
+/// Computes average distances for given beta values.
 #[pyfunction]
+#[pyo3(signature = (betas, min_threshold_wt=None))]
 pub fn avg_distances_for_betas(
     betas: Vec<f32>,
     min_threshold_wt: Option<f32>,
 ) -> PyResult<Vec<f32>> {
-    if betas.len() == 0 {
+    if betas.is_empty() {
         return Err(exceptions::PyValueError::new_err(
             "Empty iterable of betas.",
         ));
     }
     let min_threshold_wt = min_threshold_wt.unwrap_or(MIN_THRESH_WT);
-    let mut avg_distances: Vec<f32> = Vec::new();
     let distances = distances_from_betas(betas.clone(), Some(min_threshold_wt))?;
-    for (beta, distance) in betas.iter().zip(distances.iter()) {
-        if *distance <= 0 {
-            return Err(exceptions::PyValueError::new_err(
-                "Distances must be positive integers.",
-            ));
-        }
-        let auc = ((-beta * *distance as f32).exp() - 1.0) / -beta;
-        let wt = auc / *distance as f32;
-        let avg_d = -wt.ln() / beta;
-        avg_distances.push(avg_d)
-    }
+    let avg_distances: Vec<f32> = betas
+        .iter()
+        .zip(distances.iter())
+        .map(|(&beta, &distance)| {
+            if distance == 0 {
+                return Err(exceptions::PyValueError::new_err(
+                    "Distances must be positive integers.",
+                ));
+            }
+            let auc = ((-beta * distance as f32).exp() - 1.0) / -beta;
+            let wt = auc / distance as f32;
+            Ok(-wt.ln() / beta)
+        })
+        .collect::<PyResult<_>>()?;
     Ok(avg_distances)
 }
 
+/// Clips weights based on a spatial tolerance.
 #[pyfunction]
 pub fn clip_wts_curve(
     distances: Vec<u32>,
     betas: Vec<f32>,
     spatial_tolerance: u32,
 ) -> PyResult<Vec<f32>> {
-    let mut max_curve_wts: Vec<f32> = Vec::new();
-    for (dist, beta) in distances.iter().zip(betas.iter()) {
-        if spatial_tolerance > *dist {
-            return Err(exceptions::PyValueError::new_err(
-                "Clipping distance cannot be greater than the given distance threshold.",
-            ));
-        }
-        let max_curve_wt = (-beta * spatial_tolerance as f32).exp();
-        if max_curve_wt < 0.75 {}
-        max_curve_wts.push(max_curve_wt);
-    }
+    let max_curve_wts: Vec<f32> = distances
+        .iter()
+        .zip(betas.iter())
+        .map(|(&dist, &beta)| {
+            if spatial_tolerance > dist {
+                return Err(exceptions::PyValueError::new_err(
+                    "Clipping distance cannot be greater than the given distance threshold.",
+                ));
+            }
+            Ok((-beta * spatial_tolerance as f32).exp())
+        })
+        .collect::<PyResult<_>>()?;
     Ok(max_curve_wts)
 }
 
+/// Computes a clipped weight based on a beta value and maximum curve weight.
 #[pyfunction]
 pub fn clipped_beta_wt(beta: f32, max_curve_wt: f32, data_dist: f32) -> PyResult<f32> {
-    if beta < 0.0 || beta > 1.0 {
-        return Err(exceptions::PyValueError::new_err(
-            "Max curve weight must be in a range of 0 - 1.",
-        ));
-    }
-    if max_curve_wt < 0.0 || max_curve_wt > 1.0 {
+    if !(0.0..=1.0).contains(&max_curve_wt) {
         return Err(exceptions::PyValueError::new_err(
-            "Max curve weight must be in a range of 0 - 1.",
+            "Max curve weight must be in the range [0, 1].",
         ));
     }
-    // Calculates negative exponential clipped to the max_curve_wt parameter.
     let raw_wt = (-beta * data_dist).exp();
     let clipped_wt = f32::min(raw_wt, max_curve_wt) / max_curve_wt;
     Ok(clipped_wt)
diff --git a/src/data.rs b/src/data.rs
index d86fc92e..c695441f 100644
--- a/src/data.rs
+++ b/src/data.rs
@@ -74,6 +74,7 @@ pub struct DataEntry {
 #[pymethods]
 impl DataEntry {
     #[new]
+    #[pyo3(signature = (data_key, x, y, data_id=None, nearest_assign=None, next_nearest_assign=None))]
     fn new(
         data_key: String,
         x: f32,
@@ -116,6 +117,7 @@ impl DataMap {
     fn progress(&self) -> usize {
         self.progress.as_ref().load(Ordering::Relaxed)
     }
+    #[pyo3(signature = (data_key, x, y, data_id=None, nearest_assign=None, next_nearest_assign=None))]
     fn insert(
         &mut self,
         data_key: String,
@@ -179,6 +181,13 @@ impl DataMap {
             entry.next_nearest_assign = Some(assign_idx);
         }
     }
+    #[pyo3(signature = (
+        netw_src_idx,
+        network_structure,
+        max_dist,
+        jitter_scale=None,
+        angular=None
+    ))]
     fn aggregate_to_src_idx(
         &self,
         netw_src_idx: usize,
@@ -279,6 +288,18 @@ impl DataMap {
         }
         entries
     }
+    #[pyo3(signature = (
+        network_structure,
+        landuses_map,
+        accessibility_keys,
+        distances=None,
+        betas=None,
+        angular=None,
+        spatial_tolerance=None,
+        min_threshold_wt=None,
+        jitter_scale=None,
+        pbar_disabled=None
+    ))]
     fn accessibility(
         &self,
         network_structure: &NetworkStructure,
@@ -421,7 +442,21 @@ impl DataMap {
         });
         Ok(result)
     }
-
+    #[pyo3(signature = (
+        network_structure,
+        landuses_map,
+        distances=None,
+        betas=None,
+        compute_hill=None,
+        compute_hill_weighted=None,
+        compute_shannon=None,
+        compute_gini=None,
+        angular=None,
+        spatial_tolerance=None,
+        min_threshold_wt=None,
+        jitter_scale=None,
+        pbar_disabled=None
+    ))]
     fn mixed_uses(
         &self,
         network_structure: &NetworkStructure,
@@ -665,10 +700,21 @@ impl DataMap {
         });
         Ok(result)
     }
+    #[pyo3(signature = (
+        network_structure,
+        numerical_maps,
+        distances=None,
+        betas=None,
+        angular=None,
+        spatial_tolerance=None,
+        min_threshold_wt=None,
+        jitter_scale=None,
+        pbar_disabled=None
+    ))]
     fn stats(
         &self,
         network_structure: &NetworkStructure,
-        numerical_map: HashMap<String, f32>,
+        numerical_maps: Vec<HashMap<String, f32>>, // vector of numerical maps
         distances: Option<Vec<u32>>,
         betas: Option<Vec<f32>>,
         angular: Option<bool>,
@@ -677,44 +723,94 @@ impl DataMap {
         jitter_scale: Option<f32>,
         pbar_disabled: Option<bool>,
         py: Python,
-    ) -> PyResult<StatsResult> {
+    ) -> PyResult<Vec<StatsResult>> {
+        // Return a vector of StatsResult
         let (distances, betas) = pair_distances_and_betas(distances, betas, min_threshold_wt)?;
         let max_dist: u32 = distances.iter().max().unwrap().clone();
-        if numerical_map.len() != self.count() {
-            return Err(exceptions::PyValueError::new_err(
-                "The number of numerical entries must match the number of data points",
-            ));
+        // Iterate through each map in the numerical_maps
+        for (index, numerical_map) in numerical_maps.iter().enumerate() {
+            if numerical_map.len() != self.count() {
+                return Err(exceptions::PyValueError::new_err(
+                    format!(
+                        "The number of entries in numerical map {} must match the number of data points (expected: {}, found: {})",
+                        index,
+                        self.count(),
+                        numerical_map.len()
+                    ),
+                ));
+            }
         }
         let spatial_tolerance = spatial_tolerance.unwrap_or(0);
         let max_curve_wts = clip_wts_curve(distances.clone(), betas.clone(), spatial_tolerance)?;
         // track progress
         let pbar_disabled = pbar_disabled.unwrap_or(false);
         self.progress_init();
+        // collect results for all metrics
         let result = py.allow_threads(move || {
-            // prepare the containers for tracking results
-            let sum = MetricResult::new(distances.clone(), network_structure.node_count(), 0.0);
-            let sum_wt = MetricResult::new(distances.clone(), network_structure.node_count(), 0.0);
-            let count = MetricResult::new(distances.clone(), network_structure.node_count(), 0.0);
-            let count_wt =
-                MetricResult::new(distances.clone(), network_structure.node_count(), 0.0);
-            let max = MetricResult::new(
-                distances.clone(),
-                network_structure.node_count(),
-                f32::NEG_INFINITY,
-            );
-            let min = MetricResult::new(
-                distances.clone(),
-                network_structure.node_count(),
-                f32::INFINITY,
-            );
-            let mean =
-                MetricResult::new(distances.clone(), network_structure.node_count(), f32::NAN);
-            let mean_wt =
-                MetricResult::new(distances.clone(), network_structure.node_count(), f32::NAN);
-            let variance =
-                MetricResult::new(distances.clone(), network_structure.node_count(), f32::NAN);
-            let variance_wt =
-                MetricResult::new(distances.clone(), network_structure.node_count(), f32::NAN);
+            // initialize vectors to hold the metrics for each statistical result
+            let mut sum = Vec::new();
+            let mut sum_wt = Vec::new();
+            let mut count = Vec::new();
+            let mut count_wt = Vec::new();
+            let mut max = Vec::new();
+            let mut min = Vec::new();
+            let mut mean = Vec::new();
+            let mut mean_wt = Vec::new();
+            let mut variance = Vec::new();
+            let mut variance_wt = Vec::new();
+            // initialize each metric's result for all distances and node counts
+            for _ in 0..numerical_maps.len() {
+                sum.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    0.0,
+                ));
+                sum_wt.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    0.0,
+                ));
+                count.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    0.0,
+                ));
+                count_wt.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    0.0,
+                ));
+                max.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    f32::NEG_INFINITY,
+                ));
+                min.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    f32::INFINITY,
+                ));
+                mean.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    f32::NAN,
+                ));
+                mean_wt.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    f32::NAN,
+                ));
+                variance.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    f32::NAN,
+                ));
+                variance_wt.push(MetricResult::new(
+                    distances.clone(),
+                    network_structure.node_count(),
+                    f32::NAN,
+                ));
+            }
             // indices
             let node_indices: Vec<usize> = network_structure.node_indices();
             // iter
@@ -742,126 +838,163 @@ impl DataMap {
                 sort by increasing distance re: deduplication via data keys
                 because these are sorted, no need to deduplicate by respective distance thresholds
                 */
+                // iterate over each numerical map (each map contains numerical data for different metrics)
                 for (data_key, data_dist) in reachable_entries.iter() {
-                    let num = numerical_map[data_key].clone();
-                    if num.is_nan() {
-                        continue;
-                    };
-                    for i in 0..distances.len() {
-                        let d = distances[i];
-                        let b = betas[i];
-                        let mcw = max_curve_wts[i];
-                        if *data_dist <= d as f32 {
-                            let wt = clipped_beta_wt(b, mcw, *data_dist).unwrap();
-                            let num_wt = num * wt;
-                            // agg
-                            sum.metric[i][*netw_src_idx].fetch_add(num, Ordering::Relaxed);
-                            sum_wt.metric[i][*netw_src_idx].fetch_add(num_wt, Ordering::Relaxed);
-                            count.metric[i][*netw_src_idx].fetch_add(1.0, Ordering::Relaxed);
-                            count_wt.metric[i][*netw_src_idx].fetch_add(wt, Ordering::Relaxed);
-                            // not using compare_exchang in loop because only one netw_src_idx is processed at a time
-                            let current_max = max.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                            if num > current_max {
-                                max.metric[i][*netw_src_idx].store(num, Ordering::Relaxed);
-                            };
-                            let current_min = min.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                            if num < current_min {
-                                min.metric[i][*netw_src_idx].store(num, Ordering::Relaxed);
-                            };
+                    for (map_idx, numerical_map) in numerical_maps.iter().enumerate() {
+                        let num = numerical_map.get(data_key).cloned().unwrap_or(f32::NAN);
+                        if num.is_nan() {
+                            continue;
+                        }
+                        for i in 0..distances.len() {
+                            let d = distances[i];
+                            let b = betas[i];
+                            let mcw = max_curve_wts[i];
+                            if *data_dist <= d as f32 {
+                                let wt = clipped_beta_wt(b, mcw, *data_dist).unwrap();
+                                let num_wt = num * wt;
+                                // update metrics for this particular map
+                                sum[map_idx].metric[i][*netw_src_idx]
+                                    .fetch_add(num, Ordering::Relaxed);
+                                sum_wt[map_idx].metric[i][*netw_src_idx]
+                                    .fetch_add(num_wt, Ordering::Relaxed);
+                                count[map_idx].metric[i][*netw_src_idx]
+                                    .fetch_add(1.0, Ordering::Relaxed);
+                                count_wt[map_idx].metric[i][*netw_src_idx]
+                                    .fetch_add(wt, Ordering::Relaxed);
+                                // update max/min values
+                                let current_max =
+                                    max[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                                if num > current_max {
+                                    max[map_idx].metric[i][*netw_src_idx]
+                                        .store(num, Ordering::Relaxed);
+                                };
+                                let current_min =
+                                    min[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                                if num < current_min {
+                                    min[map_idx].metric[i][*netw_src_idx]
+                                        .store(num, Ordering::Relaxed);
+                                };
+                            }
                         }
                     }
                 }
-                // finalise mean calculations - this is happening for a single netw_src_idx, so fairly fast
-                for i in 0..distances.len() {
-                    let sum_val = sum.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    let count_val = count.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    mean.metric[i][*netw_src_idx].store(sum_val / count_val, Ordering::Relaxed);
-                    let sum_wt_val = sum_wt.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    let count_wt_val = count_wt.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    mean_wt.metric[i][*netw_src_idx]
-                        .store(sum_wt_val / count_wt_val, Ordering::Relaxed);
-                    // also clean up min and max - e.g. isolated locations
-                    let current_max_val = max.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    if current_max_val.is_infinite() {
-                        max.metric[i][*netw_src_idx].store(f32::NAN, Ordering::Relaxed);
-                    }
-                    let current_min_val = min.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    if current_min_val.is_infinite() {
-                        min.metric[i][*netw_src_idx].store(f32::NAN, Ordering::Relaxed);
-                    }
-                }
-                // calculate variances - counts are already computed per above
-                // weighted version is IDW by division through equivalently weighted counts above
-                for (data_key, data_dist) in reachable_entries {
-                    let num = numerical_map[&data_key].clone();
-                    if num.is_nan() {
-                        continue;
-                    };
+                // finalise mean and variance calculations for each map
+                for map_idx in 0..numerical_maps.len() {
                     for i in 0..distances.len() {
-                        let d = distances[i];
-                        if data_dist <= d as f32 {
-                            let current_var_val =
-                                variance.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                            let current_mean_val =
-                                mean.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                            let diff = num - current_mean_val;
-                            if current_var_val.is_nan() {
-                                variance.metric[i][*netw_src_idx]
-                                    .store(diff * diff, Ordering::Relaxed);
-                            } else {
-                                variance.metric[i][*netw_src_idx]
-                                    .fetch_add(diff * diff, Ordering::Relaxed);
-                            }
-                            let current_var_wt_val =
-                                variance_wt.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                            let current_mean_wt_val =
-                                mean_wt.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                            let diff_wt = num - current_mean_wt_val;
-                            if current_var_wt_val.is_nan() {
-                                variance_wt.metric[i][*netw_src_idx]
-                                    .store(diff_wt * diff_wt, Ordering::Relaxed);
-                            } else {
-                                variance_wt.metric[i][*netw_src_idx]
-                                    .fetch_add(diff_wt * diff_wt, Ordering::Relaxed);
+                        let sum_val = sum[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        let count_val =
+                            count[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        mean[map_idx].metric[i][*netw_src_idx]
+                            .store(sum_val / count_val, Ordering::Relaxed);
+                        let sum_wt_val =
+                            sum_wt[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        let count_wt_val =
+                            count_wt[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        mean_wt[map_idx].metric[i][*netw_src_idx]
+                            .store(sum_wt_val / count_wt_val, Ordering::Relaxed);
+                        // clean up min and max values
+                        let current_max_val =
+                            max[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        if current_max_val.is_infinite() {
+                            max[map_idx].metric[i][*netw_src_idx]
+                                .store(f32::NAN, Ordering::Relaxed);
+                        }
+                        let current_min_val =
+                            min[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        if current_min_val.is_infinite() {
+                            min[map_idx].metric[i][*netw_src_idx]
+                                .store(f32::NAN, Ordering::Relaxed);
+                        }
+                    }
+                    // calculate variances
+                    for (data_key, data_dist) in reachable_entries.iter() {
+                        let num = numerical_maps[map_idx]
+                            .get(data_key)
+                            .cloned()
+                            .unwrap_or(f32::NAN);
+                        if num.is_nan() {
+                            continue;
+                        };
+                        for i in 0..distances.len() {
+                            let d = distances[i];
+                            if *data_dist <= d as f32 {
+                                let current_var_val = variance[map_idx].metric[i][*netw_src_idx]
+                                    .load(Ordering::Relaxed);
+                                let current_mean_val =
+                                    mean[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                                let diff = num - current_mean_val;
+                                if current_var_val.is_nan() {
+                                    variance[map_idx].metric[i][*netw_src_idx]
+                                        .store(diff * diff, Ordering::Relaxed);
+                                } else {
+                                    variance[map_idx].metric[i][*netw_src_idx]
+                                        .fetch_add(diff * diff, Ordering::Relaxed);
+                                }
+                                let current_var_wt_val = variance_wt[map_idx].metric[i]
+                                    [*netw_src_idx]
+                                    .load(Ordering::Relaxed);
+                                let current_mean_wt_val = mean_wt[map_idx].metric[i][*netw_src_idx]
+                                    .load(Ordering::Relaxed);
+                                let diff_wt = num - current_mean_wt_val;
+                                if current_var_wt_val.is_nan() {
+                                    variance_wt[map_idx].metric[i][*netw_src_idx]
+                                        .store(diff_wt * diff_wt, Ordering::Relaxed);
+                                } else {
+                                    variance_wt[map_idx].metric[i][*netw_src_idx]
+                                        .fetch_add(diff_wt * diff_wt, Ordering::Relaxed);
+                                }
                             }
                         }
                     }
                 }
                 // finalise variance calculations
-                for i in 0..distances.len() {
-                    let current_var_val = variance.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    let current_count_val = count.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    if current_count_val == 0.0 {
-                        variance.metric[i][*netw_src_idx].store(f32::NAN, Ordering::Relaxed);
-                    } else {
-                        variance.metric[i][*netw_src_idx]
-                            .store(current_var_val / current_count_val, Ordering::Relaxed);
-                    }
-                    let current_var_wt_val =
-                        variance_wt.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    let current_count_wt_val =
-                        count_wt.metric[i][*netw_src_idx].load(Ordering::Relaxed);
-                    if current_count_wt_val == 0.0 {
-                        variance_wt.metric[i][*netw_src_idx].store(f32::NAN, Ordering::Relaxed);
-                    } else {
-                        variance_wt.metric[i][*netw_src_idx]
-                            .store(current_var_wt_val / current_count_wt_val, Ordering::Relaxed);
+                for map_idx in 0..numerical_maps.len() {
+                    for i in 0..distances.len() {
+                        let current_var_val =
+                            variance[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        let current_count_val =
+                            count[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        if current_count_val == 0.0 {
+                            variance[map_idx].metric[i][*netw_src_idx]
+                                .store(f32::NAN, Ordering::Relaxed);
+                        } else {
+                            variance[map_idx].metric[i][*netw_src_idx]
+                                .store(current_var_val / current_count_val, Ordering::Relaxed);
+                        }
+                        let current_var_wt_val =
+                            variance_wt[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        let current_count_wt_val =
+                            count_wt[map_idx].metric[i][*netw_src_idx].load(Ordering::Relaxed);
+                        if current_count_wt_val == 0.0 {
+                            variance_wt[map_idx].metric[i][*netw_src_idx]
+                                .store(f32::NAN, Ordering::Relaxed);
+                        } else {
+                            variance_wt[map_idx].metric[i][*netw_src_idx].store(
+                                current_var_wt_val / current_count_wt_val,
+                                Ordering::Relaxed,
+                            );
+                        }
                     }
                 }
             });
-            // unpack
-            StatsResult {
-                sum: sum.load(),
-                sum_wt: sum_wt.load(),
-                mean: mean.load(),
-                mean_wt: mean_wt.load(),
-                count: count.load(),
-                count_wt: count_wt.load(),
-                variance: variance.load(),
-                variance_wt: variance_wt.load(),
-                max: max.load(),
-                min: min.load(),
+
+            // create the results vector from the updated metrics
+            let mut results = Vec::new();
+            for map_idx in 0..numerical_maps.len() {
+                results.push(StatsResult {
+                    sum: sum[map_idx].load(),
+                    sum_wt: sum_wt[map_idx].load(),
+                    mean: mean[map_idx].load(),
+                    mean_wt: mean_wt[map_idx].load(),
+                    count: count[map_idx].load(),
+                    count_wt: count_wt[map_idx].load(),
+                    variance: variance[map_idx].load(),
+                    variance_wt: variance_wt[map_idx].load(),
+                    max: max[map_idx].load(),
+                    min: min[map_idx].load(),
+                });
             }
+            results
         });
         Ok(result)
     }
diff --git a/src/lib.rs b/src/lib.rs
index 86ae2a7b..2d8dea93 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -7,9 +7,7 @@ mod graph;
 mod viewshed;
 
 #[pymodule]
-fn rustalgos(_py: Python, m: &PyModule) -> PyResult<()> {
-    // let rustalgos = PyModule::new(py, "rustalgos")?;
-    // m.add_submodule(rustalgos)?;
+fn rustalgos(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<common::Coord>()?;
     m.add_function(wrap_pyfunction!(common::calculate_rotation, m)?)?;
     m.add_function(wrap_pyfunction!(common::calculate_rotation_smallest, m)?)?;
@@ -47,7 +45,5 @@ fn rustalgos(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_class::<centrality::CentralitySegmentResult>()?;
     // VIEWSHED
     m.add_class::<viewshed::Viewshed>()?;
-    // let sys = PyModule::import(py, "sys")?;
-    // sys.getattr("modules")?.set_item("rustalgos", rustalgos)?;
     Ok(())
 }
diff --git a/src/viewshed.rs b/src/viewshed.rs
index bb805603..bdb16457 100644
--- a/src/viewshed.rs
+++ b/src/viewshed.rs
@@ -9,6 +9,7 @@ use std::sync::Arc;
 pub struct Viewshed {
     pub progress: Arc<AtomicUsize>,
 }
+
 fn line_of_sight(
     raster: ArrayView2<u8>,
     start_x: usize,
@@ -49,6 +50,7 @@ fn line_of_sight(
     }
     true
 }
+
 fn calculate_visible_cells(
     raster: ArrayView2<u8>,
     start_x: usize,
@@ -86,6 +88,7 @@ fn calculate_visible_cells(
 
     (density, farness, harmonic)
 }
+
 fn calculate_viewshed(
     raster: ArrayView2<u8>,
     start_x: usize,
@@ -132,6 +135,7 @@ impl Viewshed {
     fn progress(&self) -> usize {
         self.progress.as_ref().load(Ordering::Relaxed)
     }
+    #[pyo3(signature = (bldgs_rast, view_distance, pbar_disabled=None))]
     pub fn visibility_graph(
         &self,
         bldgs_rast: PyReadonlyArray2<u8>,
@@ -186,8 +190,9 @@ impl Viewshed {
             .into_pyarray(py)
             .to_owned();
 
-        Ok((array_u32, array_f32_a, array_f32_b))
+        Ok((array_u32.into(), array_f32_a.into(), array_f32_b.into()))
     }
+
     pub fn viewshed(
         &self,
         bldgs_rast: PyReadonlyArray2<u8>,
@@ -203,6 +208,6 @@ impl Viewshed {
             .unwrap()
             .into_pyarray(py)
             .to_owned();
-        Ok(numpy_array)
+        Ok(numpy_array.into())
     }
 }
diff --git a/tests/metrics/test_layers.py b/tests/metrics/test_layers.py
index 9cf27631..db325078 100644
--- a/tests/metrics/test_layers.py
+++ b/tests/metrics/test_layers.py
@@ -182,7 +182,7 @@ def test_compute_stats(primal_graph):
     Test stats component
     """
     nodes_gdf, _edges_gdf, network_structure = io.network_structure_from_nx(primal_graph, 3395)
-    data_gdf = mock.mock_numerical_data(primal_graph, num_arrs=1)
+    data_gdf = mock.mock_numerical_data(primal_graph, num_arrs=2)
     max_assign_dist = 400
     data_map, data_gdf = layers.assign_gdf_to_network(data_gdf, network_structure, max_assign_dist)
     # test against manual implementation over underlying method
@@ -191,98 +191,99 @@ def test_compute_stats(primal_graph):
         for angular in [False, True]:
             nodes_gdf, data_gdf = layers.compute_stats(
                 data_gdf,
-                "mock_numerical_1",
+                ["mock_numerical_1", "mock_numerical_2"],
                 nodes_gdf,
                 network_structure,
                 max_assign_dist,
                 distances=distances,
                 angular=angular,
             )
-            # generate stats
             # compare to manual
-            numerical_map = data_gdf["mock_numerical_1"].to_dict()
-            stats_result = data_map.stats(
-                network_structure,
-                numerical_map=numerical_map,
-                distances=distances,
-                angular=angular,
-            )
-            for dist_key in distances:
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_sum", dist_key, angular=angular, weighted=False)],
-                    stats_result.sum[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_sum", dist_key, angular=angular, weighted=True)],
-                    stats_result.sum_wt[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_mean", dist_key, angular=angular, weighted=False)],
-                    stats_result.mean[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_mean", dist_key, angular=angular, weighted=True)],
-                    stats_result.mean_wt[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_count", dist_key, angular=angular, weighted=False)],
-                    stats_result.count[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_count", dist_key, angular=angular, weighted=True)],
-                    stats_result.count_wt[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_var", dist_key, angular=angular, weighted=False)],
-                    stats_result.variance[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_var", dist_key, angular=angular, weighted=True)],
-                    stats_result.variance_wt[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_max", dist_key, angular=angular)],
-                    stats_result.max[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
-                )
-                assert np.allclose(
-                    nodes_gdf[config.prep_gdf_key("mock_numerical_1_min", dist_key, angular=angular)],
-                    stats_result.min[dist_key],
-                    atol=config.ATOL,
-                    rtol=config.RTOL,
-                    equal_nan=True,
+            for stats_key in ["mock_numerical_1", "mock_numerical_2"]:
+                # generate stats
+                stats_results = data_map.stats(
+                    network_structure,
+                    numerical_maps=[data_gdf[stats_key].to_dict()],
+                    distances=distances,
+                    angular=angular,
                 )
+                stats_result = stats_results[0]
+                for dist_key in distances:
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_sum", dist_key, angular=angular, weighted=False)],
+                        stats_result.sum[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_sum", dist_key, angular=angular, weighted=True)],
+                        stats_result.sum_wt[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_mean", dist_key, angular=angular, weighted=False)],
+                        stats_result.mean[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_mean", dist_key, angular=angular, weighted=True)],
+                        stats_result.mean_wt[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_count", dist_key, angular=angular, weighted=False)],
+                        stats_result.count[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_count", dist_key, angular=angular, weighted=True)],
+                        stats_result.count_wt[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_var", dist_key, angular=angular, weighted=False)],
+                        stats_result.variance[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_var", dist_key, angular=angular, weighted=True)],
+                        stats_result.variance_wt[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_max", dist_key, angular=angular)],
+                        stats_result.max[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
+                    assert np.allclose(
+                        nodes_gdf[config.prep_gdf_key(f"{stats_key}_min", dist_key, angular=angular)],
+                        stats_result.min[dist_key],
+                        atol=config.ATOL,
+                        rtol=config.RTOL,
+                        equal_nan=True,
+                    )
     # check that problematic column labels are raised
     with pytest.raises(ValueError):
         layers.compute_stats(
             data_gdf,
-            "typo",
+            ["typo"],
             nodes_gdf,
             network_structure,
             distances=distances,
diff --git a/tests/rustalgos/test_data.py b/tests/rustalgos/test_data.py
index 5dc36918..b50d659d 100644
--- a/tests/rustalgos/test_data.py
+++ b/tests/rustalgos/test_data.py
@@ -308,12 +308,12 @@ def test_stats(primal_graph):
     # generate node and edge maps
     # generate node and edge maps
     _nodes_gdf, _edges_gdf, network_structure = io.network_structure_from_nx(primal_graph, 3395)
-    data_gdf = mock.mock_numerical_data(primal_graph, num_arrs=1, random_seed=13)
+    data_gdf = mock.mock_numerical_data(primal_graph, num_arrs=2, random_seed=13)
     # use a large enough distance such that simple non-weighted checks can be run for max, mean, variance
     max_assign_dist = 3200
     # don't deduplicate with data_id column otherwise below tallys won't work
     data_map, data_gdf = layers.assign_gdf_to_network(data_gdf, network_structure, max_assign_dist)
-    numerical_map = data_gdf["mock_numerical_1"].to_dict()
+    numerical_maps = [data_gdf["mock_numerical_1"].to_dict(), data_gdf["mock_numerical_2"].to_dict()]
     # for debugging
     # from cityseer.tools import plot
     # plot.plot_network_structure(network_structure, data_gdf)
@@ -328,165 +328,170 @@ def test_stats(primal_graph):
     # isolated loop = 52, 53, 54, 55 -> assigned data points = 1, 16, 24, 31, 36, 37
     isolated_nodes_idx = [52, 53, 54, 55]
     isolated_data_idx = [1, 16, 24, 31, 36, 37]
-    # numeric precision - keep fairly relaxed
-    mock_num_arr = data_gdf["mock_numerical_1"].values
     # compute - first do with no deduplication so that direct comparisons can be made to numpy methods
     # have to use a single large distance, otherwise distance cutoffs will result in limited agg
     distances = [10000]
-    stats_result = data_map.stats(
+    stats_results = data_map.stats(
         network_structure,
-        numerical_map=numerical_map,
+        numerical_maps=numerical_maps,
         distances=distances,
     )
-    for dist_key in distances:
-        # i.e. this scenarios considers all datapoints as unique (no two datapoints point to the same source)
-        # max
-        assert np.isnan(stats_result.max[dist_key][49])
-        assert np.allclose(
-            stats_result.max[dist_key][[50, 51]],
-            mock_num_arr[[33, 44]].max(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.max[dist_key][isolated_nodes_idx],
-            mock_num_arr[isolated_data_idx].max(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.max[dist_key][connected_nodes_idx],
-            mock_num_arr[connected_data_idx].max(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        # min
-        assert np.isnan(stats_result.max[dist_key][49])
-        assert np.allclose(
-            stats_result.min[dist_key][[50, 51]],
-            mock_num_arr[[33, 44]].min(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.min[dist_key][isolated_nodes_idx],
-            mock_num_arr[isolated_data_idx].min(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.min[dist_key][connected_nodes_idx],
-            mock_num_arr[connected_data_idx].min(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        # sum
-        assert np.isnan(stats_result.max[dist_key][49])
-        assert np.allclose(
-            stats_result.sum[dist_key][[50, 51]],
-            mock_num_arr[[33, 44]].sum(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.sum[dist_key][isolated_nodes_idx],
-            mock_num_arr[isolated_data_idx].sum(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.sum[dist_key][connected_nodes_idx],
-            mock_num_arr[connected_data_idx].sum(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        # mean
-        assert np.isnan(stats_result.max[dist_key][49])
-        assert np.allclose(
-            stats_result.mean[dist_key][[50, 51]],
-            mock_num_arr[[33, 44]].mean(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.mean[dist_key][isolated_nodes_idx],
-            mock_num_arr[isolated_data_idx].mean(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.mean[dist_key][connected_nodes_idx],
-            mock_num_arr[connected_data_idx].mean(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        # variance
-        assert np.isnan(stats_result.max[dist_key][49])
-        assert np.allclose(
-            stats_result.variance[dist_key][[50, 51]],
-            mock_num_arr[[33, 44]].var(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.variance[dist_key][isolated_nodes_idx],
-            mock_num_arr[isolated_data_idx].var(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
-        assert np.allclose(
-            stats_result.variance[dist_key][connected_nodes_idx],
-            mock_num_arr[connected_data_idx].var(),
-            atol=config.ATOL,
-            rtol=config.RTOL,
-        )
+    stats_result = stats_results[0]
+    for stats_result, mock_num_arr in zip(
+        stats_results, [data_gdf["mock_numerical_1"].values, data_gdf["mock_numerical_2"].values], strict=False
+    ):
+        for dist_key in distances:
+            # i.e. this scenarios considers all datapoints as unique (no two datapoints point to the same source)
+            # max
+            assert np.isnan(stats_result.max[dist_key][49])
+            assert np.allclose(
+                stats_result.max[dist_key][[50, 51]],
+                mock_num_arr[[33, 44]].max(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.max[dist_key][isolated_nodes_idx],
+                mock_num_arr[isolated_data_idx].max(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.max[dist_key][connected_nodes_idx],
+                mock_num_arr[connected_data_idx].max(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            # min
+            assert np.isnan(stats_result.max[dist_key][49])
+            assert np.allclose(
+                stats_result.min[dist_key][[50, 51]],
+                mock_num_arr[[33, 44]].min(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.min[dist_key][isolated_nodes_idx],
+                mock_num_arr[isolated_data_idx].min(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.min[dist_key][connected_nodes_idx],
+                mock_num_arr[connected_data_idx].min(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            # sum
+            assert np.isnan(stats_result.max[dist_key][49])
+            assert np.allclose(
+                stats_result.sum[dist_key][[50, 51]],
+                mock_num_arr[[33, 44]].sum(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.sum[dist_key][isolated_nodes_idx],
+                mock_num_arr[isolated_data_idx].sum(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.sum[dist_key][connected_nodes_idx],
+                mock_num_arr[connected_data_idx].sum(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            # mean
+            assert np.isnan(stats_result.max[dist_key][49])
+            assert np.allclose(
+                stats_result.mean[dist_key][[50, 51]],
+                mock_num_arr[[33, 44]].mean(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.mean[dist_key][isolated_nodes_idx],
+                mock_num_arr[isolated_data_idx].mean(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.mean[dist_key][connected_nodes_idx],
+                mock_num_arr[connected_data_idx].mean(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            # variance
+            assert np.isnan(stats_result.max[dist_key][49])
+            assert np.allclose(
+                stats_result.variance[dist_key][[50, 51]],
+                mock_num_arr[[33, 44]].var(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.variance[dist_key][isolated_nodes_idx],
+                mock_num_arr[isolated_data_idx].var(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
+            assert np.allclose(
+                stats_result.variance[dist_key][connected_nodes_idx],
+                mock_num_arr[connected_data_idx].var(),
+                atol=config.ATOL,
+                rtol=config.RTOL,
+            )
     # do deduplication - the stats should now be lower on average
     # the last five datapoints are pointing to the same source
     data_map_dedupe, data_gdf_dedupe = layers.assign_gdf_to_network(
         data_gdf, network_structure, max_assign_dist, data_id_col="data_id"
     )
-    stats_result_dedupe = data_map_dedupe.stats(
+    stats_results_dedupe = data_map_dedupe.stats(
         network_structure,
-        numerical_map=numerical_map,
+        numerical_maps=numerical_maps,
         distances=distances,
     )
-    for dist_key in distances:
-        # min and max are be the same
-        assert np.allclose(
-            stats_result.min[dist_key],
-            stats_result_dedupe.min[dist_key],
-            rtol=config.RTOL,
-            atol=config.ATOL,
-            equal_nan=True,
-        )
-        assert np.allclose(
-            stats_result.max[dist_key],
-            stats_result_dedupe.max[dist_key],
-            rtol=config.RTOL,
-            atol=config.ATOL,
-            equal_nan=True,
-        )
-        # sum should be lower when deduplicated
-        assert np.all(
-            stats_result.sum[dist_key][connected_nodes_idx] >= stats_result_dedupe.sum[dist_key][connected_nodes_idx]
-        )
-        assert np.all(
-            stats_result.sum_wt[dist_key][connected_nodes_idx]
-            >= stats_result_dedupe.sum_wt[dist_key][connected_nodes_idx]
-        )
-        # mean and variance should also be diminished
-        assert np.all(
-            stats_result.mean[dist_key][connected_nodes_idx] >= stats_result_dedupe.mean[dist_key][connected_nodes_idx]
-        )
-        assert np.all(
-            stats_result.mean_wt[dist_key][connected_nodes_idx]
-            >= stats_result_dedupe.mean_wt[dist_key][connected_nodes_idx]
-        )
-        assert np.all(
-            stats_result.variance[dist_key][connected_nodes_idx]
-            >= stats_result_dedupe.variance[dist_key][connected_nodes_idx]
-        )
-        assert np.all(
-            stats_result.variance_wt[dist_key][connected_nodes_idx]
-            >= stats_result_dedupe.variance_wt[dist_key][connected_nodes_idx]
-        )
+    for stats_result, stats_result_dedupe in zip(stats_results, stats_results_dedupe, strict=False):
+        for dist_key in distances:
+            # min and max are be the same
+            assert np.allclose(
+                stats_result.min[dist_key],
+                stats_result_dedupe.min[dist_key],
+                rtol=config.RTOL,
+                atol=config.ATOL,
+                equal_nan=True,
+            )
+            assert np.allclose(
+                stats_result.max[dist_key],
+                stats_result_dedupe.max[dist_key],
+                rtol=config.RTOL,
+                atol=config.ATOL,
+                equal_nan=True,
+            )
+            # sum should be lower when deduplicated
+            assert np.all(
+                stats_result.sum[dist_key][connected_nodes_idx]
+                >= stats_result_dedupe.sum[dist_key][connected_nodes_idx]
+            )
+            assert np.all(
+                stats_result.sum_wt[dist_key][connected_nodes_idx]
+                >= stats_result_dedupe.sum_wt[dist_key][connected_nodes_idx]
+            )
+            # mean and variance should also be diminished
+            assert np.all(
+                stats_result.mean[dist_key][connected_nodes_idx]
+                >= stats_result_dedupe.mean[dist_key][connected_nodes_idx]
+            )
+            assert np.all(
+                stats_result.mean_wt[dist_key][connected_nodes_idx]
+                >= stats_result_dedupe.mean_wt[dist_key][connected_nodes_idx]
+            )
+            assert np.all(
+                stats_result.variance[dist_key][connected_nodes_idx]
+                >= stats_result_dedupe.variance[dist_key][connected_nodes_idx]
+            )
+            assert np.all(
+                stats_result.variance_wt[dist_key][connected_nodes_idx]
+                >= stats_result_dedupe.variance_wt[dist_key][connected_nodes_idx]
+            )
diff --git a/tests/test_performance.py b/tests/test_performance.py
index 8c3aec14..f55957c7 100644
--- a/tests/test_performance.py
+++ b/tests/test_performance.py
@@ -57,11 +57,11 @@ def test_local_centrality_time(primal_graph):
     segment_cent_wrapper: 5.882331402972341 for 10000 iterations
 
     Split functions
-    dijkstra_tree_shortest_wrapper: 0.08737008203752339 for 10000 iterations
-    dijkstra_tree_simplest_wrapper: 0.09354790102224797 for 10000 iterations
-    dijkstra_tree_segment_wrapper: 0.2544154430506751 for 10000 iterations
-    node_cent_wrapper: 2.757581104990095 for 10000 iterations
-    segment_cent_wrapper: 5.390009218011983 for 10000 iterations
+    dijkstra_tree_shortest_wrapper: 0.04688391700619832 for 10000 iterations
+    dijkstra_tree_simplest_wrapper: 0.04833241600135807 for 10000 iterations
+    dijkstra_tree_segment_wrapper: 0.12999495898839086 for 10000 iterations
+    node_cent_wrapper: 3.6011423749878304 for 10000 iterations
+    segment_cent_wrapper: 4.1460652090027 for 10000 iterations
     """
 
     if "GITHUB_ACTIONS" in os.environ:
@@ -158,3 +158,5 @@ def segment_cent_wrapper():
     G_primal = mock_graph()
     G_primal = graphs.nx_simple_geoms(G_primal)
     test_local_centrality_time(G_primal)
+
+# %%