Initial commit

This commit is contained in:
RunasSudo 2023-04-17 17:50:43 +10:00
commit 46e3b189ce
Signed by: RunasSudo
GPG Key ID: 7234E476BF21C61A
5 changed files with 1593 additions and 0 deletions

2
.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
/target

938
Cargo.lock generated Normal file
View File

@ -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"

17
Cargo.toml Normal file
View File

@ -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

596
src/intcox.rs Normal file
View File

@ -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 <https://www.gnu.org/licenses/>.
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<String>;
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<f64> = 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<f64> = 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<f64>,
data_indep: DMatrix<f64>,
// Cached intermediate values
time_points: Vec<f64>,
r_star_indicator: DMatrix<f64>,
z_z_transpose: Vec<DMatrix<f64>>,
}
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<f64>, mut data_indep: DMatrix<f64>, 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<f64>;
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<DMatrix<f64>> = 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<f64> = (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<f64> = 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<f64> = profile_log_likelihood_obs(&data, beta.clone(), lambda.clone(), max_iterations, tolerance);
progress_bar.inc(1);
let pll_toggle_one: Vec<DVector<f64>> = (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<f64> = 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<f64> = 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<f64> = 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<f64>, lambda: &DVector<f64>) -> DMatrix<f64> {
// 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<f64> = 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<f64>, lambda: &DVector<f64>, time_index: usize) -> DVector<f64> {
let mut s: DVector<f64> = 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<f64>, beta: &DVector<f64>, posterior_weight: DMatrix<f64>) -> (DVector<f64>, DVector<f64>) {
// 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<f64>) -> (DVector<f64>, Vec<DVector<f64>>, Vec<DMatrix<f64>>) {
// ComputeSValues
// Compute s0
let mut s0: DVector<f64> = 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<DVector<f64>> = vec![DVector::zeros(data.num_covs()); data.num_obs()]; // Elements are DVector of len num-covariates
let mut s2_contrib: Vec<DMatrix<f64>> = 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<f64>, s0: &DVector<f64>, s1: &Vec<DVector<f64>>, s2: &Vec<DMatrix<f64>>) -> DMatrix<f64> {
// ComputeSigma
let mut sigma: DMatrix<f64> = 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<f64>, s0: &DVector<f64>, s1: &Vec<DVector<f64>>, sigma: DMatrix<f64>, beta: &DVector<f64>) -> DVector<f64> {
// ComputeNewBeta
assert!(sigma.clone().full_piv_lu().is_invertible(), "Sigma is not invertible");
let mut sum: DVector<f64> = 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) - &quotient_k);
}
}
}
let new_beta = beta.clone() - sigma.try_inverse().unwrap() * sum;
return new_beta;
}
fn m_step_compute_s0(data: &IntervalCensoredCoxData, beta: &DVector<f64>) -> DVector<f64> {
// ComputeS0
let mut s0: DVector<f64> = 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<f64>, s0: &DVector<f64>) -> DVector<f64> {
// ComputeNewLambda
let mut new_lambda: DVector<f64> = 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<f64>, lambda: &DVector<f64>, new_beta: &DVector<f64>, new_lambda: &DVector<f64>, 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<f64>, lambda: &DVector<f64>) -> DVector<f64> {
// Pre-compute exp(β^T * Z_ik)
let exp_beta_z: Matrix1xX<f64> = (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<f64>, mut lambda: DVector<f64>, max_iterations: u32, tolerance: f64) -> DVector<f64> {
for _iteration in 0..max_iterations {
// Pre-compute exp(β^T * Z_ik)
let exp_beta_z: Matrix1xX<f64> = (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<f64>,
params_se: Vec<f64>,
ll_model: f64,
ll_null: f64,
// TODO: cumulative hazard, etc.
}
fn cumulative_hazard(lambda: &DVector<f64>) -> DVector<f64> {
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<f64>, vector_new: &DVector<f64>) -> f64 {
return (vector_new - vector_old).abs().max();
}

40
src/main.rs Normal file
View File

@ -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 <https://www.gnu.org/licenses/>.
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),
}
}