From 46e3b189ce28f7b827159571bde762c30a29bdc0 Mon Sep 17 00:00:00 2001 From: RunasSudo Date: Mon, 17 Apr 2023 17:50:43 +1000 Subject: [PATCH] Initial commit --- .gitignore | 2 + Cargo.lock | 938 ++++++++++++++++++++++++++++++++++++++++++++++++++ Cargo.toml | 17 + src/intcox.rs | 596 ++++++++++++++++++++++++++++++++ src/main.rs | 40 +++ 5 files changed, 1593 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 src/intcox.rs create mode 100644 src/main.rs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..54466f5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/target + diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..9ca76c9 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,938 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "anstream" +version = "0.2.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "342258dd14006105c2b75ab1bd7543a03bdf0cfc94383303ac212a04939dff6f" +dependencies = [ + "anstyle", + "anstyle-parse", + "anstyle-wincon", + "concolor-override", + "concolor-query", + "is-terminal", + "utf8parse", +] + +[[package]] +name = "anstyle" +version = "0.3.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23ea9e81bd02e310c216d080f6223c179012256e5151c41db88d12c88a1684d2" + +[[package]] +name = "anstyle-parse" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a7d1bb534e9efed14f3e5f44e7dd1a4f709384023a4165199a4241e18dff0116" +dependencies = [ + "utf8parse", +] + +[[package]] +name = "anstyle-wincon" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3127af6145b149f3287bb9a0d10ad9c5692dba8c53ad48285e5bec4063834fa" +dependencies = [ + "anstyle", + "windows-sys 0.45.0", +] + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "bytemuck" +version = "1.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "17febce684fd15d89027105661fec94afb475cb995fbc59d2865198446ba2eea" + +[[package]] +name = "cc" +version = "1.0.79" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "50d30906286121d95be3d479533b458f87493b30a4b5f79a607db8f5d11aa91f" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "clap" +version = "4.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "046ae530c528f252094e4a77886ee1374437744b2bff1497aa898bbddbbb29b3" +dependencies = [ + "clap_builder", + "clap_derive", + "once_cell", +] + +[[package]] +name = "clap_builder" +version = "4.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "223163f58c9a40c3b0a43e1c4b50a9ce09f007ea2cb1ec258a687945b4b7929f" +dependencies = [ + "anstream", + "anstyle", + "bitflags", + "clap_lex", + "strsim", +] + +[[package]] +name = "clap_derive" +version = "4.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f9644cd56d6b87dbe899ef8b053e331c0637664e9e21a33dfcdc36093f5c5c4" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn 2.0.14", +] + +[[package]] +name = "clap_lex" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8a2dd5a6fe8c6e3502f568a6353e5273bbb15193ad9a89e457b9970798efbea1" + +[[package]] +name = "concolor-override" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a855d4a1978dc52fb0536a04d384c2c0c1aa273597f08b77c8c4d3b2eec6037f" + +[[package]] +name = "concolor-query" +version = "0.3.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "88d11d52c3d7ca2e6d0040212be9e4dbbcd78b6447f535b6b561f449427944cf" +dependencies = [ + "windows-sys 0.45.0", +] + +[[package]] +name = "console" +version = "0.15.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3d79fbe8970a77e3e34151cc13d3b3e248aa0faaecb9f6091fa07ebefe5ad60" +dependencies = [ + "encode_unicode 0.3.6", + "lazy_static", + "libc", + "unicode-width", + "windows-sys 0.42.0", +] + +[[package]] +name = "crossbeam-channel" +version = "0.5.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a33c2bf77f2df06183c3aa30d1e96c0695a313d4f9c453cc3762a6db39f99200" +dependencies = [ + "cfg-if", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-deque" +version = "0.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ce6fd6f855243022dcecf8702fef0c297d4338e226845fe067f6341ad9fa0cef" +dependencies = [ + "cfg-if", + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "46bd5f3f85273295a9d14aedfb86f6aadbff6d8f5295c4a9edb08e819dcf5695" +dependencies = [ + "autocfg", + "cfg-if", + "crossbeam-utils", + "memoffset", + "scopeguard", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c063cd8cc95f5c377ed0d4b49a4b21f632396ff690e8470c29b3359b346984b" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "csv" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b015497079b9a9d69c02ad25de6c0a6edef051ea6360a327d0bd05802ef64ad" +dependencies = [ + "csv-core", + "itoa", + "ryu", + "serde", +] + +[[package]] +name = "csv-core" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b2466559f260f48ad25fe6317b3c8dac77b5bdb5763ac7d9d6103530663bc90" +dependencies = [ + "memchr", +] + +[[package]] +name = "dirs-next" +version = "2.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b98cf8ebf19c3d1b223e151f99a4f9f0690dca41414773390fc824184ac833e1" +dependencies = [ + "cfg-if", + "dirs-sys-next", +] + +[[package]] +name = "dirs-sys-next" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ebda144c4fe02d1f7ea1a7d9641b6fc6b580adcfa024ae48797ecdeb6825b4d" +dependencies = [ + "libc", + "redox_users", + "winapi", +] + +[[package]] +name = "either" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7fcaabb2fef8c910e7f4c7ce9f67a1283a1715879a7c230ca9d6d1ae31f16d91" + +[[package]] +name = "encode_unicode" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a357d28ed41a50f9c765dbfe56cbc04a64e53e5fc58ba79fbc34c10ef3df831f" + +[[package]] +name = "encode_unicode" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34aa73646ffb006b8f5147f3dc182bd4bcb190227ce861fc4a4844bf8e3cb2c0" + +[[package]] +name = "errno" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4bcfec3a70f97c962c307b2d2c56e358cf1d00b558d74262b5f929ee8cc7e73a" +dependencies = [ + "errno-dragonfly", + "libc", + "windows-sys 0.48.0", +] + +[[package]] +name = "errno-dragonfly" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa68f1b12764fab894d2755d2518754e71b4fd80ecfb822714a1206c2aab39bf" +dependencies = [ + "cc", + "libc", +] + +[[package]] +name = "getrandom" +version = "0.2.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c85e1d9ab2eadba7e5040d4e09cbd6d072b76a557ad64e797c2cb9d4da21d7e4" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "hermit-abi" +version = "0.2.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ee512640fe35acbfb4bb779db6f0d80704c2cacfa2e39b601ef3e3f47d1ae4c7" +dependencies = [ + "libc", +] + +[[package]] +name = "hermit-abi" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fed44880c466736ef9a5c5b5facefb5ed0785676d0c02d612db14e54f0d84286" + +[[package]] +name = "hpstat" +version = "0.1.0" +dependencies = [ + "clap", + "indicatif", + "nalgebra", + "prettytable-rs", + "rayon", + "serde", + "serde_json", +] + +[[package]] +name = "indicatif" +version = "0.17.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cef509aa9bc73864d6756f0d34d35504af3cf0844373afe9b8669a5b8005a729" +dependencies = [ + "console", + "number_prefix", + "portable-atomic", + "rayon", + "unicode-width", +] + +[[package]] +name = "io-lifetimes" +version = "1.0.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9c66c74d2ae7e79a5a8f7ac924adbe38ee42a859c6539ad869eb51f0b52dc220" +dependencies = [ + "hermit-abi 0.3.1", + "libc", + "windows-sys 0.48.0", +] + +[[package]] +name = "is-terminal" +version = "0.4.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adcf93614601c8129ddf72e2d5633df827ba6551541c6d8c59520a371475be1f" +dependencies = [ + "hermit-abi 0.3.1", + "io-lifetimes", + "rustix", + "windows-sys 0.48.0", +] + +[[package]] +name = "itoa" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "453ad9f582a441959e5f0d088b02ce04cfe8d51a8eaf077f12ac6d3e94164ca6" + +[[package]] +name = "lazy_static" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" + +[[package]] +name = "libc" +version = "0.2.141" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3304a64d199bb964be99741b7a14d26972741915b3649639149b2479bb46f4b5" + +[[package]] +name = "linux-raw-sys" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d59d8c75012853d2e872fb56bc8a2e53718e2cafe1a4c823143141c6d90c322f" + +[[package]] +name = "matrixmultiply" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "add85d4dd35074e6fedc608f8c8f513a3548619a9024b751949ef0e8e45a4d84" +dependencies = [ + "rawpointer", +] + +[[package]] +name = "memchr" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2dffe52ecf27772e601905b7522cb4ef790d2cc203488bbd0e2fe85fcb74566d" + +[[package]] +name = "memoffset" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d61c719bcfbcf5d62b3a09efa6088de8c54bc0bfcd3ea7ae39fcc186108b8de1" +dependencies = [ + "autocfg", +] + +[[package]] +name = "nalgebra" +version = "0.32.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d68d47bba83f9e2006d117a9a33af1524e655516b8919caac694427a6fb1e511" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d232c68884c0c99810a5a4d333ef7e47689cfd0edc85efc9e54e1e6bf5212766" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "num-complex" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0d21255c828d6f128a1e41534206671e8c3ea0c62f32291e808dc82cff17d" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "578ede34cf02f8924ab9447f50c28075b4d3e5b269972345e7e0372b38c6cdcd" +dependencies = [ + "autocfg", +] + +[[package]] +name = "num_cpus" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fac9e2da13b5eb447a6ce3d392f23a29d8694bff781bf03a16cd9ac8697593b" +dependencies = [ + "hermit-abi 0.2.6", + "libc", +] + +[[package]] +name = "number_prefix" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "830b246a0e5f20af87141b25c173cd1b609bd7779a4617d6ec582abaf90870f3" + +[[package]] +name = "once_cell" +version = "1.17.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b7e5500299e16ebb147ae15a00a942af264cf3688f47923b8fc2cd5858f23ad3" + +[[package]] +name = "paste" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f746c4065a8fa3fe23974dd82f15431cc8d40779821001404d10d2e79ca7d79" + +[[package]] +name = "portable-atomic" +version = "0.3.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "26f6a7b87c2e435a3241addceeeff740ff8b7e76b74c13bf9acb17fa454ea00b" + +[[package]] +name = "prettytable-rs" +version = "0.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eea25e07510aa6ab6547308ebe3c036016d162b8da920dbb079e3ba8acf3d95a" +dependencies = [ + "csv", + "encode_unicode 1.0.0", + "is-terminal", + "lazy_static", + "term", + "unicode-width", +] + +[[package]] +name = "proc-macro2" +version = "1.0.56" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b63bdb0cd06f1f4dedf69b254734f9b45af66e4a031e42a7480257d9898b435" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.26" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4424af4bf778aae2051a77b60283332f386554255d722233d09fbfc7e30da2fc" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "rayon" +version = "1.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d2df5196e37bcc87abebc0053e20787d73847bb33134a69841207dd0a47f03b" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b8f95bd6966f5c87776639160a66bd8ab9895d9d4ab01ddba9fc60661aebe8d" +dependencies = [ + "crossbeam-channel", + "crossbeam-deque", + "crossbeam-utils", + "num_cpus", +] + +[[package]] +name = "redox_syscall" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb5a58c1855b4b6819d59012155603f0b22ad30cad752600aadfcb695265519a" +dependencies = [ + "bitflags", +] + +[[package]] +name = "redox_users" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b033d837a7cf162d7993aded9304e30a83213c648b6e389db233191f891e5c2b" +dependencies = [ + "getrandom", + "redox_syscall", + "thiserror", +] + +[[package]] +name = "rustix" +version = "0.37.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85597d61f83914ddeba6a47b3b8ffe7365107221c2e557ed94426489fefb5f77" +dependencies = [ + "bitflags", + "errno", + "io-lifetimes", + "libc", + "linux-raw-sys", + "windows-sys 0.48.0", +] + +[[package]] +name = "rustversion" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4f3208ce4d8448b3f3e7d168a73f5e0c43a61e32930de3bceeccedb388b6bf06" + +[[package]] +name = "ryu" +version = "1.0.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f91339c0467de62360649f8d3e185ca8de4224ff281f66000de5eb2a77a79041" + +[[package]] +name = "safe_arch" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "794821e4ccb0d9f979512f9c1973480123f9bd62a90d74ab0f9426fcf8f4a529" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "scopeguard" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" + +[[package]] +name = "serde" +version = "1.0.160" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb2f3770c8bce3bcda7e149193a069a0f4365bda1fa5cd88e03bca26afc1216c" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.160" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "291a097c63d8497e00160b166a967a4a79c64f3facdd01cbd7502231688d77df" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.14", +] + +[[package]] +name = "serde_json" +version = "1.0.96" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "057d394a50403bcac12672b2b18fb387ab6d289d957dab67dd201875391e52f1" +dependencies = [ + "itoa", + "ryu", + "serde", +] + +[[package]] +name = "simba" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "061507c94fc6ab4ba1c9a0305018408e312e17c041eb63bef8aa726fa33aceae" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + +[[package]] +name = "strsim" +version = "0.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73473c0e59e6d5812c5dfe2a064a6444949f089e20eec9a2e5506596494e4623" + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fcf316d5356ed6847742d036f8a39c3b8435cac10bd528a4bd461928a6ab34d5" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "term" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c59df8ac95d96ff9bede18eb7300b0fda5e5d8d90960e76f8e14ae765eedbf1f" +dependencies = [ + "dirs-next", + "rustversion", + "winapi", +] + +[[package]] +name = "thiserror" +version = "1.0.40" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "978c9a314bd8dc99be594bc3c175faaa9794be04a5a5e153caba6915336cebac" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.40" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f9456a42c5b0d803c8cd86e73dd7cc9edd429499f37a3550d286d5e86720569f" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.14", +] + +[[package]] +name = "typenum" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "497961ef93d974e23eb6f433eb5fe1b7930b659f06d12dec6fc44a8f554c0bba" + +[[package]] +name = "unicode-ident" +version = "1.0.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5464a87b239f13a63a501f2701565754bae92d243d4bb7eb12f6d57d2269bf4" + +[[package]] +name = "unicode-width" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c0edd1e5b14653f783770bce4a4dabb4a5108a5370a5f5d8cfe8710c361f6c8b" + +[[package]] +name = "utf8parse" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "711b9620af191e0cdc7468a8d14e709c3dcdb115b36f838e601583af800a370a" + +[[package]] +name = "wasi" +version = "0.11.0+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" + +[[package]] +name = "wide" +version = "0.7.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b689b6c49d6549434bf944e6b0f39238cf63693cb7a147e9d887507fffa3b223" +dependencies = [ + "bytemuck", + "safe_arch", +] + +[[package]] +name = "winapi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" + +[[package]] +name = "windows-sys" +version = "0.42.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a3e1820f08b8513f676f7ab6c1f99ff312fb97b553d30ff4dd86f9f15728aa7" +dependencies = [ + "windows_aarch64_gnullvm 0.42.2", + "windows_aarch64_msvc 0.42.2", + "windows_i686_gnu 0.42.2", + "windows_i686_msvc 0.42.2", + "windows_x86_64_gnu 0.42.2", + "windows_x86_64_gnullvm 0.42.2", + "windows_x86_64_msvc 0.42.2", +] + +[[package]] +name = "windows-sys" +version = "0.45.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "75283be5efb2831d37ea142365f009c02ec203cd29a3ebecbc093d52315b66d0" +dependencies = [ + "windows-targets 0.42.2", +] + +[[package]] +name = "windows-sys" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "677d2418bec65e3338edb076e806bc1ec15693c5d0104683f2efe857f61056a9" +dependencies = [ + "windows-targets 0.48.0", +] + +[[package]] +name = "windows-targets" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e5180c00cd44c9b1c88adb3693291f1cd93605ded80c250a75d472756b4d071" +dependencies = [ + "windows_aarch64_gnullvm 0.42.2", + "windows_aarch64_msvc 0.42.2", + "windows_i686_gnu 0.42.2", + "windows_i686_msvc 0.42.2", + "windows_x86_64_gnu 0.42.2", + "windows_x86_64_gnullvm 0.42.2", + "windows_x86_64_msvc 0.42.2", +] + +[[package]] +name = "windows-targets" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7b1eb6f0cd7c80c79759c929114ef071b87354ce476d9d94271031c0497adfd5" +dependencies = [ + "windows_aarch64_gnullvm 0.48.0", + "windows_aarch64_msvc 0.48.0", + "windows_i686_gnu 0.48.0", + "windows_i686_msvc 0.48.0", + "windows_x86_64_gnu 0.48.0", + "windows_x86_64_gnullvm 0.48.0", + "windows_x86_64_msvc 0.48.0", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "597a5118570b68bc08d8d59125332c54f1ba9d9adeedeef5b99b02ba2b0698f8" + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91ae572e1b79dba883e0d315474df7305d12f569b400fcf90581b06062f7e1bc" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e08e8864a60f06ef0d0ff4ba04124db8b0fb3be5776a5cd47641e942e58c4d43" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b2ef27e0d7bdfcfc7b868b317c1d32c641a6fe4629c171b8928c7b08d98d7cf3" + +[[package]] +name = "windows_i686_gnu" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c61d927d8da41da96a81f029489353e68739737d3beca43145c8afec9a31a84f" + +[[package]] +name = "windows_i686_gnu" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "622a1962a7db830d6fd0a69683c80a18fda201879f0f447f065a3b7467daa241" + +[[package]] +name = "windows_i686_msvc" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "44d840b6ec649f480a41c8d80f9c65108b92d89345dd94027bfe06ac444d1060" + +[[package]] +name = "windows_i686_msvc" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4542c6e364ce21bf45d69fdd2a8e455fa38d316158cfd43b3ac1c5b1b19f8e00" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8de912b8b8feb55c064867cf047dda097f92d51efad5b491dfb98f6bbb70cb36" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ca2b8a661f7628cbd23440e50b05d705db3686f894fc9580820623656af974b1" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "26d41b46a36d453748aedef1486d5c7a85db22e56aff34643984ea85514e94a3" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7896dbc1f41e08872e9d5e8f8baa8fdd2677f29468c4e156210174edc7f7b953" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.42.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9aec5da331524158c6d1a4ac0ab1541149c0b9505fde06423b02f5ef0106b9f0" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.48.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1a515f5799fe4961cb532f983ce2b23082366b898e52ffbce459c86f67c8378a" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..2da13fc --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,17 @@ +[package] +name = "hpstat" +version = "0.1.0" +edition = "2021" + +[dependencies] +clap = { version = "4.2.1", features = ["derive"] } +indicatif = {version = "0.17.3", features = ["rayon"]} +nalgebra = "0.32.2" +prettytable-rs = "0.10.0" +rayon = "1.7.0" +serde = { version = "1.0.160", features = ["derive"] } +serde_json = "1.0.96" + +[profile.perf] +inherits = "release" +debug = true diff --git a/src/intcox.rs b/src/intcox.rs new file mode 100644 index 0000000..716796a --- /dev/null +++ b/src/intcox.rs @@ -0,0 +1,596 @@ +// hpstat: High-performance statistics implementations +// Copyright © 2023 Lee Yingtong Li (RunasSudo) +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU Affero General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Affero General Public License for more details. +// +// You should have received a copy of the GNU Affero General Public License +// along with this program. If not, see . + +const Z_97_5: f64 = 1.959964; // This is the limit of resolution for an f64 + +use std::fs; +use std::io; + +use clap::{Args, ValueEnum}; +use indicatif::{ParallelProgressIterator, ProgressBar, ProgressStyle}; +use nalgebra::{DMatrix, DVector, Matrix1xX}; +use prettytable::{Table, format, row}; +use rayon::prelude::*; +use serde::{Serialize, Deserialize}; + +#[derive(Args)] +pub struct IntCoxArgs { + /// Path to CSV input file containing the observations + #[arg()] + input: String, + + /// Output format + #[arg(long, value_enum, default_value="text")] + output: OutputFormat, + + /// Maximum number of E-M iterations to attempt + #[arg(long, default_value="100")] + max_iterations: u32, + + /// Terminate E-M algorithm when the maximum absolute change in all parameters is less than this tolerance + #[arg(long, default_value="0.001")] + tolerance: f64, + + /// Estimate baseline hazard function using Turnbull innermost intervals + #[arg(long)] + reduced: bool, +} + +#[derive(ValueEnum, Clone)] +enum OutputFormat { + Text, + Json +} + +pub fn main(args: IntCoxArgs) { + let lines: Vec; + if args.input == "-" { + lines = io::stdin().lines().map(|l| l.unwrap()).collect(); + } else { + let contents = fs::read_to_string(args.input).unwrap(); + lines = contents.trim_end().split("\n").map(|s| s.to_string()).collect(); + } + + // Read data into matrices + + let mut data_times: DMatrix = DMatrix::zeros( + 2, // Left time, right time + lines.len() - 1 // Minus 1 row for header row + ); + + // Called "Z" in the paper and "X" in the C++ code + let mut data_indep: DMatrix = DMatrix::zeros( + lines[0].split(",").count() - 2, + lines.len() - 1 // Minus 1 row for header row + ); + + // Read header row + let indep_names: Vec<&str> = lines[0].split(",").skip(2).collect(); + + // Read data + // FIXME: Parse CSV more robustly + for (i, row) in lines.iter().skip(1).enumerate() { + for (j, item) in row.split(",").enumerate() { + let value = match item { + "inf" => f64::INFINITY, + _ => item.parse().expect("Malformed float") + }; + + if j < 2 { + data_times[(j, i)] = value; + } else { + data_indep[(j - 2, i)] = value; + } + } + } + + // Fit regression + let progress_bar = match args.output { + OutputFormat::Text => ProgressBar::new(0), + OutputFormat::Json => ProgressBar::hidden(), + }; + let result = fit_interval_censored_cox(data_times, data_indep, args.max_iterations, args.tolerance, args.reduced, progress_bar); + + // Display output + match args.output { + OutputFormat::Text => { + println!(); + println!(); + println!("LL-Model = {:.5}", result.ll_model); + println!("LL-Null = {:.5}", result.ll_null); + + let mut summary = Table::new(); + let format = format::FormatBuilder::new() + .separators(&[format::LinePosition::Top, format::LinePosition::Title, format::LinePosition::Bottom], format::LineSeparator::new('-', '+', '+', '+')) + .padding(2, 2) + .build(); + summary.set_format(format); + + summary.set_titles(row!["Parameter", c->"β", c->"Std Err.", c->"exp(β)", H2c->"(95% CI)"]); + for (i, indep_name) in indep_names.iter().enumerate() { + summary.add_row(row![ + indep_name, + r->format!("{:.5}", result.params[i]), + r->format!("{:.5}", result.params_se[i]), + r->format!("{:.5}", result.params[i].exp()), + r->format!("({:.5},", (result.params[i] - Z_97_5 * result.params_se[i]).exp()), + format!("{:.5})", (result.params[i] + Z_97_5 * result.params_se[i]).exp()), + ]); + } + summary.printstd(); + } + OutputFormat::Json => { + println!("{}", serde_json::to_string(&result).unwrap()); + } + } +} + +struct IntervalCensoredCoxData { + data_times: DMatrix, + data_indep: DMatrix, + + // Cached intermediate values + time_points: Vec, + r_star_indicator: DMatrix, + z_z_transpose: Vec>, +} + +impl IntervalCensoredCoxData { + fn num_obs(&self) -> usize { + return self.data_indep.ncols(); + } + + fn num_covs(&self) -> usize { + return self.data_indep.nrows(); + } + + fn num_times(&self) -> usize { + return self.time_points.len(); + } +} + +fn fit_interval_censored_cox(data_times: DMatrix, mut data_indep: DMatrix, max_iterations: u32, tolerance: f64, reduced: bool, progress_bar: ProgressBar) -> IntervalCensoredCoxResult { + // ---------------------- + // Prepare for regression + + // Standardise values + let indep_means = data_indep.column_mean(); + let indep_stdev = data_indep.column_variance().apply_into(|x| { *x = (*x * data_indep.ncols() as f64 / (data_indep.ncols() - 1) as f64).sqrt(); }); + for j in 0..data_indep.nrows() { + data_indep.row_mut(j).apply(|x| *x = (*x - indep_means[j]) / indep_stdev[j]); + } + + // Get time points (t_0 = 0, t_1, ..., t_m) + let mut time_points: Vec; + if reduced { + // Turnbull intervals + let mut all_time_points: Vec<(f64, bool)> = Vec::new(); // Vec of (time, is_left) + all_time_points.extend(data_times.row(0).iter().map(|t| (*t, true))); + all_time_points.extend(data_times.row(1).iter().map(|t| (*t, false))); + all_time_points.sort_by(|(t1, _), (t2, _)| t1.partial_cmp(t2).unwrap()); + + time_points = Vec::new(); + for i in 1..all_time_points.len() { + if all_time_points[i - 1].1 == true && all_time_points[i].1 == false { + time_points.push(all_time_points[i - 1].0); + time_points.push(all_time_points[i].0); + } + } + time_points.push(0.0); // Ensure 0 is in the list + time_points.retain(|t| t.is_finite()); // Remove infinity + time_points.sort_by(|a, b| a.partial_cmp(b).unwrap()); // Cannot use .sort() as f64 does not implement Ord + time_points.dedup(); + } else { + // All observed intervals + time_points = data_times.iter().copied().collect(); + time_points.push(0.0); // Ensure 0 is in the list + time_points.retain(|t| t.is_finite()); // Remove infinity + time_points.sort_by(|a, b| a.partial_cmp(b).unwrap()); // Cannot use .sort() as f64 does not implement Ord + time_points.dedup(); + } + + // Initialise β, λ + let mut beta = DVector::zeros(data_indep.nrows()); + let mut lambda = DVector::repeat(time_points.len(), 1.0 / (time_points.len() - 1) as f64); + + // Compute I(t_k <= R*_i) + // Where R*_i is R_i if R_i ≠ ∞, otherwise it is L_i + let mut r_star_indicator = DMatrix::zeros(data_indep.ncols(), time_points.len()); + for (i, observation) in data_times.column_iter().enumerate() { + let time_right_star = if observation[1].is_finite() { observation[1] } else { observation[0] }; + for (k, time) in time_points.iter().enumerate() { + if *time <= time_right_star { + // t_k <= R*_i + r_star_indicator[(i, k)] = 1.0; + } else { + r_star_indicator[(i, k)] = 0.0; + } + } + } + + // Pre-compute Z * Z^T + // Indexed by observation -> Matrix (num-covariates, num-covariates) + let mut z_z_transpose: Vec> = Vec::new(); + for i in 0..data_indep.ncols() { + let covariates = data_indep.column(i); + z_z_transpose.push(covariates * covariates.transpose()); + } + + let data = IntervalCensoredCoxData { + data_times: data_times, + data_indep: data_indep, + time_points: time_points, + r_star_indicator: r_star_indicator, + z_z_transpose: z_z_transpose, + }; + + // ------------------- + // Apply E-M algorithm + + progress_bar.set_length(u64::MAX); + progress_bar.reset(); + progress_bar.set_style(ProgressStyle::with_template("[{elapsed_precise}] {bar:40} {msg}").unwrap()); + progress_bar.println("Running E-M algorithm to fit interval-censored Cox model"); + + let mut iteration: u32 = 0; + loop { + // Pre-compute exp(β^T * Z_ik) + let exp_beta_z: Matrix1xX = (beta.transpose() * &data.data_indep).apply_into(|x| { *x = x.exp(); }); + + // Do E-step + let posterior_weight = do_e_step(&data, &exp_beta_z, &lambda); + + // Do M-step + let (new_beta, new_lambda) = do_m_step(&data, &exp_beta_z, &beta, posterior_weight); + + // Check for convergence + let (coef_change, converged) = em_check_convergence(&beta, &lambda, &new_beta, &new_lambda, tolerance); + beta = new_beta; + lambda = new_lambda; + + // Update progress bar + // Estimate progress according to either the order of magnitude of the coef_change relative to tolerance, or iteration/max_iterations + let progress1 = ((-coef_change.log10()).max(0.0) / -tolerance.log10() * u64::MAX as f64) as u64; + let progress2 = (iteration as f64 / max_iterations as f64 * u64::MAX as f64) as u64; + progress_bar.set_position(progress_bar.position().max(progress1).max(progress2)); + progress_bar.set_message(format!("Iter {} (delta = {:.6})", iteration + 1, coef_change)); + + if converged { + progress_bar.finish(); + break; + } + + iteration += 1; + if iteration >= max_iterations { + panic!("Exceeded --max-iterations"); + } + } + + // Compute log-likelihood + let ll_model = log_likelihood_obs(&data, &beta, &lambda).sum(); + + // Unstandardise betas + let mut beta_unstandardised: DVector = DVector::zeros(data.num_covs()); + for (j, beta_value) in beta.iter().enumerate() { + beta_unstandardised[j] = beta_value / indep_stdev[j]; + } + + // ------------------------- + // Compute covariance matrix + + // Compute profile log-likelihoods + let h = 5.0 / (data.num_obs() as f64).sqrt(); // "a constant of order n^(-1/2)" + + progress_bar.set_length(data.num_covs() as u64 + 2); + progress_bar.reset(); + progress_bar.set_style(ProgressStyle::with_template("[{elapsed_precise}] {bar:40} Profile LL {pos}/{len}").unwrap()); + progress_bar.println("Profiling log-likelihood to compute covariance matrix"); + + // ll_null = log-likelihood for null model + // pll_toggle_zero = log-likelihoods for each observation at final beta + // pll_toggle_one = log-likelihoods for each observation at toggled beta + + let ll_null = profile_log_likelihood_obs(&data, DVector::zeros(data.num_covs()), lambda.clone(), max_iterations, tolerance).sum(); + + let pll_toggle_zero: DVector = profile_log_likelihood_obs(&data, beta.clone(), lambda.clone(), max_iterations, tolerance); + progress_bar.inc(1); + + let pll_toggle_one: Vec> = (0..data.num_covs()).into_par_iter().map(|j| { + let mut pll_beta = beta.clone(); + pll_beta[j] += h; + profile_log_likelihood_obs(&data, pll_beta, lambda.clone(), max_iterations, tolerance) + }) + .progress_with(progress_bar.clone()) + .collect(); + progress_bar.finish(); + + let mut pll_matrix: DMatrix = DMatrix::zeros(data.num_covs(), data.num_covs()); + for i in 0..data.num_obs() { + let toggle_none_i = pll_toggle_zero[i]; + let mut ps_i: DVector = DVector::zeros(data.num_covs()); + for p in 0..data.num_covs() { + ps_i[p] = (pll_toggle_one[p][i] - toggle_none_i) / h; + } + pll_matrix += ps_i.clone() * ps_i.transpose(); + } + + let vcov = pll_matrix.try_inverse().expect("Matrix not invertible"); + + // Unstandardise SEs + let beta_se = vcov.diagonal().apply_into(|x| { *x = x.sqrt(); }); + let mut beta_se_unstandardised: DVector = DVector::zeros(data.num_covs()); + for (j, se) in beta_se.iter().enumerate() { + beta_se_unstandardised[j] = se / indep_stdev[j]; + } + + return IntervalCensoredCoxResult { + params: beta_unstandardised.data.as_vec().clone(), + params_se: beta_se_unstandardised.data.as_vec().clone(), + ll_model: ll_model, + ll_null: ll_null, + }; +} + +fn do_e_step(data: &IntervalCensoredCoxData, exp_beta_z: &Matrix1xX, lambda: &DVector) -> DMatrix { + // Compute S_L and S_R (S_i1 and S_i2 in the paper) + let s_left = e_step_compute_s(data, &exp_beta_z, lambda, 0); + let s_right = e_step_compute_s(data, &exp_beta_z, lambda, 1); + + // In the paper, consideration is given to G(x) + // But in a proportional hazards model, G(x) = x + // So we omit the details + // As a consequence, the posterior ξ_i are always 1 + + // Compute posterior weights (W_ik, "posterior mean" in C++) + let mut posterior_weight: DMatrix = DMatrix::zeros(data.num_obs(), data.num_times()); + for (i, observation) in data.data_times.column_iter().enumerate() { + let time_left = observation[0]; + let time_right = observation[1]; + + for (k, time) in data.time_points.iter().enumerate() { + if *time <= time_left { + // t_k <= L_i + posterior_weight[(i, k)] = 0.0; + } else if *time <= time_right && time_right.is_finite() { + // L_i < t_k <= R_i, with R_i < ∞ + // Assumes r = 0 + posterior_weight[(i, k)] = lambda[k] * exp_beta_z[i] / (1.0 - (s_left[i] - s_right[i]).exp()); + } else { + // None of the above circumstances + // C++ says the weight is unused in this case + // Set this to a non-NaN value so we can still do elementwise vector multiplication for masking + posterior_weight[(i, k)] = 0.0; + } + } + } + + return posterior_weight; +} + +fn e_step_compute_s(data: &IntervalCensoredCoxData, exp_beta_z: &Matrix1xX, lambda: &DVector, time_index: usize) -> DVector { + let mut s: DVector = DVector::zeros(data.num_obs()); + for (i, observation) in data.data_times.column_iter().enumerate() { + let time_cutoff = observation[time_index]; + if time_cutoff.is_infinite() { + s[i] = f64::INFINITY; + } else { + for (k, time) in data.time_points.iter().enumerate() { + if *time <= time_cutoff { + // time is t_k <= L_i, or t_k <= R_i, as applicable + s[i] += lambda[k] * exp_beta_z[i]; // Row 0, because all covariates are time-independent + } else { + break; + } + } + } + } + return s; +} + +fn do_m_step(data: &IntervalCensoredCoxData, exp_beta_z: &Matrix1xX, beta: &DVector, posterior_weight: DMatrix) -> (DVector, DVector) { + // ComputeSummandTerm + // Covariates are time-independent in this model + // And ξ_i is always 1, as discussed above + // So we can skip this step and let xi_exp_beta_z = exp_beta_z + let xi_exp_beta_z = &exp_beta_z; + + // Split these steps into functions to make profiling easier + let (mut s0, s1, s2) = m_step_compute_s_values(data, xi_exp_beta_z); + let sigma = m_step_compute_sigma(data, &posterior_weight, &s0, &s1, &s2); + let new_beta = m_step_compute_new_beta(data, &posterior_weight, &s0, &s1, sigma, beta); + s0 = m_step_compute_s0(data, beta); + let new_lambda = m_step_compute_new_lambda(data, &posterior_weight, &s0); + + return (new_beta, new_lambda); +} + +fn m_step_compute_s_values(data: &IntervalCensoredCoxData, xi_exp_beta_z: &Matrix1xX) -> (DVector, Vec>, Vec>) { + // ComputeSValues + + // Compute s0 + let mut s0: DVector = DVector::zeros(data.num_times()); // Elements are f64 + for i in 0..data.num_obs() { + let s0_contrib = xi_exp_beta_z[i]; + s0 += data.r_star_indicator.row(i).transpose() * s0_contrib; + } + + // Precompute s1, s2 contributions for each observation + let mut s1_contrib: Vec> = vec![DVector::zeros(data.num_covs()); data.num_obs()]; // Elements are DVector of len num-covariates + let mut s2_contrib: Vec> = vec![DMatrix::zeros(data.num_covs(), data.num_covs()); data.num_obs()]; // Elements are (num-covariates, num-covariates) + for i in 0..data.num_obs() { + s1_contrib[i] = xi_exp_beta_z[i] * data.data_indep.column(i); + s2_contrib[i] = xi_exp_beta_z[i] * &data.z_z_transpose[i]; // Observations are time-independent + } + + let s1 = (0..data.num_times()).into_par_iter().map(|k| { + let mut s1_k = DVector::zeros(data.num_covs()); + for i in 0..data.num_obs() { + if data.r_star_indicator[(i, k)] == 1.0 { + s1_k += &s1_contrib[i]; + } + } + s1_k + }).collect(); + + let s2 = (0..data.num_times()).into_par_iter().map(|k| { + let mut s2_k = DMatrix::zeros(data.num_covs(), data.num_covs()); + for i in 0..data.num_obs() { + if data.r_star_indicator[(i, k)] == 1.0 { + s2_k += &s2_contrib[i]; + } + } + s2_k + }).collect(); + + return (s0, s1, s2); +} + +fn m_step_compute_sigma(data: &IntervalCensoredCoxData, posterior_weight: &DMatrix, s0: &DVector, s1: &Vec>, s2: &Vec>) -> DMatrix { + // ComputeSigma + let mut sigma: DMatrix = DMatrix::zeros(data.num_covs(), data.num_covs()); + for k in 0..data.num_times() { + let factor_k = (s1[k].clone() / s0[k]) * (s1[k].transpose() / s0[k]) - (s2[k].clone() / s0[k]); + let sum_posterior_weight = data.r_star_indicator.column(k).component_mul(&posterior_weight.column(k)).sum(); + sigma += sum_posterior_weight * factor_k.clone(); + } + return sigma; +} + +fn m_step_compute_new_beta(data: &IntervalCensoredCoxData, posterior_weight: &DMatrix, s0: &DVector, s1: &Vec>, sigma: DMatrix, beta: &DVector) -> DVector { + // ComputeNewBeta + assert!(sigma.clone().full_piv_lu().is_invertible(), "Sigma is not invertible"); + + let mut sum: DVector = DVector::zeros(data.num_covs()); + for k in 0..data.num_times() { + let quotient_k = s1[k].clone() / s0[k]; + for i in 0..data.num_obs() { + if data.r_star_indicator[(i, k)] == 1.0 { + sum += posterior_weight[(i, k)] * (data.data_indep.column(i) - "ient_k); + } + } + } + + let new_beta = beta.clone() - sigma.try_inverse().unwrap() * sum; + return new_beta; +} + +fn m_step_compute_s0(data: &IntervalCensoredCoxData, beta: &DVector) -> DVector { + // ComputeS0 + let mut s0: DVector = DVector::zeros(data.num_times()); + for i in 0..data.num_obs() { + // let s0_contrib = posterior_xi[i] * self.beta.dot(&data_indep.column(i)).exp(); + let s0_contrib = beta.dot(&data.data_indep.column(i)).exp(); + s0 += data.r_star_indicator.row(i).transpose() * s0_contrib; + } + return s0; +} + +fn m_step_compute_new_lambda(data: &IntervalCensoredCoxData, posterior_weight: &DMatrix, s0: &DVector) -> DVector { + // ComputeNewLambda + let mut new_lambda: DVector = DVector::zeros(data.num_times()); + for k in 0..data.num_times() { + let mut numerator_k = 0.0; + for i in 0..data.num_obs() { + if data.r_star_indicator[(i, k)] == 1.0 { + numerator_k += posterior_weight[(i, k)]; + } + } + new_lambda[k] = numerator_k / s0[k]; + } + return new_lambda; +} + +fn em_check_convergence(beta: &DVector, lambda: &DVector, new_beta: &DVector, new_lambda: &DVector, tolerance: f64) -> (f64, bool) { + let beta_diff = max_abs_difference(beta, new_beta); + + let old_cumulative_hazard = cumulative_hazard(lambda); + let new_cumulative_hazard = cumulative_hazard(new_lambda); + let lambda_diff = max_abs_difference(&old_cumulative_hazard, &new_cumulative_hazard); + + let max_diff = beta_diff.max(lambda_diff); + + return (max_diff, max_diff < tolerance); +} + +fn log_likelihood_obs(data: &IntervalCensoredCoxData, beta: &DVector, lambda: &DVector) -> DVector { + // Pre-compute exp(β^T * Z_ik) + let exp_beta_z: Matrix1xX = (beta.transpose() * &data.data_indep).apply_into(|x| { *x = x.exp(); }); + + // Compute S_L and S_R (S_i1 and S_i2 in the paper) + let s_left = e_step_compute_s(data, &exp_beta_z, lambda, 0); + let s_right = e_step_compute_s(data, &exp_beta_z, lambda, 1); + + // Compute the log-likelihood by summing log-likelihood for each observation + // Assumes G(x) = x + let mut result = DVector::zeros(data.num_obs()); + for i in 0..data.num_obs() { + result[i] = ((-s_left[i]).exp() - (-s_right[i]).exp()).ln(); + } + return result; +} + +fn profile_log_likelihood_obs(data: &IntervalCensoredCoxData, beta: DVector, mut lambda: DVector, max_iterations: u32, tolerance: f64) -> DVector { + for _iteration in 0..max_iterations { + // Pre-compute exp(β^T * Z_ik) + let exp_beta_z: Matrix1xX = (beta.transpose() * &data.data_indep).apply_into(|x| { *x = x.exp(); }); + + // Do E-step + let posterior_weight = do_e_step(data, &exp_beta_z, &lambda); + + // Do M-step (skip expensive unnecessary steps) + let s0 = m_step_compute_s0(data, &beta); + let new_lambda = m_step_compute_new_lambda(data, &posterior_weight, &s0); + + // Check for convergence + let old_cumulative_hazard = cumulative_hazard(&lambda); + let new_cumulative_hazard = cumulative_hazard(&new_lambda); + let lambda_diff = max_abs_difference(&old_cumulative_hazard, &new_cumulative_hazard); + lambda = new_lambda; + + // TODO: Incorporate into progress bar + //println!("Profile iteration {}, estimates changed by {}", iteration + 1, lambda_diff); + + if lambda_diff < tolerance { + return log_likelihood_obs(data, &beta, &lambda); + } + } + + panic!("Exceeded --max-iterations"); +} + +#[derive(Serialize, Deserialize)] +struct IntervalCensoredCoxResult { + params: Vec, + params_se: Vec, + ll_model: f64, + ll_null: f64, + // TODO: cumulative hazard, etc. +} + +fn cumulative_hazard(lambda: &DVector) -> DVector { + let mut result = DVector::zeros(lambda.nrows()); + for (i, value) in lambda.iter().enumerate() { + if i > 0 { + result[i] += result[i - 1]; + } + result[i] += value; + } + return result; +} + +fn max_abs_difference(vector_old: &DVector, vector_new: &DVector) -> f64 { + return (vector_new - vector_old).abs().max(); +} diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..b28bb33 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,40 @@ +// hpstat: High-performance statistics implementations +// Copyright © 2023 Lee Yingtong Li (RunasSudo) +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU Affero General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Affero General Public License for more details. +// +// You should have received a copy of the GNU Affero General Public License +// along with this program. If not, see . + +use clap::{Parser, Subcommand}; + +mod intcox; + +#[derive(Parser)] +#[command(about="High-performance statistics implementations")] +struct MainArgs { + #[clap(subcommand)] + command: Command, +} + +#[derive(Subcommand)] +enum Command { + #[command(name="intcox", about="Interval-censored Cox regression", long_about="Fit a Cox proportional hazards model on time-independent interval-censored observations, using an E-M algorithm by Zeng, Mao & Lin (Biometrika 2016;103(2):253-71)")] + IntCox(intcox::IntCoxArgs), +} + +fn main() { + let args = MainArgs::parse(); + + match args.command { + Command::IntCox(intcox_args) => intcox::main(intcox_args), + } +}