mirror of
https://github.com/dnlbauer/WHAM.git
synced 2026-09-10 22:25:31 +00:00
Compare commits
19 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
b130596fe8 | ||
|
|
05451a6cb2 | ||
|
|
bf8bf39319 | ||
|
6c63eda769 | ||
|
|
1bbafeba70 | ||
|
|
64f468e636 | ||
|
|
2f93b77879 | ||
|
|
b148998b11 | ||
|
|
87f3cceb6a | ||
|
|
7b39d603bd | ||
|
|
465f46b280 | ||
|
|
d221fc7195 | ||
|
|
dc218df051 | ||
|
|
1e27cfd457 | ||
|
|
7116142c2a | ||
|
|
fbdb451f6b | ||
|
|
d4119b2fae | ||
|
|
b9f6b24b20 | ||
|
|
5390b089ac |
3
.gitignore
vendored
3
.gitignore
vendored
@@ -1,4 +1,5 @@
|
||||
target
|
||||
reference
|
||||
.idea
|
||||
*bench.*
|
||||
*bench.*
|
||||
*.sublime*
|
||||
288
Cargo.lock
generated
288
Cargo.lock
generated
@@ -11,6 +11,21 @@ dependencies = [
|
||||
"pkg-config",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "addr2line"
|
||||
version = "0.13.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "1b6a2d3371669ab3ca9797670853d61402b03d0b4b9ebf33d677dfa720203072"
|
||||
dependencies = [
|
||||
"gimli",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "adler"
|
||||
version = "0.2.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ee2a4ec343196209d6594e19543ae87a39f96d5534d7174822a3ad825dd6ed7e"
|
||||
|
||||
[[package]]
|
||||
name = "ansi_term"
|
||||
version = "0.11.0"
|
||||
@@ -20,6 +35,12 @@ dependencies = [
|
||||
"winapi",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "assert_approx_eq"
|
||||
version = "1.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "3c07dab4369547dbe5114677b33fbbf724971019f3818172d59a97a61c774ffd"
|
||||
|
||||
[[package]]
|
||||
name = "atty"
|
||||
version = "0.2.14"
|
||||
@@ -33,32 +54,24 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "autocfg"
|
||||
version = "1.0.0"
|
||||
version = "1.0.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "f8aac770f1885fd7e387acedd76065302551364496e46b3dd00860b2f8359b9d"
|
||||
checksum = "cdb031dd78e28731d87d56cc8ffef4a8f36ca26c38fe2de700543e627f8a464a"
|
||||
|
||||
[[package]]
|
||||
name = "backtrace"
|
||||
version = "0.3.45"
|
||||
version = "0.3.53"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ad235dabf00f36301792cfe82499880ba54c6486be094d1047b02bacb67c14e8"
|
||||
checksum = "707b586e0e2f247cbde68cdd2c3ce69ea7b7be43e1c5b426e37c9319c4b9838e"
|
||||
dependencies = [
|
||||
"backtrace-sys",
|
||||
"cfg-if",
|
||||
"addr2line",
|
||||
"cfg-if 1.0.0",
|
||||
"libc",
|
||||
"miniz_oxide",
|
||||
"object",
|
||||
"rustc-demangle",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "backtrace-sys"
|
||||
version = "0.1.33"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e17b52e737c40a7d75abca20b29a19a0eb7ba9fc72c5a72dd282a0a3c2c0dc35"
|
||||
dependencies = [
|
||||
"cc",
|
||||
"libc",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "bitflags"
|
||||
version = "1.2.1"
|
||||
@@ -71,12 +84,6 @@ version = "1.0.12"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "aa9e1d9f7d49e289f36f19effbf3d5a5e30163ecf9c7a3c9be94d5374dec5b9a"
|
||||
|
||||
[[package]]
|
||||
name = "cc"
|
||||
version = "1.0.50"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "95e28fa049fda1c330bcf9d723be7663a899c4679724b34c81e9f5a326aab8cd"
|
||||
|
||||
[[package]]
|
||||
name = "cfg-if"
|
||||
version = "0.1.10"
|
||||
@@ -84,10 +91,16 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "4785bdd1c96b2a846b2bd7cc02e86b6b3dbf14e7e53446c4f54c92a361040822"
|
||||
|
||||
[[package]]
|
||||
name = "clap"
|
||||
version = "2.33.0"
|
||||
name = "cfg-if"
|
||||
version = "1.0.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "5067f5bb2d80ef5d68b4c87db81601f0b75bca627bc2ef76b141d7b846a3c6d9"
|
||||
checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
|
||||
|
||||
[[package]]
|
||||
name = "clap"
|
||||
version = "2.33.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "37e58ac78573c40708d45522f0d80fa2f01cc4f9b4e2bf749807255454312002"
|
||||
dependencies = [
|
||||
"ansi_term",
|
||||
"atty",
|
||||
@@ -100,88 +113,96 @@ dependencies = [
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "cloudabi"
|
||||
version = "0.0.3"
|
||||
name = "const_fn"
|
||||
version = "0.4.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ddfc5b9aa5d4507acaf872de71051dfd0e309860e88966e1051e462a077aac4f"
|
||||
checksum = "ce90df4c658c62f12d78f7508cf92f9173e5184a539c10bfe54a3107b3ffd0f2"
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-channel"
|
||||
version = "0.5.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "dca26ee1f8d361640700bde38b2c37d8c22b3ce2d360e1fc1c74ea4b0aa7d775"
|
||||
dependencies = [
|
||||
"bitflags",
|
||||
"cfg-if 1.0.0",
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-deque"
|
||||
version = "0.7.3"
|
||||
version = "0.8.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "9f02af974daeee82218205558e51ec8768b48cf524bd01d550abe5573a608285"
|
||||
checksum = "94af6efb46fef72616855b036a624cf27ba656ffc9be1b9a3c931cfc7749a9a9"
|
||||
dependencies = [
|
||||
"cfg-if 1.0.0",
|
||||
"crossbeam-epoch",
|
||||
"crossbeam-utils",
|
||||
"maybe-uninit",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-epoch"
|
||||
version = "0.8.2"
|
||||
version = "0.9.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "058ed274caafc1f60c4997b5fc07bf7dc7cca454af7c6e81edffe5f33f70dace"
|
||||
checksum = "ec0f606a85340376eef0d6d8fec399e6d4a544d648386c6645eb6d0653b27d9f"
|
||||
dependencies = [
|
||||
"autocfg",
|
||||
"cfg-if",
|
||||
"cfg-if 1.0.0",
|
||||
"const_fn",
|
||||
"crossbeam-utils",
|
||||
"lazy_static",
|
||||
"maybe-uninit",
|
||||
"memoffset",
|
||||
"scopeguard",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-queue"
|
||||
version = "0.2.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "c695eeca1e7173472a32221542ae469b3e9aac3a4fc81f7696bcad82029493db"
|
||||
dependencies = [
|
||||
"cfg-if",
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-utils"
|
||||
version = "0.7.2"
|
||||
version = "0.8.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "c3c7c73a2d1e9fc0886a08b93e98eb643461230d5f1925e4036204d5f2e261a8"
|
||||
checksum = "ec91540d98355f690a86367e566ecad2e9e579f230230eb7c21398372be73ea5"
|
||||
dependencies = [
|
||||
"autocfg",
|
||||
"cfg-if",
|
||||
"cfg-if 1.0.0",
|
||||
"const_fn",
|
||||
"lazy_static",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "either"
|
||||
version = "1.5.3"
|
||||
version = "1.6.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "bb1f6b1ce1c140482ea30ddd3335fc0024ac7ee112895426e0a629a6c20adfe3"
|
||||
checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457"
|
||||
|
||||
[[package]]
|
||||
name = "error-chain"
|
||||
version = "0.12.2"
|
||||
version = "0.12.4"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d371106cc88ffdfb1eabd7111e432da544f16f3e2d7bf1dfe8bf575f1df045cd"
|
||||
checksum = "2d2f06b9cac1506ece98fe3231e3cc9c4410ec3d5b1f24ae1c8946f0742cdefc"
|
||||
dependencies = [
|
||||
"backtrace",
|
||||
"version_check",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "fuchsia-cprng"
|
||||
version = "0.1.1"
|
||||
name = "getrandom"
|
||||
version = "0.1.15"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "a06f77d526c1a601b7c4cdd98f54b5eaabffc14d5f2f0296febdc7f357c6d3ba"
|
||||
checksum = "fc587bc0ec293155d5bfa6b9891ec18a1e330c234f896ea47fbada4cadbe47e6"
|
||||
dependencies = [
|
||||
"cfg-if 0.1.10",
|
||||
"libc",
|
||||
"wasi",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "gimli"
|
||||
version = "0.22.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "aaf91faf136cb47367fa430cd46e37a788775e7fa104f8b4bcb3861dc389b724"
|
||||
|
||||
[[package]]
|
||||
name = "hermit-abi"
|
||||
version = "0.1.8"
|
||||
version = "0.1.17"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "1010591b26bbfe835e9faeabeb11866061cc7dcebffd56ad7d0942d0e61aefd8"
|
||||
checksum = "5aca5565f760fb5b220e499d72710ed156fdb74e631659e99377d9ebfbd13ae8"
|
||||
dependencies = [
|
||||
"libc",
|
||||
]
|
||||
@@ -194,75 +215,105 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
|
||||
|
||||
[[package]]
|
||||
name = "libc"
|
||||
version = "0.2.67"
|
||||
version = "0.2.79"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "eb147597cdf94ed43ab7a9038716637d2d1bf2bc571da995d0028dec06bd3018"
|
||||
|
||||
[[package]]
|
||||
name = "maybe-uninit"
|
||||
version = "2.0.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "60302e4db3a61da70c0cb7991976248362f30319e88850c487b9b95bbf059e00"
|
||||
checksum = "2448f6066e80e3bfc792e9c98bf705b4b0fc6e8ef5b43e5889aff0eaa9c58743"
|
||||
|
||||
[[package]]
|
||||
name = "memoffset"
|
||||
version = "0.5.3"
|
||||
version = "0.5.6"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "75189eb85871ea5c2e2c15abbdd541185f63b408415e5051f5cac122d8c774b9"
|
||||
checksum = "043175f069eda7b85febe4a74abbaeff828d9f8b448515d3151a14a3542811aa"
|
||||
dependencies = [
|
||||
"rustc_version",
|
||||
"autocfg",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "miniz_oxide"
|
||||
version = "0.4.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "0f2d26ec3309788e423cfbf68ad1800f061638098d76a83681af979dc4eda19d"
|
||||
dependencies = [
|
||||
"adler",
|
||||
"autocfg",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num_cpus"
|
||||
version = "1.12.0"
|
||||
version = "1.13.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "46203554f085ff89c235cd12f7075f3233af9b11ed7c9e16dfe2560d03313ce6"
|
||||
checksum = "05499f3756671c15885fee9034446956fff3f243d6077b91e5767df161f766b3"
|
||||
dependencies = [
|
||||
"hermit-abi",
|
||||
"libc",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "pkg-config"
|
||||
version = "0.3.17"
|
||||
name = "object"
|
||||
version = "0.21.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "05da548ad6865900e60eaba7f589cc0783590a92e940c26953ff81ddbab2d677"
|
||||
checksum = "37fd5004feb2ce328a52b0b3d01dbf4ffff72583493900ed15f22d4111c51693"
|
||||
|
||||
[[package]]
|
||||
name = "pkg-config"
|
||||
version = "0.3.19"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "3831453b3449ceb48b6d9c7ad7c96d5ea673e9b470a1dc578c2ce6521230884c"
|
||||
|
||||
[[package]]
|
||||
name = "ppv-lite86"
|
||||
version = "0.2.9"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "c36fa947111f5c62a733b652544dd0016a43ce89619538a8ef92724a6f501a20"
|
||||
|
||||
[[package]]
|
||||
name = "rand"
|
||||
version = "0.5.6"
|
||||
version = "0.7.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "c618c47cd3ebd209790115ab837de41425723956ad3ce2e6a7f09890947cacb9"
|
||||
checksum = "6a6b1679d49b24bbfe0c803429aa1874472f50d9b363131f0e89fc356b544d03"
|
||||
dependencies = [
|
||||
"cloudabi",
|
||||
"fuchsia-cprng",
|
||||
"getrandom",
|
||||
"libc",
|
||||
"rand_core 0.3.1",
|
||||
"winapi",
|
||||
"rand_chacha",
|
||||
"rand_core",
|
||||
"rand_hc",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand_core"
|
||||
version = "0.3.1"
|
||||
name = "rand_chacha"
|
||||
version = "0.2.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "7a6fdeb83b075e8266dcc8762c22776f6877a63111121f5f8c7411e5be7eed4b"
|
||||
checksum = "f4c8ed856279c9737206bf725bf36935d8666ead7aa69b52be55af369d193402"
|
||||
dependencies = [
|
||||
"rand_core 0.4.2",
|
||||
"ppv-lite86",
|
||||
"rand_core",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand_core"
|
||||
version = "0.4.2"
|
||||
version = "0.5.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "9c33a3c44ca05fa6f1807d8e6743f3824e8509beca625669633be0acbdf509dc"
|
||||
checksum = "90bde5296fc891b0cef12a6d03ddccc162ce7b2aff54160af9338f8d40df6d19"
|
||||
dependencies = [
|
||||
"getrandom",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand_hc"
|
||||
version = "0.2.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ca3129af7b92a17112d59ad498c6f81eaf463253766b90396d39ea7a39d6613c"
|
||||
dependencies = [
|
||||
"rand_core",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rayon"
|
||||
version = "1.3.0"
|
||||
version = "1.5.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "db6ce3297f9c85e16621bb8cca38a06779ffc31bb8184e1be4bed2be4678a098"
|
||||
checksum = "8b0d8e0819fadc20c74ea8373106ead0600e3a67ef1fe8da56e39b9ae7275674"
|
||||
dependencies = [
|
||||
"autocfg",
|
||||
"crossbeam-deque",
|
||||
"either",
|
||||
"rayon-core",
|
||||
@@ -270,12 +321,12 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "rayon-core"
|
||||
version = "1.7.0"
|
||||
version = "1.9.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "08a89b46efaf957e52b18062fb2f4660f8b8a4dde1807ca002690868ef2c85a9"
|
||||
checksum = "9ab346ac5921dc62ffa9f89b7a773907511cdfa5490c572ae9be1be33e8afa4a"
|
||||
dependencies = [
|
||||
"crossbeam-channel",
|
||||
"crossbeam-deque",
|
||||
"crossbeam-queue",
|
||||
"crossbeam-utils",
|
||||
"lazy_static",
|
||||
"num_cpus",
|
||||
@@ -283,18 +334,9 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "rustc-demangle"
|
||||
version = "0.1.16"
|
||||
version = "0.1.18"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "4c691c0e608126e00913e33f0ccf3727d5fc84573623b8d65b2df340b5201783"
|
||||
|
||||
[[package]]
|
||||
name = "rustc_version"
|
||||
version = "0.2.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "138e3e0acb6c9fb258b19b67cb8abd63c00679d2851805ea151465464fe9030a"
|
||||
dependencies = [
|
||||
"semver",
|
||||
]
|
||||
checksum = "6e3bad0ee36814ca07d7968269dd4b7ec89ec2da10c4bb613928d3077083c232"
|
||||
|
||||
[[package]]
|
||||
name = "scopeguard"
|
||||
@@ -302,21 +344,6 @@ version = "1.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd"
|
||||
|
||||
[[package]]
|
||||
name = "semver"
|
||||
version = "0.9.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "1d7eb9ef2c18661902cc47e535f9bc51b78acd254da71d375c2f6720d9a40403"
|
||||
dependencies = [
|
||||
"semver-parser",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "semver-parser"
|
||||
version = "0.7.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "388a1df253eca08550bef6c72392cfe7c30914bf41df5269b68cbd6ff8f570a3"
|
||||
|
||||
[[package]]
|
||||
name = "strsim"
|
||||
version = "0.8.0"
|
||||
@@ -334,27 +361,34 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "unicode-width"
|
||||
version = "0.1.7"
|
||||
version = "0.1.8"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "caaa9d531767d1ff2150b9332433f32a24622147e5ebb1f26409d5da67afd479"
|
||||
checksum = "9337591893a19b88d8d87f2cec1e73fad5cdfd10e5a6f349f498ad6ea2ffb1e3"
|
||||
|
||||
[[package]]
|
||||
name = "vec_map"
|
||||
version = "0.8.1"
|
||||
version = "0.8.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "05c78687fb1a80548ae3250346c3db86a80a7cdd77bda190189f2d0a0987c81a"
|
||||
checksum = "f1bddf1187be692e79c5ffeab891132dfb0f236ed36a43c7ed39f1165ee20191"
|
||||
|
||||
[[package]]
|
||||
name = "version_check"
|
||||
version = "0.9.1"
|
||||
version = "0.9.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "078775d0255232fb988e6fccf26ddc9d1ac274299aaedcedce21c6f72cc533ce"
|
||||
checksum = "b5a972e5669d67ba988ce3dc826706fb0a8b01471c088cb0b6110b805cc36aed"
|
||||
|
||||
[[package]]
|
||||
name = "wasi"
|
||||
version = "0.9.0+wasi-snapshot-preview1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519"
|
||||
|
||||
[[package]]
|
||||
name = "wham"
|
||||
version = "0.9.6"
|
||||
version = "0.9.9"
|
||||
dependencies = [
|
||||
"GSL",
|
||||
"assert_approx_eq",
|
||||
"clap",
|
||||
"error-chain",
|
||||
"rand",
|
||||
@@ -363,9 +397,9 @@ dependencies = [
|
||||
|
||||
[[package]]
|
||||
name = "winapi"
|
||||
version = "0.3.8"
|
||||
version = "0.3.9"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "8093091eeb260906a183e6ae1abdba2ef5ef2257a21801128899c3fc699229c6"
|
||||
checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419"
|
||||
dependencies = [
|
||||
"winapi-i686-pc-windows-gnu",
|
||||
"winapi-x86_64-pc-windows-gnu",
|
||||
|
||||
18
Cargo.toml
18
Cargo.toml
@@ -1,20 +1,24 @@
|
||||
[package]
|
||||
name = "wham"
|
||||
version = "0.9.6"
|
||||
version = "0.9.9"
|
||||
authors = ["Daniel Bauer <bauer@cbs.tu-darmstadt.de>"]
|
||||
description = "An implementation of the weighted histogram analysis method"
|
||||
license = "GPL-3.0"
|
||||
repository = "https://github.com/danijoo/WHAM"
|
||||
readme = "README.md"
|
||||
categories = ["science", "command-line-utilities", "molecular-dynamics", "algorithms"]
|
||||
keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dynamics"]
|
||||
|
||||
[dependencies]
|
||||
clap = {version="2.32.0", features=['yaml']}
|
||||
error-chain = "0.12.0"
|
||||
rand = "0.5.5"
|
||||
GSL = "1.0.0"
|
||||
rand = "0.7.*"
|
||||
GSL = "1.1"
|
||||
rayon = "1.0.3"
|
||||
|
||||
[profile.release]
|
||||
opt-level = 2
|
||||
[dev-dependencies]
|
||||
assert_approx_eq = "1.1.0"
|
||||
|
||||
[profile.release]
|
||||
opt-level = 3
|
||||
|
||||
[features]
|
||||
default = ["GSL/v2"]
|
||||
|
||||
55
README.md
55
README.md
@@ -38,6 +38,56 @@ Usage
|
||||
wham has a convenient command line interface. You can see all options with
|
||||
```wham -h```:
|
||||
|
||||
```
|
||||
wham 0.9.9
|
||||
D. Bauer <bauer@bio.tu-darmstadt.de>
|
||||
wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports
|
||||
potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
|
||||
|
||||
Metadata file format:
|
||||
/path/to/timeseries_file1 x_1 x_2 x_N fc_1 fc_2 fc_N
|
||||
/path/to/timeseries_file2 x_1 x_2 x_N fc_1 fc_2 fc_N
|
||||
/path/to/timeseries_file3 x_1 x_2 x_N fc_1 fc_2 fc_N
|
||||
The first column is a path to a timeseries file _relative_ to the metadata file (see below). This is followed by the
|
||||
position of the umbrella potential x in N dimensions and the force constant fc in each dimension. Lines starting with a
|
||||
# are treated as comments and will not be parsed.
|
||||
|
||||
Timeseries file format:
|
||||
time x_1 x_2 x_N
|
||||
time x_1 x_2 x_N
|
||||
time x_1 x_2 x_N
|
||||
The first column will be ignored and is followed by N reaction coordinates x.
|
||||
|
||||
Shipped under the GPLv3 license.
|
||||
|
||||
USAGE:
|
||||
wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature>
|
||||
|
||||
FLAGS:
|
||||
-c, --cyclic For periodic reaction coordinates. If this is set, the first and last coordinate bin in each
|
||||
dimension are treated as neighbors for the bias calculation.
|
||||
-h, --help Prints help information
|
||||
-V, --version Prints version information
|
||||
-v, --verbose Enables verbose output.
|
||||
|
||||
OPTIONS:
|
||||
-b, --bins <BINS> Number of histogram bins (comma separated).
|
||||
--bt <bootstrap> Number of bayesian bootstrapping runs for error analysis by assigning random
|
||||
weights (defaults to 0).
|
||||
--seed <bootstrap_seed> Random seed for bootstrapping runs.
|
||||
--end <end> Skip rows in timeseries with an index larger than this value (defaults to 1e+20)
|
||||
-i, --iterations <ITERATIONS> Stop WHAM after this many iterations without convergence (defaults to 100,000).
|
||||
--max <HIST_MAX> Histogram maxima (comma separated). Also accepts "pi".
|
||||
-f, --file <METADATA> Path to the metadata file.
|
||||
--min <HIST_MIN> Histogram minima (comma separated for multiple dimensions). Also accepts "pi".
|
||||
-o, --output <output> Free energy output file (defaults to wham.out).
|
||||
--start <start> Skip rows in timeseries with an index smaller than this value (defaults to 0)
|
||||
-T, --temperature <temperature> WHAM temperature in Kelvin.
|
||||
-t, --tolerance <TOLERANCE> Abortion criteria for WHAM calculation. WHAM stops if abs(F_new - F_old) <
|
||||
tolerance (defaults to 0.000001).
|
||||
|
||||
```
|
||||
|
||||
To run the two dimensional example (simulation of dialanine phi and psi angle):
|
||||
```bash
|
||||
wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example/2d/metadata.dat
|
||||
@@ -81,8 +131,8 @@ one can estimate the error as standard deviation between the N bootstrapping run
|
||||
Autocorrelation Estimates, JCTC, 6(12), 3713-3720*.
|
||||
|
||||
To perform bayesian bootstrapping in WHAM, use the ```-bt <RUNS>``` flag to perform <RUNS> individual bootstrapping
|
||||
runs. The error estimates of bin probabilities and free energy will be given as separate column (+/-) in the output file.
|
||||
If no error analysis is performed, these columns are set to 0.0.
|
||||
runs. The error estimates of bin probabilities and free energy will be given as standard error (SE) in a
|
||||
separate column (+/-) in the output file. If no error analysis is performed, these columns are set to 0.0.
|
||||
|
||||
Examples
|
||||
---
|
||||
@@ -94,7 +144,6 @@ The example folder contains input and output files for two simple test systems:
|
||||
|
||||
TODO
|
||||
---
|
||||
- better error estimation
|
||||
- Autocorrelation
|
||||
- Replica exchange
|
||||
|
||||
|
||||
@@ -1,101 +1,101 @@
|
||||
#coord1 Free Energy +/- Probability +/-
|
||||
-3.110177 7.158102 0.476037 0.003494 0.000667
|
||||
-3.047345 5.365727 0.470738 0.007168 0.001353
|
||||
-2.984513 3.873190 0.433120 0.013039 0.002264
|
||||
-2.921681 2.953162 0.441412 0.018855 0.003337
|
||||
-2.858849 1.949554 0.443549 0.028195 0.005014
|
||||
-2.796017 1.391747 0.407397 0.035261 0.005759
|
||||
-2.733186 1.128270 0.380789 0.039189 0.005983
|
||||
-2.670354 0.839970 0.370008 0.043991 0.006526
|
||||
-2.607522 0.624769 0.381081 0.047955 0.007326
|
||||
-2.544690 0.663757 0.389082 0.047211 0.007364
|
||||
-2.481858 1.051932 0.379999 0.040407 0.006156
|
||||
-2.419026 1.463048 0.394914 0.034267 0.005425
|
||||
-2.356194 1.990616 0.335141 0.027734 0.003726
|
||||
-2.293363 2.190692 0.196497 0.025597 0.002016
|
||||
-2.230531 2.553036 0.177045 0.022136 0.001571
|
||||
-2.167699 2.572522 0.208798 0.021964 0.001839
|
||||
-2.104867 2.472360 0.168929 0.022863 0.001548
|
||||
-2.042035 2.517562 0.166774 0.022453 0.001501
|
||||
-1.979203 2.469778 0.182308 0.022887 0.001673
|
||||
-1.916372 2.223125 0.405294 0.025266 0.004105
|
||||
-1.853540 2.080157 0.340376 0.026756 0.003651
|
||||
-1.790708 1.793841 0.394974 0.030011 0.004752
|
||||
-1.727876 1.458784 0.335956 0.034326 0.004623
|
||||
-1.665044 0.914441 0.336724 0.042697 0.005764
|
||||
-1.602212 0.303813 0.320695 0.054540 0.007012
|
||||
-1.539380 0.268335 0.330495 0.055321 0.007330
|
||||
-1.476549 0.000000 0.338143 0.061604 0.008351
|
||||
-1.413717 0.537061 0.369220 0.049671 0.007352
|
||||
-1.350885 1.439940 0.381006 0.034586 0.005283
|
||||
-1.288053 2.391838 0.407306 0.023614 0.003856
|
||||
-1.225221 3.779470 0.413629 0.013538 0.002245
|
||||
-1.162389 5.685555 0.432485 0.006305 0.001093
|
||||
-1.099557 7.661896 0.446994 0.002855 0.000512
|
||||
-1.036726 9.946594 0.489990 0.001142 0.000224
|
||||
-0.973894 12.359408 0.544847 0.000434 0.000095
|
||||
-0.911062 14.954655 0.622373 0.000153 0.000038
|
||||
-0.848230 17.742242 0.747233 0.000050 0.000015
|
||||
-0.785398 20.555785 0.921773 0.000016 0.000006
|
||||
-0.722566 22.811160 1.179343 0.000007 0.000003
|
||||
-0.659734 25.178094 1.484376 0.000003 0.000002
|
||||
-0.596903 26.442288 1.565063 0.000002 0.000001
|
||||
-0.534071 27.897565 1.418232 0.000001 0.000000
|
||||
-0.471239 29.062473 1.246738 0.000001 0.000000
|
||||
-0.408407 30.384521 1.231662 0.000000 0.000000
|
||||
-0.345575 31.638454 1.209295 0.000000 0.000000
|
||||
-0.282743 32.817727 1.291637 0.000000 0.000000
|
||||
-0.219911 33.770369 1.323362 0.000000 0.000000
|
||||
-0.157080 34.505503 1.396454 0.000000 0.000000
|
||||
-0.094248 35.431659 1.662548 0.000000 0.000000
|
||||
-0.031416 35.615810 1.766523 0.000000 0.000000
|
||||
0.031416 35.561946 1.799603 0.000000 0.000000
|
||||
0.094248 35.382089 1.374764 0.000000 0.000000
|
||||
0.157080 34.934827 1.659669 0.000000 0.000000
|
||||
0.219911 33.673460 1.317449 0.000000 0.000000
|
||||
0.282743 32.731563 1.191975 0.000000 0.000000
|
||||
0.345575 31.261855 1.199834 0.000000 0.000000
|
||||
0.408407 29.717377 1.339899 0.000000 0.000000
|
||||
0.471239 28.076620 1.303458 0.000001 0.000000
|
||||
0.534071 26.481097 1.281459 0.000002 0.000001
|
||||
0.596903 24.487358 1.251139 0.000003 0.000002
|
||||
0.659734 22.344251 1.233967 0.000008 0.000004
|
||||
0.722566 20.241543 1.251834 0.000018 0.000009
|
||||
0.785398 18.341869 1.241684 0.000039 0.000020
|
||||
0.848230 16.261582 1.281159 0.000091 0.000047
|
||||
0.911062 14.301801 1.294427 0.000199 0.000103
|
||||
0.973894 12.603788 1.270069 0.000394 0.000200
|
||||
1.036726 11.249601 1.273360 0.000678 0.000346
|
||||
1.099557 10.087886 1.280795 0.001079 0.000554
|
||||
1.162389 9.443303 1.266983 0.001398 0.000710
|
||||
1.225221 9.152799 1.255009 0.001570 0.000790
|
||||
1.288053 9.331937 1.256822 0.001462 0.000736
|
||||
1.350885 9.905546 1.251237 0.001161 0.000583
|
||||
1.413717 11.042050 1.263816 0.000736 0.000373
|
||||
1.476549 12.598167 1.275325 0.000395 0.000202
|
||||
1.539380 14.520167 1.280735 0.000183 0.000094
|
||||
1.602212 16.569783 1.295384 0.000080 0.000042
|
||||
1.665044 18.687390 1.309518 0.000034 0.000018
|
||||
1.727876 20.775408 1.356414 0.000015 0.000008
|
||||
1.790708 22.905200 1.308568 0.000006 0.000003
|
||||
1.853540 24.643852 1.425870 0.000003 0.000002
|
||||
1.916372 26.301740 1.362138 0.000002 0.000001
|
||||
1.979203 27.372071 1.391716 0.000001 0.000001
|
||||
2.042035 28.697726 1.504676 0.000001 0.000000
|
||||
2.104867 29.417901 1.369019 0.000000 0.000000
|
||||
2.167699 30.008351 1.409081 0.000000 0.000000
|
||||
2.230531 30.406016 1.240128 0.000000 0.000000
|
||||
2.293363 30.171275 1.179444 0.000000 0.000000
|
||||
2.356194 29.884646 1.219365 0.000000 0.000000
|
||||
2.419026 29.428153 1.177281 0.000000 0.000000
|
||||
2.481858 28.546982 1.351057 0.000001 0.000000
|
||||
2.544690 27.757520 1.238695 0.000001 0.000000
|
||||
2.607522 26.505787 1.308430 0.000001 0.000001
|
||||
2.670354 24.491866 1.039660 0.000003 0.000001
|
||||
2.733186 22.320664 0.927893 0.000008 0.000003
|
||||
2.796017 20.052723 0.894483 0.000020 0.000007
|
||||
2.858849 17.655650 0.893580 0.000052 0.000019
|
||||
2.921681 15.471590 0.930897 0.000125 0.000047
|
||||
2.984513 13.138167 0.837026 0.000318 0.000107
|
||||
3.047345 11.092386 0.674785 0.000722 0.000195
|
||||
3.110177 9.065722 0.521230 0.001626 0.000340
|
||||
-3.110177 7.158102 0.071438 0.003494 0.000068
|
||||
-3.047345 5.365727 0.069449 0.007168 0.000135
|
||||
-2.984513 3.873190 0.067495 0.013039 0.000232
|
||||
-2.921681 2.953162 0.067589 0.018855 0.000343
|
||||
-2.858849 1.949554 0.064640 0.028195 0.000480
|
||||
-2.796017 1.391747 0.063570 0.035261 0.000584
|
||||
-2.733186 1.128270 0.061710 0.039189 0.000620
|
||||
-2.670354 0.839970 0.060445 0.043991 0.000667
|
||||
-2.607522 0.624769 0.060762 0.047955 0.000739
|
||||
-2.544690 0.663757 0.060786 0.047211 0.000731
|
||||
-2.481858 1.051932 0.059774 0.040407 0.000617
|
||||
-2.419026 1.463048 0.060422 0.034267 0.000527
|
||||
-2.356194 1.990616 0.055242 0.027734 0.000360
|
||||
-2.293363 2.190692 0.044458 0.025597 0.000210
|
||||
-2.230531 2.553036 0.042144 0.022136 0.000159
|
||||
-2.167699 2.572522 0.043313 0.021964 0.000165
|
||||
-2.104867 2.472360 0.039799 0.022863 0.000157
|
||||
-2.042035 2.517562 0.036792 0.022453 0.000203
|
||||
-1.979203 2.469778 0.036324 0.022887 0.000237
|
||||
-1.916372 2.223125 0.049496 0.025266 0.000565
|
||||
-1.853540 2.080157 0.039205 0.026756 0.000453
|
||||
-1.790708 1.793841 0.041488 0.030011 0.000558
|
||||
-1.727876 1.458784 0.042362 0.034326 0.000565
|
||||
-1.665044 0.914441 0.037561 0.042697 0.000682
|
||||
-1.602212 0.303813 0.028701 0.054540 0.000761
|
||||
-1.539380 0.268335 0.025765 0.055321 0.000771
|
||||
-1.476549 0.000000 0.023371 0.061604 0.000866
|
||||
-1.413717 0.537061 0.026370 0.049671 0.000715
|
||||
-1.350885 1.439940 0.025620 0.034586 0.000547
|
||||
-1.288053 2.391838 0.027653 0.023614 0.000393
|
||||
-1.225221 3.779470 0.027093 0.013538 0.000229
|
||||
-1.162389 5.685555 0.027762 0.006305 0.000111
|
||||
-1.099557 7.661896 0.028996 0.002855 0.000051
|
||||
-1.036726 9.946594 0.032118 0.001142 0.000021
|
||||
-0.973894 12.359408 0.042801 0.000434 0.000009
|
||||
-0.911062 14.954655 0.052443 0.000153 0.000004
|
||||
-0.848230 17.742242 0.063261 0.000050 0.000002
|
||||
-0.785398 20.555785 0.068133 0.000016 0.000001
|
||||
-0.722566 22.811160 0.075518 0.000007 0.000000
|
||||
-0.659734 25.178094 0.086226 0.000003 0.000000
|
||||
-0.596903 26.442288 0.091783 0.000002 0.000000
|
||||
-0.534071 27.897565 0.093175 0.000001 0.000000
|
||||
-0.471239 29.062473 0.096027 0.000001 0.000000
|
||||
-0.408407 30.384521 0.096855 0.000000 0.000000
|
||||
-0.345575 31.638454 0.099956 0.000000 0.000000
|
||||
-0.282743 32.817727 0.107689 0.000000 0.000000
|
||||
-0.219911 33.770369 0.108946 0.000000 0.000000
|
||||
-0.157080 34.505503 0.112902 0.000000 0.000000
|
||||
-0.094248 35.431659 0.120856 0.000000 0.000000
|
||||
-0.031416 35.615810 0.122046 0.000000 0.000000
|
||||
0.031416 35.561946 0.118225 0.000000 0.000000
|
||||
0.094248 35.382089 0.108652 0.000000 0.000000
|
||||
0.157080 34.934827 0.118384 0.000000 0.000000
|
||||
0.219911 33.673460 0.119195 0.000000 0.000000
|
||||
0.282743 32.731563 0.121592 0.000000 0.000000
|
||||
0.345575 31.261855 0.129698 0.000000 0.000000
|
||||
0.408407 29.717377 0.139789 0.000000 0.000000
|
||||
0.471239 28.076620 0.137289 0.000001 0.000000
|
||||
0.534071 26.481097 0.135712 0.000002 0.000000
|
||||
0.596903 24.487358 0.135156 0.000003 0.000000
|
||||
0.659734 22.344251 0.131027 0.000008 0.000000
|
||||
0.722566 20.241543 0.133145 0.000018 0.000001
|
||||
0.785398 18.341869 0.131952 0.000039 0.000002
|
||||
0.848230 16.261582 0.134676 0.000091 0.000006
|
||||
0.911062 14.301801 0.134417 0.000199 0.000013
|
||||
0.973894 12.603788 0.131018 0.000394 0.000024
|
||||
1.036726 11.249601 0.131276 0.000678 0.000043
|
||||
1.099557 10.087886 0.132498 0.001079 0.000070
|
||||
1.162389 9.443303 0.132990 0.001398 0.000092
|
||||
1.225221 9.152799 0.132250 0.001570 0.000104
|
||||
1.288053 9.331937 0.133178 0.001462 0.000099
|
||||
1.350885 9.905546 0.133357 0.001161 0.000078
|
||||
1.413717 11.042050 0.133807 0.000736 0.000051
|
||||
1.476549 12.598167 0.132722 0.000395 0.000027
|
||||
1.539380 14.520167 0.131816 0.000183 0.000012
|
||||
1.602212 16.569783 0.131773 0.000080 0.000005
|
||||
1.665044 18.687390 0.132601 0.000034 0.000002
|
||||
1.727876 20.775408 0.133867 0.000015 0.000001
|
||||
1.790708 22.905200 0.129265 0.000006 0.000000
|
||||
1.853540 24.643852 0.128894 0.000003 0.000000
|
||||
1.916372 26.301740 0.130266 0.000002 0.000000
|
||||
1.979203 27.372071 0.128620 0.000001 0.000000
|
||||
2.042035 28.697726 0.133263 0.000001 0.000000
|
||||
2.104867 29.417901 0.133513 0.000000 0.000000
|
||||
2.167699 30.008351 0.130925 0.000000 0.000000
|
||||
2.230531 30.406016 0.124398 0.000000 0.000000
|
||||
2.293363 30.171275 0.122493 0.000000 0.000000
|
||||
2.356194 29.884646 0.130278 0.000000 0.000000
|
||||
2.419026 29.428153 0.130844 0.000000 0.000000
|
||||
2.481858 28.546982 0.148114 0.000001 0.000000
|
||||
2.544690 27.757520 0.133045 0.000001 0.000000
|
||||
2.607522 26.505787 0.134364 0.000001 0.000000
|
||||
2.670354 24.491866 0.115061 0.000003 0.000000
|
||||
2.733186 22.320664 0.110036 0.000008 0.000000
|
||||
2.796017 20.052723 0.107894 0.000020 0.000001
|
||||
2.858849 17.655650 0.105964 0.000052 0.000002
|
||||
2.921681 15.471590 0.107005 0.000125 0.000005
|
||||
2.984513 13.138167 0.099378 0.000318 0.000012
|
||||
3.047345 11.092386 0.087770 0.000722 0.000022
|
||||
3.110177 9.065722 0.077092 0.001626 0.000038
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
name: wham
|
||||
version: "0.9.0"
|
||||
version: "0.9.9"
|
||||
author: D. Bauer <bauer@bio.tu-darmstadt.de>
|
||||
about: |
|
||||
wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
|
||||
|
||||
@@ -1,7 +1,7 @@
|
||||
use rand::{SeedableRng, StdRng, Rng};
|
||||
use rand::prelude::*;
|
||||
use super::histogram::{Dataset};
|
||||
use super::perform_wham;
|
||||
use super::Config;
|
||||
use super::{Config,calc_free_energy};
|
||||
use rgsl::statistics;
|
||||
|
||||
// returns a set of num_windows continious weights by
|
||||
@@ -20,7 +20,7 @@ fn generate_random_weights(num_windows: usize, rng: &mut StdRng) -> Vec<f64> {
|
||||
for i in 0..num_windows {
|
||||
weights[i] = rnds[i+1] - rnds[i]
|
||||
}
|
||||
return weights
|
||||
weights
|
||||
}
|
||||
|
||||
// Generate a random weighted dataset from the given dataset by changing the weights
|
||||
@@ -32,7 +32,7 @@ fn generate_random_weighted_dataset(ds: Dataset, rng: &mut StdRng) -> Dataset {
|
||||
// Perform bootstrap error analysis. This runs the WHAM analysis num_runs times on random weighted
|
||||
// datasets. The standard deviation is calculated on the bootstrapped probabilities of each bin. The
|
||||
// standard deviation of the free eneergy is then deduced by error propagation (A_std = kT*1/P*P_std)
|
||||
pub fn run_bootstrap(cfg: &Config, ds: Dataset, P: &[f64], num_runs: usize) -> (Vec<f64>,Vec<f64>) {
|
||||
pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec<f64>,Vec<f64>) {
|
||||
// seed the rng
|
||||
let mut rng: StdRng = SeedableRng::seed_from_u64(cfg.bootstrap_seed);
|
||||
|
||||
@@ -43,16 +43,27 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, P: &[f64], num_runs: usize) -> (
|
||||
perform_wham(cfg, &rnd_weighted_dataset).unwrap().0
|
||||
}).collect();
|
||||
|
||||
// Evaulate standard deviation of P per bin
|
||||
let mut P_std = vec![0.0; ds.num_bins];
|
||||
// Standard error (SE) of P per bin
|
||||
// SE = SD/sqrt(n)
|
||||
let mut P_se = vec![0.0; ds.num_bins];
|
||||
for bin in 0..ds.num_bins {
|
||||
let Ps = bootstrapped_Ps.iter().map(|window| window[bin]).collect::<Vec<f64>>();
|
||||
P_std[bin] = statistics::sd(&Ps, 1, num_runs);
|
||||
P_se[bin] = statistics::sd(&Ps, 1, num_runs)/(num_runs as f64).sqrt();
|
||||
}
|
||||
|
||||
// A_std by error propagation
|
||||
let A_std = P_std.iter().zip(P.iter()).map(|(std,P)| ds.kT*1.0/P*std).collect();
|
||||
(P_std, A_std)
|
||||
// SE of A
|
||||
let bootstrapped_As: Vec<Vec<f64>> = (0..num_runs).map(|x| {
|
||||
let run_Ps = &bootstrapped_Ps[x];
|
||||
calc_free_energy(&ds, run_Ps)
|
||||
}).collect();
|
||||
|
||||
let mut A_se = vec![0.0; ds.num_bins];
|
||||
for bin in 0..ds.num_bins {
|
||||
let As = bootstrapped_As.iter().map(|window| window[bin]).collect::<Vec<f64>>();
|
||||
A_se[bin] = statistics::sd(&As, 1, num_runs)/(num_runs as f64).sqrt();
|
||||
}
|
||||
|
||||
(P_se, A_se)
|
||||
}
|
||||
|
||||
#[cfg(tests)]
|
||||
|
||||
@@ -61,7 +61,9 @@ pub struct Dataset {
|
||||
|
||||
impl Dataset {
|
||||
|
||||
pub fn new(num_bins: usize, dimens_lengths: Vec<usize>, bin_width: Vec<f64>, hist_min: Vec<f64>, hist_max: Vec<f64>, bias_pos: Vec<f64>, bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset {
|
||||
pub fn new(num_bins: usize, dimens_lengths: Vec<usize>, bin_width: Vec<f64>,
|
||||
hist_min: Vec<f64>, hist_max: Vec<f64>, bias_pos: Vec<f64>,
|
||||
bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset {
|
||||
let num_windows = histograms.len();
|
||||
let bias: Vec<f64> = vec![0.0; num_bins*num_windows];
|
||||
let weights = vec![1.0; num_windows];
|
||||
@@ -92,7 +94,7 @@ impl Dataset {
|
||||
|
||||
pub fn new_weighted(ds: Dataset, weights: Vec<f64>) -> Dataset {
|
||||
Dataset {
|
||||
weights: weights,
|
||||
weights,
|
||||
..ds
|
||||
}
|
||||
}
|
||||
@@ -105,9 +107,9 @@ impl Dataset {
|
||||
let mut tmp = bin;
|
||||
let mut idx = vec![0; lengths.len()];
|
||||
for dimen in (1..lengths.len()).rev() {
|
||||
let denom = lengths.iter().take(dimen).fold(1, |s,&x| s*x);
|
||||
let denom: usize = lengths.iter().take(dimen).product();
|
||||
idx[dimen] = tmp / denom;
|
||||
tmp = tmp % denom;
|
||||
tmp %= denom;
|
||||
}
|
||||
idx[0] = tmp;
|
||||
idx
|
||||
@@ -151,8 +153,7 @@ impl Dataset {
|
||||
// store exp(U/kT) for better performance
|
||||
bias_sum += 0.5 * bias_fc[i] * dist * dist
|
||||
}
|
||||
let bias_sum = (-bias_sum/self.kT).exp();
|
||||
bias_sum
|
||||
(-bias_sum/self.kT).exp()
|
||||
}
|
||||
}
|
||||
|
||||
@@ -205,13 +206,13 @@ mod tests {
|
||||
let ds = build_hist_set(); // k = 10
|
||||
|
||||
// 3th element -> x=3.5, x0=3.5
|
||||
assert_delta!(0.134722337796, ds.calc_bias(3, 0), 0.00000001);
|
||||
assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01);
|
||||
|
||||
// 8th element -> x=8.5, x0=3.5
|
||||
assert_delta!(1.0, ds.calc_bias(4,0), 0.00000001);
|
||||
assert_delta!(1.0, ds.calc_bias(4,0), 0.000_000_01);
|
||||
|
||||
// 1st element -> x=0.5, x0=3.5. non-cyclic!
|
||||
assert_delta!(0.0, ds.calc_bias(0,0), 0.0000001);
|
||||
assert_delta!(0.0, ds.calc_bias(0,0), 0.000_000_1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -220,28 +221,28 @@ mod tests {
|
||||
ds.cyclic = true;
|
||||
|
||||
// 7th element -> x=3.5, x0=3.5
|
||||
assert_delta!(0.134722337796, ds.calc_bias(3, 0), 0.00000001);
|
||||
assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01);
|
||||
|
||||
// 8th element -> x=4.5, x0=3.5
|
||||
assert_delta!(1.0, ds.calc_bias(4, 0), 0.00000001);
|
||||
assert_delta!(1.0, ds.calc_bias(4, 0), 0.000_000_01);
|
||||
|
||||
|
||||
// 1th element -> x=0.5, x0=3.5
|
||||
// cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2
|
||||
assert_delta!(0.0000000000000117769, ds.calc_bias(0, 0), 0.00000001);
|
||||
assert_delta!(0.000_000_000_000_011_776_9, ds.calc_bias(0, 0), 0.000_000_01);
|
||||
|
||||
// 2nd element -> x=1.5, x0=3.5
|
||||
assert_delta!(0.00000001, ds.calc_bias(1, 0), 0.00000001);
|
||||
assert_delta!(0.000_000_01, ds.calc_bias(1, 0), 0.000_000_01);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn get_x_for_bin() {
|
||||
let ds = build_hist_set();
|
||||
let expected: Vec<f64> = vec![0,1,2,3,4,5,6,7,8].iter()
|
||||
.map(|x| *x as f64 + 0.5).collect();
|
||||
for i in 0..9 {
|
||||
assert_eq!(expected[i], ds.get_coords_for_bin(i)[0]);
|
||||
}
|
||||
.map(|x| *x as f64 + 0.5).collect();
|
||||
expected.iter().enumerate().for_each(|(i, exp)| {
|
||||
assert_approx_eq!(exp, &ds.get_coords_for_bin(i)[0]);
|
||||
})
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -258,10 +259,10 @@ mod tests {
|
||||
vec![build_hist(), build_hist()], // hists
|
||||
false // cyclic
|
||||
);
|
||||
assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.0000000001);
|
||||
assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.0000000001);
|
||||
assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.0000000001);
|
||||
assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.0000000001);
|
||||
assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.0000000001);
|
||||
assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.000_000_000_1);
|
||||
assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.000_000_000_1);
|
||||
assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.000_000_000_1);
|
||||
assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.000_000_000_1);
|
||||
assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.000_000_000_1);
|
||||
}
|
||||
}
|
||||
57
src/io.rs
57
src/io.rs
@@ -35,7 +35,7 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
|
||||
let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
|
||||
(cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
|
||||
}).collect();
|
||||
let num_bins = cfg.num_bins.iter().fold(1, |state, &bins| state*bins);
|
||||
let num_bins = cfg.num_bins.iter().product();
|
||||
let dimens_length = cfg.num_bins.clone();
|
||||
|
||||
let f = File::open(&cfg.metadata_file).chain_err(|| "Failed to open metadata file")?;
|
||||
@@ -46,7 +46,7 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
|
||||
let line = l.chain_err(|| "Failed to read line")?;
|
||||
|
||||
// skip comments and empty lines
|
||||
if line.starts_with("#") || line.len() == 0 {
|
||||
if line.starts_with('#') || line.is_empty() {
|
||||
continue;
|
||||
}
|
||||
|
||||
@@ -67,19 +67,19 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
|
||||
histograms.last().unwrap().num_points), cfg.verbose);
|
||||
|
||||
// parse bias force constants and positions
|
||||
for i in 1..cfg.dimens+1 {
|
||||
let pos = split[i].parse()
|
||||
for val in split.iter().skip(1).take(cfg.dimens) {
|
||||
let pos = val.parse()
|
||||
.chain_err(|| format!("Failed to read bias position in line {} of metadata file", line_num+1))?;
|
||||
bias_pos.push(pos);
|
||||
}
|
||||
for i in (1+cfg.dimens)..(1+2*cfg.dimens) {
|
||||
let fc = split[i].parse()
|
||||
for val in split.iter().skip(1+cfg.dimens).take(cfg.dimens) {
|
||||
let fc = val.parse()
|
||||
.chain_err(|| format!("Failed to read bias fc in line {} of metadata file", line_num+1))?;
|
||||
bias_fc.push(fc);
|
||||
}
|
||||
}
|
||||
|
||||
if histograms.len() > 0 {
|
||||
if !histograms.is_empty() {
|
||||
Ok(Dataset::new(num_bins, dimens_length, bin_width, cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos, bias_fc, kT, histograms, cfg.cyclic))
|
||||
} else {
|
||||
bail!("Histogram has no datapoints.")
|
||||
@@ -91,13 +91,10 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
|
||||
// lengths: length of the matrix in each dimension
|
||||
// returns an index if the matrix is flattened to a one dimensional vector
|
||||
// example for 3 dimensions N,M,O: idx = i_O + l_O*l_M*i_M + l_O*l_M*l_N*i_N
|
||||
fn flat_index(indeces: &Vec<usize>, lengths: &Vec<usize>) -> usize {
|
||||
let mut idx = 0;
|
||||
for i in 0..indeces.len() {
|
||||
idx += indeces[i]*lengths[0..i].iter()
|
||||
.fold(1, |state, &l| { state * l });
|
||||
}
|
||||
idx
|
||||
fn flat_index(indeces: &[usize], lengths: &[usize]) -> usize {
|
||||
indeces.iter().enumerate().map(|(i, idx)| {
|
||||
idx * lengths.iter().take(i).product::<usize>()
|
||||
}).sum()
|
||||
}
|
||||
|
||||
// returns true if the values are inside the histogram boundaries defined by cfg
|
||||
@@ -118,14 +115,14 @@ fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool {
|
||||
false
|
||||
}
|
||||
|
||||
// parse a timeseries file into a histogram
|
||||
// parse a time series file into a histogram
|
||||
fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
|
||||
let f = File::open(window_file)
|
||||
.chain_err(|| format!("Failed to open sample data file {}.", window_file))?;
|
||||
let mut buf = BufReader::new(&f);
|
||||
|
||||
// total number of bins is the product of all dimensions length
|
||||
let total_bins = cfg.num_bins.iter().fold(1, |s, &x| { s*x });
|
||||
let total_bins = cfg.num_bins.iter().product();
|
||||
let mut hist = vec![0.0; total_bins];
|
||||
|
||||
// bin width for each dimension: (max-min)/bins
|
||||
@@ -139,7 +136,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
|
||||
while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 {
|
||||
linecount += 1;
|
||||
// skip comments and empty lines
|
||||
if line.starts_with("#") || line.starts_with("@") || line.len() == 0 {
|
||||
if line.starts_with('#') || line.starts_with('@') || line.is_empty() {
|
||||
line.clear();
|
||||
continue;
|
||||
}
|
||||
@@ -158,7 +155,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
|
||||
}
|
||||
|
||||
if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) {
|
||||
let bin_indeces = (0..cfg.dimens).map(|dimen: usize| {
|
||||
let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
|
||||
let val = values[dimen+1];
|
||||
((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
|
||||
}).collect();
|
||||
@@ -174,20 +171,22 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
|
||||
}
|
||||
|
||||
// Write WHAM calculation results to out_file.
|
||||
pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, free_std: &Vec<f64>, prob: &Vec<f64>, prob_std: &Vec<f64>) -> Result<()> {
|
||||
pub fn write_results(out_file: &str, ds: &Dataset, free: &[f64],
|
||||
free_std: &[f64], prob: &[f64], prob_std: &[f64]) -> Result<()> {
|
||||
let output = File::create(out_file)
|
||||
.chain_err(|| format!("Failed to create file with path {}", out_file))?;
|
||||
let mut buf = BufWriter::new(output);
|
||||
|
||||
let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1))
|
||||
.collect::<Vec<String>>().join(" ");
|
||||
writeln!(buf, "#{} {} {} {} {}", header, "Free Energy", "+/-", "Probability", "+/-");
|
||||
writeln!(buf, "#{} Free Energy +/- Probability +/-", header).unwrap();
|
||||
|
||||
for bin in 0..free.len() {
|
||||
let coords = ds.get_coords_for_bin(bin);
|
||||
let coords_str: String = coords.iter().map(|c| {format!("{:8.6} ", c)})
|
||||
.collect::<Vec<String>>().join("\t");
|
||||
writeln!(buf, "{}{:8.6} {:8.6} {:8.6} {:8.6}", coords_str, free[bin], free_std[bin], prob[bin], prob_std[bin])
|
||||
writeln!(buf, "{}{:8.6} {:8.6} {:8.6} {:8.6}", coords_str,
|
||||
free[bin], free_std[bin], prob[bin], prob_std[bin])
|
||||
.chain_err(|| "Failed to write to file.")?;
|
||||
}
|
||||
Ok(())
|
||||
@@ -196,10 +195,12 @@ pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, free_std: &V
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use assert_approx_eq::assert_approx_eq;
|
||||
|
||||
fn cfg() -> Config {
|
||||
Config {
|
||||
metadata_file: "example/1d_cyclic/metadata.dat".to_string(),
|
||||
|
||||
hist_min: vec![-3.14],
|
||||
hist_max: vec![3.14],
|
||||
num_bins: vec![10],
|
||||
@@ -224,12 +225,12 @@ mod tests {
|
||||
let h = super::read_window_file(&f, &cfg).unwrap();
|
||||
println!("{:?}", h);
|
||||
assert_eq!(5000, h.num_points);
|
||||
assert_eq!(0.0, h.bins[2]);
|
||||
assert_eq!(11.0, h.bins[3]);
|
||||
assert_eq!(2236.0, h.bins[4]);
|
||||
assert_eq!(2714.0, h.bins[5]);
|
||||
assert_eq!(39.0, h.bins[6]);
|
||||
assert_eq!(0.0, h.bins[7]);
|
||||
assert_approx_eq!(0.0, h.bins[2]);
|
||||
assert_approx_eq!(11.0, h.bins[3]);
|
||||
assert_approx_eq!(2236.0, h.bins[4]);
|
||||
assert_approx_eq!(2714.0, h.bins[5]);
|
||||
assert_approx_eq!(39.0, h.bins[6]);
|
||||
assert_approx_eq!(0.0, h.bins[7]);
|
||||
}
|
||||
|
||||
|
||||
@@ -241,7 +242,7 @@ mod tests {
|
||||
assert_eq!(25, ds.num_windows);
|
||||
assert_eq!(cfg.num_bins.len(), ds.dimens_lengths.len());
|
||||
assert_eq!(cfg.num_bins[0], ds.dimens_lengths[0]);
|
||||
assert_eq!(cfg.temperature * k_B, ds.kT);
|
||||
assert_approx_eq!(cfg.temperature * k_B, ds.kT);
|
||||
assert_eq!(25, ds.histograms.len())
|
||||
}
|
||||
|
||||
|
||||
149
src/lib.rs
149
src/lib.rs
@@ -5,6 +5,10 @@ extern crate error_chain;
|
||||
extern crate rand;
|
||||
extern crate rgsl;
|
||||
extern crate rayon;
|
||||
#[cfg(test)]
|
||||
#[macro_use]
|
||||
extern crate assert_approx_eq;
|
||||
|
||||
|
||||
pub mod io;
|
||||
pub mod histogram;
|
||||
@@ -21,7 +25,7 @@ pub mod errors { error_chain!{} }
|
||||
use errors::*;
|
||||
|
||||
#[allow(non_upper_case_globals)]
|
||||
static k_B: f64 = 0.0083144621; // kJ/mol*K
|
||||
static k_B: f64 = 0.008_314_462_1; // kJ/mol*K
|
||||
|
||||
// Application config
|
||||
#[derive(Debug)]
|
||||
@@ -45,7 +49,9 @@ pub struct Config {
|
||||
|
||||
impl fmt::Display for Config {
|
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
|
||||
write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?} verbose={}, tolerance={}, iterations={}, temperature={}, cyclic={:?}, bootstrap={:?}, seed={:?}",
|
||||
write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?},
|
||||
verbose={}, tolerance={}, iterations={}, temperature={},
|
||||
cyclic={:?}, bootstrap={:?}, seed={:?}",
|
||||
self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
|
||||
self.verbose, self.tolerance, self.max_iterations, self.temperature,
|
||||
self.cyclic, self.bootstrap, self.bootstrap_seed)
|
||||
@@ -56,25 +62,34 @@ impl fmt::Display for Config {
|
||||
// converged if the maximal difference for the calculated bias offsets is
|
||||
// smaller then a tolerance value.
|
||||
fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool {
|
||||
!new_F.iter().zip(old_F.iter())
|
||||
.map(|x| { (x.0-x.1).abs() })
|
||||
.any(|diff| { diff > tolerance })
|
||||
// calculates abs diff between every old and new F and checks if any
|
||||
// is larger than tolerance
|
||||
!new_F.iter()
|
||||
.zip(old_F.iter())
|
||||
.map(|x| { (x.0-x.1).abs() })
|
||||
.any(|diff| { diff > tolerance })
|
||||
}
|
||||
|
||||
// estimate the probability of a bin of the histogram set based on given bias offsets (F)
|
||||
// This evaluates the first WHAM equation for each bin.
|
||||
// estimate the probability of a bin of the histogram set based on given bias
|
||||
// offsets (F). This evaluates the first WHAM equation for each bin:
|
||||
// P(x) = \frac {\sum_{i=1}^N{n_i(x)}}
|
||||
// {\sum_{i=1}^N{ N_i exp(\beta [F_i - U_{bias,i}(x)])}}
|
||||
fn calc_bin_probability(bin: usize, dataset: &Dataset, F: &[f64]) -> f64 {
|
||||
let mut denom_sum: f64 = 0.0;
|
||||
let bin_count: f64 = dataset.get_weighted_bin_count(bin);
|
||||
for (window, h) in dataset.histograms.iter().enumerate() {
|
||||
let bias = dataset.get_bias(bin, window);
|
||||
denom_sum += (dataset.weights[window] * h.num_points as f64) * bias * F[window];
|
||||
denom_sum += (dataset.weights[window] * h.num_points as f64)
|
||||
* bias * F[window];
|
||||
}
|
||||
bin_count / denom_sum
|
||||
}
|
||||
|
||||
// estimate the bias offset F of the histogram based on given probabilities.
|
||||
// This evaluates the second WHAM equation for each window and returns exp(F/kT)
|
||||
// This evaluates the second WHAM equation for each window and returns exp(F/kT).
|
||||
// exp(F/kT) is not required in intermediate steps so we save some time by not
|
||||
// calculating it for every iteration.
|
||||
// F_i = - 1/\beta ln[\sum_{X_{bins}}{P(x)exp(-\beta U_{bias,i}(x))}]
|
||||
fn calc_window_F(window: usize, dataset: &Dataset, P: &[f64]) -> f64 {
|
||||
let f: f64 = (0..dataset.num_bins).zip(P.iter()) // zip bins and P
|
||||
.map(|bin_and_prob: (usize, &f64)| {
|
||||
@@ -84,16 +99,18 @@ fn calc_window_F(window: usize, dataset: &Dataset, P: &[f64]) -> f64 {
|
||||
1.0/f
|
||||
}
|
||||
|
||||
// One full WHAM iteration includes calculation of new probabilities P and
|
||||
// new bias offsets F based on previous bias offsets F_prev. This updates
|
||||
// the values in vectors F and P
|
||||
// One full WHAM iteration: calculation of new probabilities P and new bias
|
||||
// offsets F based on previous bias offsets F_prev. This updates the values in
|
||||
// vectors F and P.
|
||||
fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut Vec<f64>, P: &mut Vec<f64>) {
|
||||
// evaluate first WHAM equation for each bin to
|
||||
// estimage probabilities based on previous offsets (F_prev))
|
||||
(0..dataset.num_bins).into_par_iter()
|
||||
// Update P
|
||||
// evaluate first WHAM equation for each bin to
|
||||
// estimate probabilities based on previous offsets (F_prev))
|
||||
(0..dataset.num_bins).into_par_iter()
|
||||
.map(|bin| { calc_bin_probability(bin, dataset, F_prev) })
|
||||
.collect_into_vec(P);
|
||||
|
||||
// Update F
|
||||
// evaluate second WHAM equation for each window to
|
||||
// estimate new bias offsets from propabilities
|
||||
(0..dataset.num_windows).into_par_iter()
|
||||
@@ -101,12 +118,20 @@ fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut Vec<f64>, P
|
||||
.collect_into_vec(F);
|
||||
}
|
||||
|
||||
pub fn perform_wham(cfg: &Config, dataset: &Dataset) -> Result<(Vec<f64>, Vec<f64>, Vec<f64>)> {
|
||||
// Full WHAM calculation. Calls `perform_wham_iteration` until convergence
|
||||
// criteria are met or max iterations reached.
|
||||
pub fn perform_wham(cfg: &Config, dataset: &Dataset)
|
||||
-> Result<(Vec<f64>, Vec<f64>, Vec<f64>)> {
|
||||
// allocate required vectors.
|
||||
let mut P: Vec<f64> = vec![f64::NAN; dataset.num_bins]; // bin probability
|
||||
let mut F: Vec<f64> = vec![1.0; dataset.num_windows]; // bias offset exp(F/kT)
|
||||
let mut F_prev: Vec<f64> = vec![f64::NAN; dataset.num_windows]; // previous bias offset
|
||||
let mut F_tmp: Vec<f64> = vec![f64::NAN; dataset.num_windows]; // temp storage for F
|
||||
|
||||
// bin probability
|
||||
let mut P: Vec<f64> = vec![f64::NAN; dataset.num_bins];
|
||||
// bias offset exp(F/kT)
|
||||
let mut F: Vec<f64> = vec![1.0; dataset.num_windows];
|
||||
// previous bias offset
|
||||
let mut F_prev: Vec<f64> = vec![f64::NAN; dataset.num_windows];
|
||||
// temp storage for F
|
||||
let mut F_tmp: Vec<f64> = vec![f64::NAN; dataset.num_windows];
|
||||
|
||||
let mut iteration = 0;
|
||||
let mut converged = false;
|
||||
@@ -118,33 +143,32 @@ pub fn perform_wham(cfg: &Config, dataset: &Dataset) -> Result<(Vec<f64>, Vec<f6
|
||||
// store F values before the next iteration
|
||||
F_prev.copy_from_slice(&F);
|
||||
|
||||
// perform wham iteration (this updates F and P)
|
||||
// perform wham iteration (this updates F and P).
|
||||
perform_wham_iteration(&dataset, &F_prev, &mut F, &mut P);
|
||||
|
||||
// convergence check
|
||||
if iteration % 10 == 0 {
|
||||
// This backups exp(F/kT) in a temporary vector and calculates true F and F_prev for
|
||||
// convergence. Finally, F is restored. F_prev does not need to be restored because
|
||||
// its overwritten for the next iteration.
|
||||
// This backups exp(F/kT) in a temporary vector and calculates
|
||||
// true F and F_prev for convergence. Finally, F is restored.
|
||||
// F_prev does not need to be restored because its overwritten
|
||||
// for the next iteration.
|
||||
F_tmp.copy_from_slice(&F);
|
||||
for f in F.iter_mut() { *f = -dataset.kT * f.ln() }
|
||||
for f in F_prev.iter_mut() { *f = -dataset.kT * f.ln() }
|
||||
converged = is_converged(&F_prev, &F, cfg.tolerance);
|
||||
|
||||
println!("Iteration {}: dF={}", &iteration, &diff_avg(&F_prev, &F));
|
||||
if cfg.verbose {
|
||||
println!("Iteration {}: dF={}", &iteration, &diff_avg(&F_prev, &F));
|
||||
}
|
||||
F.copy_from_slice(&F_tmp);
|
||||
}
|
||||
|
||||
// Dump free energy and bias offsets
|
||||
//if iteration % 100 == 0 {
|
||||
// free_energy(&histograms, &mut P, &mut A);
|
||||
// dump_state(&histograms, &F, &F_prev, &P, &A);
|
||||
//}
|
||||
}
|
||||
|
||||
// Normalize P to sum(P) = 1.0
|
||||
let P_sum: f64 = P.iter().sum();
|
||||
P.iter_mut().map(|p| *p /= P_sum).count();
|
||||
for p in P.iter_mut() {
|
||||
*p /= P_sum;
|
||||
}
|
||||
|
||||
if iteration == cfg.max_iterations {
|
||||
bail!("WHAM not converged! (max iterations reached)");
|
||||
@@ -161,17 +185,14 @@ pub fn run(cfg: &Config) -> Result<()>{
|
||||
println!("{}", &dataset);
|
||||
|
||||
let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
|
||||
println!("WHAM converged.");
|
||||
|
||||
let P_std: Vec<f64>;
|
||||
let free_energy_std: Vec<f64>;
|
||||
if cfg.bootstrap > 0 {
|
||||
let error_est = error_analysis::run_bootstrap(&cfg, dataset.clone(), &P, cfg.bootstrap);
|
||||
P_std = error_est.0;
|
||||
free_energy_std = error_est.1;
|
||||
let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
|
||||
println!("Bootstrapping..");
|
||||
error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
|
||||
} else {
|
||||
P_std = vec![0.0; P.len()];
|
||||
free_energy_std = vec![0.0; P.len()];
|
||||
}
|
||||
(vec![0.0; P.len()], vec![0.0; P.len()])
|
||||
};
|
||||
|
||||
// calculate free energy and dump state
|
||||
println!("Finished. Dumping final PMF");
|
||||
@@ -214,19 +235,23 @@ fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec<f64> {
|
||||
free_energy
|
||||
}
|
||||
|
||||
fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], P_std: &[f64], A: &[f64], A_std: &[f64]) {
|
||||
// Print the current WHAM iteration state. Dumps the PMF and associated vectors
|
||||
fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64],
|
||||
P_std: &[f64], A: &[f64], A_std: &[f64]) {
|
||||
// TODO fix output of F/F_prev
|
||||
let out = std::io::stdout();
|
||||
let mut lock = out.lock();
|
||||
writeln!(lock, "# PMF");
|
||||
writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-");
|
||||
writeln!(lock, "# PMF").unwrap();
|
||||
writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-").unwrap();
|
||||
for bin in 0..dataset.num_bins {
|
||||
writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}", bin, A[bin], A_std[bin], P[bin], P_std[bin]);
|
||||
writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}",
|
||||
bin, A[bin], A_std[bin], P[bin], P_std[bin]).unwrap();
|
||||
}
|
||||
writeln!(lock, "# Bias offsets");
|
||||
writeln!(lock, "#Window\t\tF\t\tF_prev");
|
||||
writeln!(lock, "# Bias offsets").unwrap();
|
||||
writeln!(lock, "#Window\t\tF\t\tF_prev").unwrap();
|
||||
for window in 0..dataset.num_windows {
|
||||
writeln!(lock, "{}\t{:9.5}\t{:8.8}", window, F[window], (F[window]-F_prev[window]).abs());
|
||||
writeln!(lock, "{}\t{:9.5}\t{:8.8}",
|
||||
window, F[window], (F[window]-F_prev[window]).abs()).unwrap();
|
||||
}
|
||||
}
|
||||
|
||||
@@ -268,23 +293,23 @@ mod tests {
|
||||
fn calc_bin_probability() {
|
||||
let dataset = create_test_dataset();
|
||||
let F = vec![1.0; dataset.num_bins] ;
|
||||
let expected = vec!(0.0, 0.0825296687031316, 40.92355847097493,
|
||||
124226.70003377, 2308526035.5283747);
|
||||
for b in 0..dataset.num_bins {
|
||||
let p = super::calc_bin_probability(b, &dataset, &F);
|
||||
assert_delta!(expected[b], p, 0.0000001);
|
||||
}
|
||||
}
|
||||
let expected = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93,
|
||||
124_226.700_033_77, 2_308_526_035.528_374_7);
|
||||
expected.iter().enumerate().for_each(|(i, exp)| {
|
||||
let p = super::calc_bin_probability(i, &dataset, &F);
|
||||
assert_delta!(exp, p, 0.000_000_1);
|
||||
})
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn calc_bias_offset() {
|
||||
let dataset = create_test_dataset();
|
||||
let probability = vec!(0.0, 0.1, 0.2, 0.3, 0.4);
|
||||
let expected = vec!(15.927477169990633, 15.927477169990633);
|
||||
for window in 0..dataset.num_windows {
|
||||
let F = super::calc_window_F(window, &dataset, &probability);
|
||||
assert_delta!(expected[window], F, 0.0000001);
|
||||
}
|
||||
let expected = vec!(15.927_477_169_990_633, 15.927_477_169_990_633);
|
||||
expected.iter().enumerate().for_each(|(i, exp)| {
|
||||
let F = super::calc_window_F(i, &dataset, &probability);
|
||||
assert_delta!(exp, F, 0.000_000_1);
|
||||
})
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -295,8 +320,8 @@ mod tests {
|
||||
let mut P = vec![f64::NAN; dataset.num_bins];
|
||||
super::perform_wham_iteration(&dataset, &prev_F, &mut F, &mut P);
|
||||
let expected_F = vec!(1.0, 1.0);
|
||||
let expected_P = vec!(0.0, 0.0825296687031316, 40.92355847097493,
|
||||
124226.70003377, 2308526035.5283747);
|
||||
let expected_P = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93,
|
||||
124_226.700_033_77, 2_308_526_035.528_374_7);
|
||||
for bin in 0..dataset.num_bins {
|
||||
assert_delta!(expected_P[bin], P[bin], 0.01)
|
||||
}
|
||||
|
||||
@@ -3,6 +3,7 @@ extern crate wham;
|
||||
extern crate clap;
|
||||
extern crate rand;
|
||||
|
||||
use rand::prelude::*;
|
||||
use clap::App;
|
||||
use wham::Config;
|
||||
use wham::errors::*;
|
||||
@@ -49,7 +50,6 @@ fn cli() -> Result<Config> {
|
||||
.chain_err(|| "Cannot parse bootstrap iteration.")?;
|
||||
let bootstrap_seed: u64 = matches.value_of("bootstrap_seed")
|
||||
.unwrap_or({
|
||||
use rand::Rng;
|
||||
let mut rng = rand::thread_rng();
|
||||
&rng.gen::<u32>().to_string()
|
||||
}).parse()
|
||||
|
||||
@@ -2,8 +2,6 @@ mod command;
|
||||
|
||||
#[cfg(test)]
|
||||
mod integration {
|
||||
|
||||
use std::process::Command;
|
||||
use super::command::get_command;
|
||||
|
||||
#[test]
|
||||
|
||||
Reference in New Issue
Block a user