diff --git a/Cargo.lock b/Cargo.lock index 9d6bab0..a7bad1b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1,12 +1,12 @@ # This file is automatically @generated by Cargo. # It is not intended for manual editing. -version = 3 +version = 4 [[package]] -name = "adler" -version = "1.0.2" +name = "adler2" +version = "2.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f26201604c87b1e01bd3d98f8d5d9a8fcbb815e8cedb41ffccbeb4bf593a35fe" +checksum = "512761e0bb2578dd7380c6baaa0f4ce03e84f95e960231d1dec8bf4d7d6e2627" [[package]] name = "aho-corasick" @@ -107,33 +107,33 @@ checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" [[package]] name = "base16ct" -version = "0.1.1" +version = "0.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "349a06037c7bf932dd7e7d1f653678b2038b9ad46a74102f1fc7bd7872678cce" +checksum = "4c7f02d4ea65f2c1853089ffd8d2787bdbc63de2f0d29dedbcf8ccdfa0ccd4cf" [[package]] name = "base64" -version = "0.13.1" +version = "0.22.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9e1b586273c5702936fe7b7d6896644d8be71e6314cfe09d3167c95f712589e8" +checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" [[package]] -name = "base64" -version = "0.21.2" +name = "bitflags" +version = "1.3.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "604178f6c5c21f02dc555784810edfb88d34ac2c73b2eae109655649ee73ce3d" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" [[package]] name = "bitflags" -version = "1.3.2" +version = "2.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" +checksum = "5c8214115b7bf84099f1309324e63141d4c5d7cc26862f97a0a857dbefe165bd" [[package]] name = "blas" -version = "0.22.0" +version = "0.23.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ae980f75c3215bfe8203c349b28149b0f4130a262e072913ccb55f877cd239dc" +checksum = "280f3f48ac4edf086d59c7aa71fb28b78df2b14db8f591fb4eb29b88b77ac6ad" dependencies = [ "blas-sys", "libc", @@ -142,9 +142,9 @@ dependencies = [ [[package]] name = "blas-sys" -version = "0.7.1" +version = "0.8.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "13b1b279ceb25d7c4faaea95a5f7addbe7d8c34f9462044bd8e630cebcfc2440" +checksum = "b49cc48acc7ead5c33d2e3cac9c47e55bfd04c929abd1456459084a8c8cdf2cb" dependencies = [ "libc", ] @@ -172,9 +172,12 @@ checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" [[package]] name = "cc" -version = "1.0.79" +version = "1.2.16" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "50d30906286121d95be3d479533b458f87493b30a4b5f79a607db8f5d11aa91f" +checksum = "be714c154be609ec7f5dad223a33bf1482fff90472de28f7362806e6d4832b8c" +dependencies = [ + "shlex", +] [[package]] name = "cfg-if" @@ -184,17 +187,16 @@ checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" [[package]] name = "chrono" -version = "0.4.26" +version = "0.4.40" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ec837a71355b28f6556dbd569b37b3f363091c0bd4b2e735674521b4c5fd9bc5" +checksum = "1a7964611d71df112cb1730f2ee67324fcf4d0fc6606acbbe9bfe06df124637c" dependencies = [ "android-tzdata", "iana-time-zone", "js-sys", "num-traits", - "time", "wasm-bindgen", - "winapi", + "windows-link", ] [[package]] @@ -240,7 +242,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "acd4f3c17c83b0ba34ffbc4f8bbd74f079413f747f84a6f89292f138057e36ab" dependencies = [ "anstyle", - "bitflags", + "bitflags 1.3.2", "clap_lex", ] @@ -256,6 +258,15 @@ version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0b6a852b24ab71dffc585bcb46eaf7959d175cb865a7152e35b348d1b2960422" +[[package]] +name = "colored" +version = "3.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fde0e0ec90c9dfb3b4b1a0891a7dcd0e2bffde2f7efed5fe7c9bb00e5bfb915e" +dependencies = [ + "windows-sys 0.52.0", +] + [[package]] name = "console" version = "0.15.8" @@ -277,15 +288,15 @@ checksum = "6245d59a3e82a7fc217c5828a6692dbc6dfb63a0c8c90495621f7b9d79704a0e" [[package]] name = "core-foundation-sys" -version = "0.8.4" +version = "0.8.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e496a50fda8aacccc86d7529e2c1e0892dbd0f898a6b5645b5561b89c3210efa" +checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" [[package]] name = "cpufeatures" -version = "0.2.9" +version = "0.2.17" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a17b76ff3a4162b0b27f354a0c87015ddad39d35f9c0c36607a3bdd175dde1f1" +checksum = "59ed5838eebb26a2bb2e58f6d5b5316989ae9d08bab10e0e6d103e656d1b0280" dependencies = [ "libc", ] @@ -401,9 +412,9 @@ dependencies = [ [[package]] name = "darling" -version = "0.14.4" +version = "0.20.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7b750cb3417fd1b327431a470f388520309479ab0bf5e323505daf0290cd3850" +checksum = "6f63b86c8a8826a49b8c21f08a2d07338eec8d900540f8630dc76284be802989" dependencies = [ "darling_core", "darling_macro", @@ -411,58 +422,58 @@ dependencies = [ [[package]] name = "darling_core" -version = "0.14.4" +version = "0.20.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "109c1ca6e6b7f82cc233a97004ea8ed7ca123a9af07a8230878fcfda9b158bf0" +checksum = "95133861a8032aaea082871032f5815eb9e98cef03fa916ab4500513994df9e5" dependencies = [ "fnv", "ident_case", "proc-macro2", "quote", "strsim", - "syn 1.0.109", + "syn 2.0.99", ] [[package]] name = "darling_macro" -version = "0.14.4" +version = "0.20.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a4aab4dbc9f7611d8b55048a3a16d2d010c2c8334e46304b40ac1cc14bf3b48e" +checksum = "d336a2a514f6ccccaa3e09b02d41d35330c07ddf03a62165fcec10bb561c7806" dependencies = [ "darling_core", "quote", - "syn 1.0.109", + "syn 2.0.99", ] [[package]] name = "derive_builder" -version = "0.11.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d07adf7be193b71cc36b193d0f5fe60b918a3a9db4dad0449f57bcfd519704a3" +checksum = "507dfb09ea8b7fa618fcf76e953f4f5e192547945816d5358edffe39f6f94947" dependencies = [ "derive_builder_macro", ] [[package]] name = "derive_builder_core" -version = "0.11.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1f91d4cfa921f1c05904dc3c57b4a32c38aed3340cce209f3a6fd1478babafc4" +checksum = "2d5bcf7b024d6835cfb3d473887cd966994907effbe9227e8c8219824d06c4e8" dependencies = [ "darling", "proc-macro2", "quote", - "syn 1.0.109", + "syn 2.0.99", ] [[package]] name = "derive_builder_macro" -version = "0.11.2" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8f0314b72bed045f3a68671b3c86328386762c93f82d98c65c3cb5e5f573dd68" +checksum = "ab63b0e2bf4d5928aff72e83a7dace85d7bba5fe12dcc3c5a572d78caffd3f3c" dependencies = [ "derive_builder_core", - "syn 1.0.109", + "syn 2.0.99", ] [[package]] @@ -490,22 +501,34 @@ dependencies = [ [[package]] name = "directories" -version = "4.0.1" +version = "5.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f51c5d4ddabd36886dd3e1438cb358cdcb0d7c499cb99cb4ac2e38e18b5cb210" +checksum = "9a49173b84e034382284f27f1af4dcbbd231ffa358c0fe316541a7337f376a35" dependencies = [ "dirs-sys", ] [[package]] name = "dirs-sys" -version = "0.3.7" +version = "0.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1b1d1d91c932ef41c0f2663aa8b0ca0342d444d842c06914aa0a7e352d0bada6" +checksum = "520f05a5cbd335fae5a99ff7a6ab8627577660ee5cfd6a94a6a929b52ff0321c" dependencies = [ "libc", + "option-ext", "redox_users", - "winapi", + "windows-sys 0.48.0", +] + +[[package]] +name = "displaydoc" +version = "0.2.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "97369cbbc041bc366949bc74d34658d6cda5621039731c6310521892a3a20ae0" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.99", ] [[package]] @@ -544,24 +567,19 @@ dependencies = [ ] [[package]] -name = "errno" -version = "0.3.1" +name = "equivalent" +version = "1.0.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4bcfec3a70f97c962c307b2d2c56e358cf1d00b558d74262b5f929ee8cc7e73a" -dependencies = [ - "errno-dragonfly", - "libc", - "windows-sys 0.48.0", -] +checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f" [[package]] -name = "errno-dragonfly" -version = "0.1.2" +name = "errno" +version = "0.3.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "aa68f1b12764fab894d2755d2518754e71b4fd80ecfb822714a1206c2aab39bf" +checksum = "33d852cb9b869c2a9b3df2f71a3074817f01e1844f839a144f5fcef059a4eb5d" dependencies = [ - "cc", "libc", + "windows-sys 0.52.0", ] [[package]] @@ -578,9 +596,9 @@ dependencies = [ [[package]] name = "flate2" -version = "1.0.26" +version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3b9429470923de8e8cbd4d2dc513535400b4b3fef0319fb5c4e1f520a7bef743" +checksum = "11faaf5a5236997af9848be0bef4db95824b1d534ebc64d0f0c6cf3e67bd38dc" dependencies = [ "crc32fast", "miniz_oxide", @@ -594,9 +612,9 @@ checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" [[package]] name = "form_urlencoded" -version = "1.2.0" +version = "1.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a62bc1cf6f830c2ec14a513a9fb124d0a213a629668a4186f329db21fe045652" +checksum = "e13624c2627564efccf4934284bdd98cbaa14e79b0b5a141218e507b3a823456" dependencies = [ "percent-encoding", ] @@ -622,16 +640,28 @@ dependencies = [ "wasi 0.11.0+wasi-snapshot-preview1", ] +[[package]] +name = "getrandom" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "43a49c392881ce6d5c3b8cb70f98717b7c07aabbdff06687b9030dbfbe2725f8" +dependencies = [ + "cfg-if", + "libc", + "wasi 0.13.3+wasi-0.2.2", + "windows-targets 0.52.5", +] + [[package]] name = "getset" -version = "0.1.2" +version = "0.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e45727250e75cc04ff2846a66397da8ef2b3db8e40e0cef4df67950a07621eb9" +checksum = "f3586f256131df87204eb733da72e3d3eb4f343c639f4b7be279ac7c48baeafe" dependencies = [ - "proc-macro-error", + "proc-macro-error2", "proc-macro2", "quote", - "syn 1.0.109", + "syn 2.0.99", ] [[package]] @@ -640,12 +670,24 @@ version = "1.8.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "eabb4a44450da02c90444cf74558da904edde8fb4e9035a9a6a4e15445af0bd7" +[[package]] +name = "hashbrown" +version = "0.15.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bf151400ff0baff5465007dd2f3e717f3fe502074ca563069ce3a6629d07b289" + [[package]] name = "heck" version = "0.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + [[package]] name = "hermit-abi" version = "0.3.1" @@ -660,16 +702,16 @@ checksum = "9a3a5bfb195931eeb336b2a7b4d761daec841b97f947d34394601737a7bba5e4" [[package]] name = "iana-time-zone" -version = "0.1.57" +version = "0.1.61" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2fad5b825842d2b38bd206f3e81d6957625fd7f0a361e345c30e01a0ae2dd613" +checksum = "235e081f3925a06703c2d0117ea8b91f042756fd6e7a6e5d901e8ca1a996b220" dependencies = [ "android_system_properties", "core-foundation-sys", "iana-time-zone-haiku", "js-sys", "wasm-bindgen", - "windows", + "windows-core", ] [[package]] @@ -681,6 +723,124 @@ dependencies = [ "cc", ] +[[package]] +name = "icu_collections" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "db2fa452206ebee18c4b5c2274dbf1de17008e874b4dc4f0aea9d01ca79e4526" +dependencies = [ + "displaydoc", + "yoke", + "zerofrom", + "zerovec", +] + +[[package]] +name = "icu_locid" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "13acbb8371917fc971be86fc8057c41a64b521c184808a698c02acc242dbf637" +dependencies = [ + "displaydoc", + "litemap", + "tinystr", + "writeable", + "zerovec", +] + +[[package]] +name = "icu_locid_transform" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01d11ac35de8e40fdeda00d9e1e9d92525f3f9d887cdd7aa81d727596788b54e" +dependencies = [ + "displaydoc", + "icu_locid", + "icu_locid_transform_data", + "icu_provider", + "tinystr", + "zerovec", +] + +[[package]] +name = "icu_locid_transform_data" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fdc8ff3388f852bede6b579ad4e978ab004f139284d7b28715f773507b946f6e" + +[[package]] +name = "icu_normalizer" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19ce3e0da2ec68599d193c93d088142efd7f9c5d6fc9b803774855747dc6a84f" +dependencies = [ + "displaydoc", + "icu_collections", + "icu_normalizer_data", + "icu_properties", + "icu_provider", + "smallvec", + "utf16_iter", + "utf8_iter", + "write16", + "zerovec", +] + +[[package]] +name = "icu_normalizer_data" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8cafbf7aa791e9b22bec55a167906f9e1215fd475cd22adfcf660e03e989516" + +[[package]] +name = "icu_properties" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93d6020766cfc6302c15dbbc9c8778c37e62c14427cb7f6e601d849e092aeef5" +dependencies = [ + "displaydoc", + "icu_collections", + "icu_locid_transform", + "icu_properties_data", + "icu_provider", + "tinystr", + "zerovec", +] + +[[package]] +name = "icu_properties_data" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "67a8effbc3dd3e4ba1afa8ad918d5684b8868b3b26500753effea8d2eed19569" + +[[package]] +name = "icu_provider" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ed421c8a8ef78d3e2dbc98a973be2f3770cb42b606e3ab18d6237c4dfde68d9" +dependencies = [ + "displaydoc", + "icu_locid", + "icu_provider_macros", + "stable_deref_trait", + "tinystr", + "writeable", + "yoke", + "zerofrom", + "zerovec", +] + +[[package]] +name = "icu_provider_macros" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ec89e9337638ecdc08744df490b221a7399bf8d164eb52a665454e60e075ad6" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.99", +] + [[package]] name = "ident_case" version = "1.0.1" @@ -689,12 +849,23 @@ checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" [[package]] name = "idna" -version = "0.4.0" +version = "1.0.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7d20d6b07bfbc108882d88ed8e37d39636dcc260e15e30c45e6ba089610b917c" +checksum = "686f825264d630750a544639377bae737628043f20d38bbc029e8f29ea968a7e" dependencies = [ - "unicode-bidi", - "unicode-normalization", + "idna_adapter", + "smallvec", + "utf8_iter", +] + +[[package]] +name = "idna_adapter" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "daca1df1c957320b2cf139ac61e7bd64fed304c5040df000a745aa1de3b4ef71" +dependencies = [ + "icu_normalizer", + "icu_properties", ] [[package]] @@ -704,6 +875,7 @@ dependencies = [ "assert", "blas", "blas-sys", + "colored", "criterion", "csv", "derive_more", @@ -722,6 +894,16 @@ dependencies = [ "rayon", ] +[[package]] +name = "indexmap" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8c9c992b02b5b4c94ea26e32fe5bccb7aa7d9f390ab5c1221ff895bc7ea8b652" +dependencies = [ + "equivalent", + "hashbrown", +] + [[package]] name = "indicatif" version = "0.17.8" @@ -791,7 +973,7 @@ checksum = "adcf93614601c8129ddf72e2d5633df827ba6551541c6d8c59520a371475be1f" dependencies = [ "hermit-abi", "io-lifetimes", - "rustix", + "rustix 0.37.20", "windows-sys 0.48.0", ] @@ -818,18 +1000,19 @@ checksum = "453ad9f582a441959e5f0d088b02ce04cfe8d51a8eaf077f12ac6d3e94164ca6" [[package]] name = "js-sys" -version = "0.3.64" +version = "0.3.77" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c5f195fe497f702db0f318b07fdd68edb16955aed830df8363d837542f8f935a" +checksum = "1cfaf33c695fc6e08064efbc1f72ec937429614f25eef83af942d0e227c3a28f" dependencies = [ + "once_cell", "wasm-bindgen", ] [[package]] name = "lapack" -version = "0.19.0" +version = "0.20.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ad676a6b4df7e76a9fd80a0c50c619a3948d6105b62a0ab135f064d99c51d207" +checksum = "508f7f48955e5e4ec3feb771a4d981c6e59ff2400d750cf41733a3350fdb2dc4" dependencies = [ "lapack-sys", "libc", @@ -838,9 +1021,9 @@ dependencies = [ [[package]] name = "lapack-sys" -version = "0.14.0" +version = "0.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "447f56c85fb410a7a3d36701b2153c1018b1d2b908c5fbaf01c1b04fac33bcbe" +checksum = "314b879030845b68571809a6978e52d3b67ac5fba07e77b1b317b484092e2fb5" dependencies = [ "libc", ] @@ -853,9 +1036,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" [[package]] name = "libc" -version = "0.2.146" +version = "0.2.170" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f92be4933c13fd498862a9e02a3055f8a8d9c039ce33db97306fd5a6caa7f29b" +checksum = "875b3680cb2f8f71bdcf9a30f38d48282f5d3c95cbf9b3fa57269bb5d5c06828" [[package]] name = "linux-raw-sys" @@ -863,6 +1046,18 @@ version = "0.3.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ef53942eb7bf7ff43a617b3e2c1c4a5ecf5944a7c1bc12d7ee39bbb15e5c1519" +[[package]] +name = "linux-raw-sys" +version = "0.9.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6db9c683daf087dc577b7506e9695b3d556a9f3849903fa28186283afd6809e9" + +[[package]] +name = "litemap" +version = "0.7.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23fb14cb19457329c82206317a5663005a4d404783dc74f4252769b0d5f42856" + [[package]] name = "lock_api" version = "0.4.11" @@ -881,9 +1076,9 @@ checksum = "90ed8c1e510134f979dbc4f070f87d4313098b704861a105fe34231c70a3901c" [[package]] name = "memchr" -version = "2.5.0" +version = "2.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2dffe52ecf27772e601905b7522cb4ef790d2cc203488bbd0e2fe85fcb74566d" +checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" [[package]] name = "memoffset" @@ -896,11 +1091,11 @@ dependencies = [ [[package]] name = "miniz_oxide" -version = "0.7.1" +version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e7810e0be55b428ada41041c41f32c9f1a42817901b4ccf45fa3d4b6561e74c7" +checksum = "8e3e04debbb59698c15bacbb6d93584a8c0ca9cc3213cb423d31f760d8843ce5" dependencies = [ - "adler", + "adler2", ] [[package]] @@ -987,25 +1182,29 @@ checksum = "830b246a0e5f20af87141b25c173cd1b609bd7779a4617d6ec582abaf90870f3" [[package]] name = "oci-spec" -version = "0.5.8" +version = "0.6.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "98135224dd4faeb24c05a2fac911ed53ea6b09ecb09d7cada1cb79963ab2ee34" +checksum = "3f5a3fe998d50101ae009351fec56d88a69f4ed182e11000e711068c2f5abf72" dependencies = [ "derive_builder", "getset", + "once_cell", + "regex", "serde", "serde_json", + "strum", + "strum_macros", "thiserror", ] [[package]] name = "ocipkg" -version = "0.2.8" +version = "0.2.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60cf01280832705a4e4c4d046d9e67a54b900297c69191457a8fc6d198ddefa2" +checksum = "9bb3293021f06540803301af45e7ab81693d50e89a7398a3420bdab139e7ba5e" dependencies = [ "base16ct", - "base64 0.13.1", + "base64", "chrono", "directories", "flate2", @@ -1027,9 +1226,9 @@ dependencies = [ [[package]] name = "once_cell" -version = "1.18.0" +version = "1.19.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" +checksum = "3fdb12b2476b595f9358c5161aa467c2438859caa136dec86c26fdd2efe17b92" [[package]] name = "oorandom" @@ -1037,6 +1236,12 @@ version = "11.1.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0ab1bc2a289d34bd04a330323ac98a1b4bc82c9d9fcb1e66b63caa84da26b575" +[[package]] +name = "option-ext" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "04744f49eae99ab78e0d5c0b603ab218f515ea8cfe5a456d7629ad883a3b6e7d" + [[package]] name = "parking_lot" version = "0.12.1" @@ -1062,9 +1267,9 @@ dependencies = [ [[package]] name = "percent-encoding" -version = "2.3.0" +version = "2.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9b2a4787296e9989611394c33f193f676704af1686e70b8f8033ab5ba9a35a94" +checksum = "e3148f5046208a5d56bcfc03053e3ca6334e51da8dfb19b6cdc8b306fae3283e" [[package]] name = "pfapack" @@ -1128,34 +1333,32 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de" [[package]] -name = "proc-macro-error" -version = "1.0.4" +name = "proc-macro-error-attr2" +version = "2.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "da25490ff9892aab3fcf7c36f08cfb902dd3e71ca0f9f9517bea02a73a5ce38c" +checksum = "96de42df36bb9bba5542fe9f1a054b8cc87e172759a1868aa05c1f3acc89dfc5" dependencies = [ - "proc-macro-error-attr", "proc-macro2", "quote", - "syn 1.0.109", - "version_check", ] [[package]] -name = "proc-macro-error-attr" -version = "1.0.4" +name = "proc-macro-error2" +version = "2.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a1be40180e52ecc98ad80b184934baf3d0d29f979574e439af5a55274b35f869" +checksum = "11ec05c52be0a07b08061f7dd003e7d7092e0472bc731b4af7bb1ef876109802" dependencies = [ + "proc-macro-error-attr2", "proc-macro2", "quote", - "version_check", + "syn 2.0.99", ] [[package]] name = "proc-macro2" -version = "1.0.66" +version = "1.0.94" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "18fb31db3f9bddb2ea821cde30a9f70117e3f119938b5ee630b7403aa6e2ead9" +checksum = "a31971752e70b8b2686d7e46ec17fb38dad4051d94024c88df49b667caea9c84" dependencies = [ "unicode-ident", ] @@ -1207,7 +1410,7 @@ dependencies = [ "proc-macro2", "pyo3-macros-backend", "quote", - "syn 2.0.28", + "syn 2.0.99", ] [[package]] @@ -1216,18 +1419,18 @@ version = "0.21.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "08260721f32db5e1a5beae69a55553f56b99bd0e1c3e6e0a5e8851a9d0f5a85c" dependencies = [ - "heck", + "heck 0.4.1", "proc-macro2", "pyo3-build-config", "quote", - "syn 2.0.28", + "syn 2.0.99", ] [[package]] name = "quote" -version = "1.0.32" +version = "1.0.39" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "50f3b39ccfb720540debaa0164757101c08ecb8d326b15358ce76a62c7e85965" +checksum = "c1f1914ce909e1658d9907913b4b91947430c7d9be598b15a1912935b8c04801" dependencies = [ "proc-macro2", ] @@ -1259,7 +1462,7 @@ version = "0.6.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" dependencies = [ - "getrandom", + "getrandom 0.2.10", ] [[package]] @@ -1297,7 +1500,7 @@ version = "0.2.16" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "fb5a58c1855b4b6819d59012155603f0b22ad30cad752600aadfcb695265519a" dependencies = [ - "bitflags", + "bitflags 1.3.2", ] [[package]] @@ -1306,7 +1509,7 @@ version = "0.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4722d768eff46b75989dd134e5c353f0d6296e5aaa3132e776cbdb56be7731aa" dependencies = [ - "bitflags", + "bitflags 1.3.2", ] [[package]] @@ -1315,16 +1518,28 @@ version = "0.4.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b033d837a7cf162d7993aded9304e30a83213c648b6e389db233191f891e5c2b" dependencies = [ - "getrandom", + "getrandom 0.2.10", "redox_syscall 0.2.16", "thiserror", ] [[package]] name = "regex" -version = "1.8.4" +version = "1.10.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d0ab3ca65655bb1e41f2a8c8cd662eb4fb035e67c3f78da1d61dffe89d07300f" +checksum = "4219d74c6b67a3654a9fbebc4b419e22126d13d2f3c4a07ee0cb61ff79a79619" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "809e8dc61f6de73b46c85f4c96486310fe304c434cfa43669d7b40f711150908" dependencies = [ "aho-corasick", "memchr", @@ -1333,23 +1548,22 @@ dependencies = [ [[package]] name = "regex-syntax" -version = "0.7.2" +version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "436b050e76ed2903236f032a59761c1eb99e1b0aead2c257922771dab1fc8c78" +checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" [[package]] name = "ring" -version = "0.16.20" +version = "0.17.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3053cf52e236a3ed746dfc745aa9cacf1b791d846bdaf412f60a8d7d6e17c8fc" +checksum = "70ac5d832aa16abd7d1def883a8545280c20a60f523a370aa3a9617c2b8550ee" dependencies = [ "cc", + "cfg-if", + "getrandom 0.2.10", "libc", - "once_cell", - "spin", "untrusted", - "web-sys", - "winapi", + "windows-sys 0.52.0", ] [[package]] @@ -1367,46 +1581,65 @@ version = "0.37.20" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b96e891d04aa506a6d1f318d2771bcb1c7dfda84e126660ace067c9b474bb2c0" dependencies = [ - "bitflags", + "bitflags 1.3.2", "errno", "io-lifetimes", "libc", - "linux-raw-sys", + "linux-raw-sys 0.3.8", "windows-sys 0.48.0", ] [[package]] -name = "rustls" -version = "0.21.6" +name = "rustix" +version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1d1feddffcfcc0b33f5c6ce9a29e341e4cd59c3f78e7ee45f4a40c038b1d6cbb" +checksum = "dade4812df5c384711475be5fcd8c162555352945401aed22a35bffeab61f657" dependencies = [ - "log", - "ring", - "rustls-webpki 0.101.2", - "sct", + "bitflags 2.9.0", + "errno", + "libc", + "linux-raw-sys 0.9.2", + "windows-sys 0.52.0", ] [[package]] -name = "rustls-webpki" -version = "0.100.1" +name = "rustls" +version = "0.23.23" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d6207cd5ed3d8dca7816f8f3725513a34609c0c765bf652b8c3cb4cfd87db46b" +checksum = "47796c98c480fce5406ef69d1c76378375492c3b0a0de587be0c1d9feb12f395" dependencies = [ + "log", + "once_cell", "ring", - "untrusted", + "rustls-pki-types", + "rustls-webpki", + "subtle", + "zeroize", ] +[[package]] +name = "rustls-pki-types" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "917ce264624a4b4db1c364dcc35bfca9ded014d0a958cd47ad3e960e988ea51c" + [[package]] name = "rustls-webpki" -version = "0.101.2" +version = "0.102.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "513722fd73ad80a71f72b61009ea1b584bcfa1483ca93949c8f290298837fa59" +checksum = "64ca1bc8749bd4cf37b5ce386cc146580777b4e8572c7b97baf22c83f444bee9" dependencies = [ "ring", + "rustls-pki-types", "untrusted", ] +[[package]] +name = "rustversion" +version = "1.0.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eded382c5f5f786b989652c49544c4877d9f015cc22e145a5ea8ea66c2921cd2" + [[package]] name = "ryu" version = "1.0.13" @@ -1428,16 +1661,6 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" -[[package]] -name = "sct" -version = "0.7.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d53dcdb7c9f8158937a7981b48accfd39a43af418591a5d008c7b22b5e1b7ca4" -dependencies = [ - "ring", - "untrusted", -] - [[package]] name = "semver" version = "1.0.17" @@ -1446,46 +1669,62 @@ checksum = "bebd363326d05ec3e2f532ab7660680f3b02130d780c299bca73469d521bc0ed" [[package]] name = "serde" -version = "1.0.180" +version = "1.0.218" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0ea67f183f058fe88a4e3ec6e2788e003840893b91bac4559cabedd00863b3ed" +checksum = "e8dfc9d19bdbf6d17e22319da49161d5d0108e4188e8b680aef6299eed22df60" dependencies = [ "serde_derive", ] [[package]] name = "serde_derive" -version = "1.0.180" +version = "1.0.218" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "24e744d7782b686ab3b73267ef05697159cc0e5abbed3f47f9933165e5219036" +checksum = "f09503e191f4e797cb8aac08e9a4a4695c5edf6a2e70e376d961ddd5c969f82b" dependencies = [ "proc-macro2", "quote", - "syn 2.0.28", + "syn 2.0.99", ] [[package]] name = "serde_json" -version = "1.0.104" +version = "1.0.140" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "076066c5f1078eac5b722a31827a8832fe108bed65dfa75e233c89f8206e976c" +checksum = "20068b6e96dc6c9bd23e01df8827e6c7e1f2fddd43c21810382803c136b99373" dependencies = [ "itoa", + "memchr", "ryu", "serde", ] +[[package]] +name = "serde_spanned" +version = "0.6.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "87607cb1398ed59d48732e575a4c28a7a8ebf2454b964fe3f224f2afc07909e1" +dependencies = [ + "serde", +] + [[package]] name = "sha2" -version = "0.10.7" +version = "0.10.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "479fb9d862239e610720565ca91403019f2f00410f1864c5aa7479b950a76ed8" +checksum = "793db75ad2bcafc3ffa7c68b215fee268f537982cd901d132f89c6343f3a3dc8" dependencies = [ "cfg-if", "cpufeatures", "digest", ] +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + [[package]] name = "smallvec" version = "1.13.2" @@ -1493,16 +1732,41 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" [[package]] -name = "spin" -version = "0.5.2" +name = "stable_deref_trait" +version = "1.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6e63cff320ae2c57904679ba7cb63280a3dc4613885beafb148ee7bf9aa9042d" +checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3" [[package]] name = "strsim" -version = "0.10.0" +version = "0.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f" + +[[package]] +name = "strum" +version = "0.26.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "73473c0e59e6d5812c5dfe2a064a6444949f089e20eec9a2e5506596494e4623" +checksum = "8fec0f0aef304996cf250b31b5a10dee7980c85da9d759361292b8bca5a18f06" + +[[package]] +name = "strum_macros" +version = "0.26.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4c6bee85a5a24955dc440386795aa378cd9cf82acd5f764469152d2270e581be" +dependencies = [ + "heck 0.5.0", + "proc-macro2", + "quote", + "rustversion", + "syn 2.0.99", +] + +[[package]] +name = "subtle" +version = "2.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "13c2bddecc57b384dee18652358fb23172facb8a2c51ccc10d74c157bdea3292" [[package]] name = "syn" @@ -1517,20 +1781,31 @@ dependencies = [ [[package]] name = "syn" -version = "2.0.28" +version = "2.0.99" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "04361975b3f5e348b2189d8dc55bc942f278b2d482a6a0365de5bdd62d351567" +checksum = "e02e925281e18ffd9d640e234264753c43edc62d64b2d4cf898f1bc5e75f3fc2" dependencies = [ "proc-macro2", "quote", "unicode-ident", ] +[[package]] +name = "synstructure" +version = "0.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8af7666ab7b6390ab78131fb5b0fce11d6b7a6951602017c35fa82800708971" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.99", +] + [[package]] name = "tar" -version = "0.4.39" +version = "0.4.44" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ec96d2ffad078296368d46ff1cb309be1c23c513b4ab0e22a45de0185275ac96" +checksum = "1d863878d212c87a19c1a610eb53bb01fe12951c0501cf5a0d65f724914a667a" dependencies = [ "filetime", "libc", @@ -1545,33 +1820,32 @@ checksum = "e1fc403891a21bcfb7c37834ba66a547a8f402146eba7265b5a6d88059c9ff2f" [[package]] name = "thiserror" -version = "1.0.44" +version = "1.0.69" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "611040a08a0439f8248d1990b111c95baa9c704c805fa1f62104b39655fd7f90" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" dependencies = [ "thiserror-impl", ] [[package]] name = "thiserror-impl" -version = "1.0.44" +version = "1.0.69" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "090198534930841fab3a5d1bb637cde49e339654e606195f8d9c76eeb081dc96" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" dependencies = [ "proc-macro2", "quote", - "syn 2.0.28", + "syn 2.0.99", ] [[package]] -name = "time" -version = "0.1.45" +name = "tinystr" +version = "0.7.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1b797afad3f312d1c66a56d11d0316f916356d11bd158fbc6ca6389ff6bf805a" +checksum = "9117f5d4db391c1cf6927e7bea3db74b9a1c1add8f7eda9ffd5364f40f57b82f" dependencies = [ - "libc", - "wasi 0.10.0+wasi-snapshot-preview1", - "winapi", + "displaydoc", + "zerovec", ] [[package]] @@ -1585,40 +1859,44 @@ dependencies = [ ] [[package]] -name = "tinyvec" -version = "1.6.0" +name = "toml" +version = "0.8.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "87cc5ceb3875bb20c2890005a4e226a4651264a5c75edb2421b52861a0a0cb50" +checksum = "cd87a5cdd6ffab733b2f74bc4fd7ee5fff6634124999ac278c35fc78c6120148" dependencies = [ - "tinyvec_macros", + "serde", + "serde_spanned", + "toml_datetime", + "toml_edit", ] [[package]] -name = "tinyvec_macros" -version = "0.1.1" +name = "toml_datetime" +version = "0.6.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20" +checksum = "0dd7358ecb8fc2f8d014bf86f6f638ce72ba252a2c3a2572f2a795f1d23efb41" +dependencies = [ + "serde", +] [[package]] -name = "toml" -version = "0.5.11" +name = "toml_edit" +version = "0.22.24" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f4f7f0dd8d50a853a531c426359045b1998f04219d88799810762cd4ad314234" +checksum = "17b4795ff5edd201c7cd6dca065ae59972ce77d1b80fa0a84d94950ece7d1474" dependencies = [ + "indexmap", "serde", + "serde_spanned", + "toml_datetime", + "winnow", ] [[package]] name = "typenum" -version = "1.16.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "497961ef93d974e23eb6f433eb5fe1b7930b659f06d12dec6fc44a8f554c0bba" - -[[package]] -name = "unicode-bidi" -version = "0.3.13" +version = "1.18.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "92888ba5573ff080736b3648696b70cafad7d250551175acbaa4e0385b3e1460" +checksum = "1dccffe3ce07af9386bfd29e80c0ab1a8205a2fc34e4bcd40364df902cfa8f3f" [[package]] name = "unicode-ident" @@ -1626,15 +1904,6 @@ version = "1.0.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b15811caf2415fb889178633e7724bad2509101cde276048e013b9def5e51fa0" -[[package]] -name = "unicode-normalization" -version = "0.1.22" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5c5713f0fc4b5db668a2ac63cdb7bb4469d8c9fed047b1d0292cc7b0ce2ba921" -dependencies = [ - "tinyvec", -] - [[package]] name = "unicode-width" version = "0.1.14" @@ -1649,22 +1918,22 @@ checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" [[package]] name = "untrusted" -version = "0.7.1" +version = "0.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a156c684c91ea7d62626509bce3cb4e1d9ed5c4d978f7b4352658f96a4c26b4a" +checksum = "8ecb6da28b8a351d773b68d5825ac39017e680750f980f3a1a85cd8dd28a47c1" [[package]] name = "ureq" -version = "2.7.1" +version = "2.12.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0b11c96ac7ee530603dcdf68ed1557050f374ce55a5a07193ebf8cbc9f8927e9" +checksum = "02d1a66277ed75f640d608235660df48c8e3c19f3b4edb6a263315626cc3c01d" dependencies = [ - "base64 0.21.2", + "base64", "flate2", "log", "once_cell", "rustls", - "rustls-webpki 0.100.1", + "rustls-pki-types", "serde", "serde_json", "url", @@ -1673,15 +1942,27 @@ dependencies = [ [[package]] name = "url" -version = "2.4.0" +version = "2.5.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "50bff7831e19200a85b17131d085c25d7811bc4e186efdaf54bbd132994a88cb" +checksum = "32f8b686cadd1473f4bd0117a5d28d36b1ade384ea9b5069a1c40aefed7fda60" dependencies = [ "form_urlencoded", "idna", "percent-encoding", ] +[[package]] +name = "utf16_iter" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8232dd3cdaed5356e0f716d285e4b40b932ac434100fe9b7e0e8e935b9e6246" + +[[package]] +name = "utf8_iter" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6c140620e7ffbb22c2dee59cafe6084a59b5ffc27a8859a5f0d494b5d52b6be" + [[package]] name = "utf8parse" version = "0.2.2" @@ -1690,24 +1971,24 @@ checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" [[package]] name = "uuid" -version = "1.4.1" +version = "1.15.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "79daa5ed5740825c40b389c5e50312b9c86df53fccd33f281df655642b43869d" +checksum = "e0f540e3240398cce6128b64ba83fdbdd86129c16a3aa1a3a252efd66eb3d587" dependencies = [ - "getrandom", + "getrandom 0.3.1", ] [[package]] name = "version_check" -version = "0.9.4" +version = "0.9.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "49874b5167b65d7193b8aba1567f5c7d93d001cafc34600cee003eda787e483f" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "walkdir" -version = "2.3.3" +version = "2.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "36df944cda56c7d8d8b7496af378e6b16de9284591917d307c9b4d313c44e698" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" dependencies = [ "same-file", "winapi-util", @@ -1715,46 +1996,50 @@ dependencies = [ [[package]] name = "wasi" -version = "0.10.0+wasi-snapshot-preview1" +version = "0.11.0+wasi-snapshot-preview1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1a143597ca7c7793eff794def352d41792a93c481eb1042423ff7ff72ba2c31f" +checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" [[package]] name = "wasi" -version = "0.11.0+wasi-snapshot-preview1" +version = "0.13.3+wasi-0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" +checksum = "26816d2e1a4a36a2940b96c5296ce403917633dff8f3440e9b236ed6f6bacad2" +dependencies = [ + "wit-bindgen-rt", +] [[package]] name = "wasm-bindgen" -version = "0.2.87" +version = "0.2.100" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7706a72ab36d8cb1f80ffbf0e071533974a60d0a308d01a5d0375bf60499a342" +checksum = "1edc8929d7499fc4e8f0be2262a241556cfc54a0bea223790e71446f2aab1ef5" dependencies = [ "cfg-if", + "once_cell", + "rustversion", "wasm-bindgen-macro", ] [[package]] name = "wasm-bindgen-backend" -version = "0.2.87" +version = "0.2.100" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5ef2b6d3c510e9625e5fe6f509ab07d66a760f0885d858736483c32ed7809abd" +checksum = "2f0a0651a5c2bc21487bde11ee802ccaf4c51935d0d3d42a6101f98161700bc6" dependencies = [ "bumpalo", "log", - "once_cell", "proc-macro2", "quote", - "syn 2.0.28", + "syn 2.0.99", "wasm-bindgen-shared", ] [[package]] name = "wasm-bindgen-macro" -version = "0.2.87" +version = "0.2.100" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dee495e55982a3bd48105a7b947fd2a9b4a8ae3010041b9e0faab3f9cd028f1d" +checksum = "7fe63fc6d09ed3792bd0897b314f53de8e16568c2b3f7982f468c0bf9bd0b407" dependencies = [ "quote", "wasm-bindgen-macro-support", @@ -1762,22 +2047,25 @@ dependencies = [ [[package]] name = "wasm-bindgen-macro-support" -version = "0.2.87" +version = "0.2.100" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "54681b18a46765f095758388f2d0cf16eb8d4169b639ab575a8f5693af210c7b" +checksum = "8ae87ea40c9f689fc23f209965b6fb8a99ad69aeeb0231408be24920604395de" dependencies = [ "proc-macro2", "quote", - "syn 2.0.28", + "syn 2.0.99", "wasm-bindgen-backend", "wasm-bindgen-shared", ] [[package]] name = "wasm-bindgen-shared" -version = "0.2.87" +version = "0.2.100" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ca6ad05a4870b2bf5fe995117d3728437bd27d7cd5f06f13c17443ef369775a1" +checksum = "1a05d73b933a847d6cccdda8f838a22ff101ad9bf93e33684f39c1f5f0eece3d" +dependencies = [ + "unicode-ident", +] [[package]] name = "web-sys" @@ -1791,11 +2079,11 @@ dependencies = [ [[package]] name = "webpki-roots" -version = "0.23.1" +version = "0.26.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b03058f88386e5ff5310d9111d53f48b17d732b401aeb83a8d5190f2ac459338" +checksum = "2210b291f7ea53617fbafcc4939f10914214ec15aace5ba62293a668f322c5c9" dependencies = [ - "rustls-webpki 0.100.1", + "rustls-pki-types", ] [[package]] @@ -1830,14 +2118,20 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" [[package]] -name = "windows" -version = "0.48.0" +name = "windows-core" +version = "0.52.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e686886bc078bc1b0b600cac0147aadb815089b6e4da64016cbd754b6342700f" +checksum = "33ab640c8d7e35bf8ba19b884ba838ceb4fba93a4e8c65a9059d08afcfc683d9" dependencies = [ - "windows-targets 0.48.0", + "windows-targets 0.52.5", ] +[[package]] +name = "windows-link" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6dccfd733ce2b1753b03b6d3c65edf020262ea35e20ccdf3e288043e6dd620e3" + [[package]] name = "windows-sys" version = "0.48.0" @@ -1977,11 +2271,115 @@ version = "0.52.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bec47e5bfd1bff0eeaf6d8b485cc1074891a197ab4225d504cb7a1ab88b02bf0" +[[package]] +name = "winnow" +version = "0.7.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0e7f4ea97f6f78012141bcdb6a216b2609f0979ada50b20ca5b52dde2eac2bb1" +dependencies = [ + "memchr", +] + +[[package]] +name = "wit-bindgen-rt" +version = "0.33.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3268f3d866458b787f390cf61f4bbb563b922d091359f9608842999eaee3943c" +dependencies = [ + "bitflags 2.9.0", +] + +[[package]] +name = "write16" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d1890f4022759daae28ed4fe62859b1236caebfc61ede2f63ed4e695f3f6d936" + +[[package]] +name = "writeable" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e9df38ee2d2c3c5948ea468a8406ff0db0b29ae1ffde1bcf20ef305bcc95c51" + [[package]] name = "xattr" -version = "0.2.3" +version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6d1526bbe5aaeb5eb06885f4d987bcdfa5e23187055de9b83fe00156a821fabc" +checksum = "0d65cbf2f12c15564212d48f4e3dfb87923d25d611f2aed18f4cb23f0413d89e" dependencies = [ "libc", + "rustix 1.0.1", +] + +[[package]] +name = "yoke" +version = "0.7.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "120e6aef9aa629e3d4f52dc8cc43a015c7724194c97dfaf45180d2daf2b77f40" +dependencies = [ + "serde", + "stable_deref_trait", + "yoke-derive", + "zerofrom", +] + +[[package]] +name = "yoke-derive" +version = "0.7.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2380878cad4ac9aac1e2435f3eb4020e8374b5f13c296cb75b4620ff8e229154" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.99", + "synstructure", +] + +[[package]] +name = "zerofrom" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "50cc42e0333e05660c3587f3bf9d0478688e15d870fab3346451ce7f8c9fbea5" +dependencies = [ + "zerofrom-derive", +] + +[[package]] +name = "zerofrom-derive" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d71e5d6e06ab090c67b5e44993ec16b72dcbaabc526db883a360057678b48502" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.99", + "synstructure", +] + +[[package]] +name = "zeroize" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ced3678a2879b30306d323f4542626697a464a97c0a07c9aebf7ebca65cd4dde" + +[[package]] +name = "zerovec" +version = "0.10.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa2b893d79df23bfb12d5461018d408ea19dfafe76c2c7ef6d4eba614f8ff079" +dependencies = [ + "yoke", + "zerofrom", + "zerovec-derive", +] + +[[package]] +name = "zerovec-derive" +version = "0.10.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6eafa6dfb17584ea3e2bd6e76e0cc15ad7af12b09abdd1ca55961bed9b1063c6" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.99", ] diff --git a/Cargo.toml b/Cargo.toml index e3d7cf3..b5ca202 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,17 +15,17 @@ csv = "1.2.2" derive_more = "0.99.16" num = "0.4.0" # Wrappers crates -lapack = "0.19.0" -blas = "0.22.0" +lapack = "0.20.0" +blas = "0.23.0" intel-mkl-tool = "0.8.1" # Bindings crates -lapack-sys = "0.14.0" -blas-sys = "0.7.1" +lapack-sys = "0.15.0" +blas-sys = "0.8.0" # Source crates # Intel implementation. intel-mkl-src = {version = "0.8.1", features = ["mkl-dynamic-lp64-iomp"]} # Openblas implementation. -# openblas-src = "0.10.8" +#openblas-src = "0.10.8" # For python interface pyo3 = {version = "0.21.2", features = ["extension-module"], optional = true} rayon = "1.10.0" @@ -33,6 +33,8 @@ rayon = "1.10.0" rand_mt = {version = "4.2.2", features = ["rand_core"]} # Progress bar indicatif = "0.17.8" +# Colored terminal text +colored = "3.0.0" [dependencies.rand] version = "0.8.5" diff --git a/src/bin/boboche2000.rs b/src/bin/boboche2000.rs index 057b3d0..14d1995 100644 --- a/src/bin/boboche2000.rs +++ b/src/bin/boboche2000.rs @@ -79,6 +79,7 @@ fn main() { nvij: NVIJ, ngi: NGI, mcsample_interval: 1, + nbootstrap: 1, transfert_matrix: &HOPPINGS, hopping_bitmask: &bitmask, clean_update_frequency: CLEAN_UPDATE_FREQUENCY, @@ -86,6 +87,7 @@ fn main() { nmcsample: NMCSAMP, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERANCE_SINGULARITY, + pair_wavefunction: false, }; log_system_parameters(&sys); @@ -170,7 +172,7 @@ fn main() { info!("Initial Nelec: {}, {}", state.spin_down.count_ones(), state.spin_up.count_ones()); info!("Nsites: {}", state.n_sites); - let (energy, _, _) = compute_mean_energy(&mut rng, state, ¶meters, &sys, &mut der); + let (energy, _, _, _) = compute_mean_energy(&mut rng, state, ¶meters, &sys, &mut der); println!("energy: {}", energy); //close(energy, -0.35, MONTE_CARLO_CONVERGENCE_TOLERANCE); //close(energy, mean_energy_analytic_2sites(¶meters, &sys), MONTE_CARLO_CONVERGENCE_TOLERANCE); diff --git a/src/bin/dvmc_fastupdate.rs b/src/bin/dvmc_fastupdate.rs index 02f9e1e..f979ab5 100644 --- a/src/bin/dvmc_fastupdate.rs +++ b/src/bin/dvmc_fastupdate.rs @@ -5,20 +5,21 @@ use indicatif::{ProgressBar, ProgressStyle}; use std::fs::File; use std::io::Write; -use impurity::optimisation::conjugate_gradiant; +use impurity::optimisation::{exact_overlap_inverse, spread_eigenvalues}; use impurity::{generate_bitmask, DerivativeOperator, FockState, RandomStateGeneration, SysParams, VarParams}; use impurity::monte_carlo::compute_mean_energy; -const SEED: u64 = 1434; +const SEED: u64 = 14345; const SIZE: usize = 2; const NFIJ: usize = 4*SIZE*SIZE; -const NVIJ: usize = SIZE*SIZE; +const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; const NMCSAMP: usize = 10_000; +const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; -const MCSAMPLE_INTERVAL: usize = 2; +const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 6; const CLEAN_UPDATE_FREQUENCY: usize = 2; const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; @@ -26,9 +27,11 @@ const TOLERENCE_SINGULARITY: f64 = 1e-12; const CONS_U: f64 = 1.0; const CONS_T: f64 = -1.0; const EPSILON_CG: f64 = 1e-16; -const EPSILON_SPREAD: f64 = 0.0; +const EPSILON_SPREAD: f64 = 0.01; const OPTIMISATION_TIME_STEP: f64 = 1e-2; const NOPTITER: usize = 1_000; +const _KMAX: usize = NMCSAMP; +const PARAMTHRESHOLD: f64 = 0.5; pub const HOPPINGS: [f64; SIZE*SIZE] = [ 0.0, 1.0, 1.0, 0.0, @@ -57,7 +60,7 @@ fn norm(par: &VarParams) -> f64 { let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; let g0 = par.gi[0]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let a = ::exp(2.0 * g0 - 2.0 * v)*sq(::abs(f00ud - f00du)); let b = ::exp(2.0 * g1 - 2.0 * v)*sq(::abs(f11ud - f11du)); let c = sq(::abs(f01uu - f10uu)); @@ -74,17 +77,17 @@ fn mean_energy_analytic_2sites(par: &VarParams, _sys: &SysParams) -> f64 { - par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; let b = (par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) - * ::exp(par.gi[0]-par.vij[1]); + * ::exp(par.gi[0]-par.vij[0]); let c = (par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) - * ::exp(par.gi[1]-par.vij[1]); + * ::exp(par.gi[1]-par.vij[0]); let d = 2.0 * CONS_T * (b + c) * a; let e = sq(par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) - * ::exp(2.0*par.gi[1]-2.0*par.vij[1]) * CONS_U; + * ::exp(2.0*par.gi[1]-2.0*par.vij[0]) * CONS_U; let f = sq(par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) - * ::exp(2.0*par.gi[0]-2.0*par.vij[1]) * CONS_U; + * ::exp(2.0*par.gi[0]-2.0*par.vij[0]) * CONS_U; (d + e + f) / norm(par) } @@ -108,7 +111,7 @@ fn log_system_parameters(e: f64, ae: f64, corr_time: f64, fp: &mut File, params: let a = format!("{:+>.05e} ", fij[i]); params.push_str(&a); } - for i in 0..SIZE*SIZE { + for i in 0..NVIJ { let a = format!("{:+>.05e} ", vij[i]); params.push_str(&a); } @@ -136,6 +139,7 @@ fn zero_out_derivatives(der: &mut DerivativeOperator) { fn main() { let mut fp = File::create("params").unwrap(); + writeln!(fp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); let mut statesfp = File::create("states").unwrap(); let mut save: bool = true; // Initialize logger @@ -154,14 +158,20 @@ fn main() { hopping_bitmask: &bitmask, clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcsample: NMCSAMP, + nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, - tolerance_singularity: TOLERENCE_SINGULARITY + tolerance_singularity: TOLERENCE_SINGULARITY, + pair_wavefunction: false, }; let mut rng = Mt64::new(SEED); //let parameters = generate_random_params(&mut rng); + //let mut all_params = Vec::with_capacity(NGI + NVIJ + NFIJ); + //for i in 0..(NGI + NVIJ + NFIJ) { + // all_params.push(rng.gen()); + //} //let mut all_params = vec![ // -1.0, -1.0, -1.0, -1.0, // 0.0, 2.0, 2.0, 2.0, @@ -227,10 +237,7 @@ fn main() { let mut all_params = vec![ 2.992823494391085859e-01, -8.118842052166648227e-01, - 0.0, -5.126018557775564588e-01, - -5.126018557775564588e-01, - 0.0, //0.000000000000000000e+00, 0.0, 0.0, 0.0, 0.0, 1.085729148576013436e-01, @@ -292,7 +299,7 @@ fn main() { .progress_chars("##-")); for opt_iter in 0..NOPTITER { - let (mean_energy, accumulated_states, correlation_time) = { + let (mean_energy, accumulated_states, correlation_time, _) = { compute_mean_energy(&mut rng, state, ¶meters, &system_params, &mut derivative) }; if save { @@ -313,21 +320,19 @@ fn main() { unsafe { let incx = 1; let incy = 1; - dscal(derivative.n, 1.0 / (NMCSAMP as f64), derivative.ho, incx); - dscal(derivative.n, 1.0 / (NMCSAMP as f64), derivative.expval_o, incx); daxpy(derivative.n, -mean_energy, derivative.expval_o, incx, derivative.ho, incy); dcopy(derivative.n, derivative.ho, incx, &mut b, incy); } - //spread_eigenvalues(&mut derivative); - conjugate_gradiant(&derivative, &mut b, &mut x0, EPSILON_CG, 4, NPARAMS as i32); + spread_eigenvalues(&mut derivative); + exact_overlap_inverse(&derivative, &mut b, EPSILON_CG, NPARAMS as i32, PARAMTHRESHOLD); info!("Need to update parameters with: {:?}", x0); unsafe { let incx = 1; let incy = 1; let alpha = - OPTIMISATION_TIME_STEP; - daxpy(NGI as i32, alpha, &x0, incx, &mut parameters.gi, incy); - daxpy(NVIJ as i32, alpha, &x0[NGI..NPARAMS], incx, &mut parameters.vij, incy); - daxpy(NFIJ as i32, alpha, &x0[NGI + NVIJ..NPARAMS], incx, &mut parameters.fij, incy); + daxpy(NGI as i32, alpha, &b, incx, &mut parameters.gi, incy); + daxpy(NVIJ as i32, alpha, &b[NGI..NPARAMS], incx, &mut parameters.vij, incy); + daxpy(NFIJ as i32, alpha, &b[NGI + NVIJ..NPARAMS], incx, &mut parameters.fij, incy); } info!("Correctly finished optimisation iteration {}", opt_iter); // Sorella Louche stuff @@ -361,7 +366,6 @@ fn main() { parameters.vij[i] -= shift; } // HARD CODE vij = vji - parameters.vij[1] = parameters.vij[2]; // Slater Rescaling unsafe { let incx = 1; diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs new file mode 100644 index 0000000..3e0af27 --- /dev/null +++ b/src/bin/dvmc_pairwf.rs @@ -0,0 +1,460 @@ +use blas::{daxpy, dcopy, dgemm, dnrm2, dscal, idamax}; +use log::{debug, error, info}; +use rand_mt::Mt64; +use indicatif::{ProgressBar, ProgressStyle}; +use std::fs::File; +use std::io::Write; + +use impurity::optimisation::{conjugate_gradiant, exact_overlap_inverse}; +use impurity::{generate_bitmask, mapto_pairwf, DerivativeOperator, FockState, RandomStateGeneration, SysParams, VarParams}; +use impurity::monte_carlo::{compute_mean_energy, compute_mean_energy_exact}; + +const SEED: u64 = 14; +const SIZE: usize = 2; +const NFIJ: usize = SIZE*SIZE; +const NVIJ: usize = SIZE*(SIZE - 1) / 2; +const NGI: usize = SIZE; +const NPARAMS: usize = NFIJ + NGI + NVIJ; +const NELEC: usize = SIZE; +const NMCSAMP: usize = 1_000; +const NBOOTSTRAP: usize = 1; +const NMCWARMUP: usize = 1000; +const MCSAMPLE_INTERVAL: usize = 1; +const _NTHREADS: usize = 1; +const CLEAN_UPDATE_FREQUENCY: usize = 32; +const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; +const TOLERENCE_SINGULARITY: f64 = 1e-12; +const CONS_U: f64 = 8.0; +const CONS_T: f64 = -1.0; +const EPSILON_CG: f64 = 1e-16; +const EPSILON_SHIFT: f64 = 1e-2; +const OPTIMISATION_TIME_STEP: f64 = 1e-2; +const OPTIMISATION_DECAY: f64 = 0.0; +const NOPTITER: usize = 1000; +const KMAX: usize = NPARAMS; +const PARAM_THRESHOLD: f64 = ::MIN; +//const PARAM_THRESHOLD: f64 = 0.0; +const OPTIMISE: bool = true; +const OPTIMISE_GUTZ: bool = true; +const OPTIMISE_JAST: bool = true; +const OPTIMISE_ORB: bool = true; +const SET_EXPVALO_ZERO: bool = false; +const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo; +const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ExactInverse; + +pub enum EnergyOptimisationMethod { + ExactInverse, + ConjugateGradiant, +} + +pub enum EnergyComputationMethod { + MonteCarlo, + ExactSum, +} + +pub const HOPPINGS: [f64; SIZE*SIZE] = [ + 0.0, 1.0, 1.0, 0.0, + //0.0, 1.0, 1.0, 0.0, + //1.0, 0.0, 0.0, 1.0, + //1.0, 0.0, 0.0, 1.0, + //0.0, 1.0, 1.0, 0.0 +]; + +fn _print_delta_alpha(da: &[f64], ngi: usize, nvij: usize, nfij: usize) { + let mut outstr = "".to_owned(); + for i in 0..ngi { + outstr.push_str(&format!(" G_{} = {}", i, da[i])); + } + for i in ngi..ngi+nvij { + outstr.push_str(&format!(" V_{} = {}", i, da[i])); + } + for i in ngi+nvij..ngi+nvij+nfij { + outstr.push_str(&format!(" F_{} = {}", i, da[i])); + } + println!("{}", outstr); +} + +fn _save_otilde(fp: &mut File, der: &DerivativeOperator) { + let width = 16; + let mut o_tilde = "".to_owned(); + for mu in 0..(der.mu + 1) as usize { + for n in 0..der.n as usize { + o_tilde.push_str(&format!("{:>width$.04e}", der.o_tilde[n + mu * der.n as usize])); + } + o_tilde.push_str("\n"); + } + fp.write(&o_tilde.as_bytes()).unwrap(); + let mut c = vec![0.0; (der.n * der.n) as usize]; + println!("dim = {}", der.n * der.n); + unsafe { + dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, der.o_tilde, der.n, der.o_tilde, der.n, 0.0, &mut c, der.n); + } + let mut outstr = "".to_owned(); + outstr.push_str(&format!(" = ")); + for i in 0..der.n as usize { + outstr.push_str(&format!("\n ")); + for j in 0..der.n as usize { + outstr.push_str(&format!("{:>width$.04e}", c[i + der.n as usize * j])); + } + } + println!("{}", outstr); +} + +fn sq(a: f64) -> f64 { + ::abs(a) * ::abs(a) +} + +fn norm(par: &VarParams) -> f64 { + let f00ud = par.fij[0 + 0 * SIZE + SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + SIZE * SIZE]; + let f01ud = par.fij[0 + 1 * SIZE + SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = ::exp(2.0 * g0 - 2.0 * v)*sq(::abs(f00ud)); + let b = ::exp(2.0 * g1 - 2.0 * v)*sq(::abs(f11ud)); + let e = sq(::abs(f10ud)); + let f = sq(::abs(f01ud)); + a + b + e + f +} + +fn mean_energy_analytic_2sites(par: &VarParams, _sys: &SysParams) -> f64 { + let a = par.fij[1 + 0 * SIZE + SIZE * SIZE] + + par.fij[0 + 1 * SIZE + SIZE * SIZE]; + let b = par.fij[0 + 0 * SIZE + SIZE * SIZE] + * ::exp(par.gi[0]-par.vij[0]); + let c = par.fij[1 + 1 * SIZE + SIZE * SIZE] + * ::exp(par.gi[1]-par.vij[0]); + let d = 2.0 * CONS_T * (b + c) * a; + let e = sq(par.fij[1 + 1 * SIZE + SIZE * SIZE]) + * ::exp(2.0*par.gi[1]-2.0*par.vij[0]) * CONS_U; + let f = sq(par.fij[0 + 0 * SIZE + SIZE * SIZE]) + * ::exp(2.0*par.gi[0]-2.0*par.vij[0]) * CONS_U; + (d + e + f) / norm(par) +} + +fn log_system_parameters(e: f64, ae: f64, deltae: f64, corr_time: f64, fp: &mut File, params: &VarParams, sys: &SysParams, dpar: &[f64]) { + let fij = ¶ms.fij; + let vij = ¶ms.vij; + let gi = ¶ms.gi; + info!("System parameter SIZE = {}", SIZE); + info!("System parameter NELEC = {}", NELEC); + info!("System parameter NMCSAMP = {}", NMCSAMP); + info!("System parameter NMCWARMUP = {}", NMCWARMUP); + info!("System parameter CONS_U = {}", sys.cons_u); + info!("System parameter CONS_T = {}", sys.cons_t); + debug!("System parameter CLEAN_UPDATE_FREQUENCY = {}", CLEAN_UPDATE_FREQUENCY); + debug!("System parameter TOLERENCE_SHERMAN_MORRISSON = {}", TOLERENCE_SHERMAN_MORRISSON); + debug!("\n{}", params); + let mut params = format!("{:+>.05e} ", e).to_owned(); + params.push_str(&format!("{:+>.05e} ", ae).to_owned()); + params.push_str(&format!("{:+>.05e} ", deltae).to_owned()); + params.push_str(&format!("{:+>.05e} ", corr_time).to_owned()); + for i in 0..SIZE { + let a = format!("{:+>.05e} ", gi[i]); + params.push_str(&a); + } + for i in 0..NVIJ { + let a = format!("{:+>.05e} ", vij[i]); + params.push_str(&a); + } + for i in 0..4 { + let a = format!("{:+>.05e} ", fij[4+i]); + params.push_str(&a); + } + for i in 0..7 { + let a = format!("{:+>.05e} ", dpar[i]); + params.push_str(&a); + } + params.push_str(&"\n"); + fp.write(params.as_bytes()).unwrap(); +} + +fn zero_out_derivatives(der: &mut DerivativeOperator) { + for i in 0.. (NFIJ + NVIJ + NGI) * NMCSAMP { + der.o_tilde[i] = 0.0; + } + for i in 0..NFIJ + NVIJ + NGI { + der.expval_o[i] = 0.0; + der.ho[i] = 0.0; + } + for i in 0..NMCSAMP { + der.visited[i] = 0; + } + der.mu = -1; +} + +fn main() { + let mut fp = File::create("params").unwrap(); + writeln!(fp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); + let mut _save: bool = true; + // Initialize logger + env_logger::init(); + let bitmask = generate_bitmask(&HOPPINGS, SIZE); + let system_params = SysParams { + size: SIZE, + nelec: NELEC, + array_size: (SIZE + 7) / 8, + cons_t: CONS_T, + cons_u: CONS_U, + nfij: NFIJ, + nvij: NVIJ, + ngi: NGI, + transfert_matrix: &HOPPINGS, + hopping_bitmask: &bitmask, + clean_update_frequency: CLEAN_UPDATE_FREQUENCY, + nmcsample: NMCSAMP, + nbootstrap: NBOOTSTRAP, + nmcwarmup: NMCWARMUP, + mcsample_interval: MCSAMPLE_INTERVAL, + tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, + tolerance_singularity: TOLERENCE_SINGULARITY, + pair_wavefunction: true, + }; + + let mut rng = Mt64::new(SEED); + //let parameters = generate_random_params(&mut rng); + //let mut all_params = Vec::with_capacity(NGI + NVIJ + NFIJ); + //for _ in 0..(NGI + NVIJ + NFIJ) { + // all_params.push(rng.gen()); + //} + //let (gi, params) = all_params.split_at_mut(NGI); + //let (vij, mut fij) = params.split_at_mut(NVIJ); + // 0.000000000000000000e+00 5.012713072533996646e-09 + let mut fij = [0.5, 0.5, 0.5, 0.5]; + let mut general_fij = [ + 0.0, 0.0, 0.0, 0.0, + //1.093500753438337580e-01, 3.768419990611672210e-01, 3.769186909982900624e-01, 3.322533463612635796e-01, + 0.5, 0.5, 0.5, 0.5, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + ]; + let mut vij = [0.0]; + let mut gi = [0.0, 0.0]; + //let mut general_fij: Vec = vec![0.0; 4*NFIJ]; + let mut parameters = VarParams { + fij: &mut general_fij, + gi: &mut gi, + vij: &mut vij, + size: SIZE + }; + unsafe { + dcopy(NFIJ as i32, &fij[0..NFIJ], 1, &mut parameters.fij[NFIJ..2*NFIJ], 1); + } + //println!("{}", mean_energy_analytic_2sites(¶meters, &system_params)); + //log_system_parameters(0.0, &mut fp, ¶meters, &system_params); + + let state: FockState = { + let mut tmp: FockState = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + while tmp.spin_up.count_ones() != tmp.spin_down.count_ones() { + tmp = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + } + tmp + }; + + let mut otilde: Vec = vec![0.0; (4*NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; + let mut work_otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; + let mut expvalo: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; + let mut work_expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut expval_ho: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; + let mut work_expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut visited: Vec = vec![0; NMCSAMP + 1]; + let mut work_visited: Vec = vec![0; NMCSAMP + 1]; + let mut derivative = DerivativeOperator { + o_tilde: &mut otilde, + expval_o: &mut expvalo, + ho: &mut expval_ho, + n: (4*NFIJ + NVIJ + NGI) as i32, + nsamp: match COMPUTE_ENERGY_METHOD { + EnergyComputationMethod::ExactSum => 1.0, + EnergyComputationMethod::MonteCarlo => NMCSAMP as f64, + }, + nsamp_int: MCSAMPLE_INTERVAL, + mu: -1, + visited: &mut visited, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + epsilon: EPSILON_SHIFT, + }; + let mut work_derivative = DerivativeOperator { + o_tilde: &mut work_otilde, + expval_o: &mut work_expvalo, + ho: &mut work_expval_ho, + n: (NFIJ + NVIJ + NGI) as i32, + nsamp: match COMPUTE_ENERGY_METHOD { + EnergyComputationMethod::ExactSum => 1.0, + EnergyComputationMethod::MonteCarlo => NMCSAMP as f64, + }, + nsamp_int: MCSAMPLE_INTERVAL, + mu: -1, + visited: &mut work_visited, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + epsilon: EPSILON_SHIFT, + }; + + info!("Initial State: {}", state); + info!("Initial Nelec: {}, {}", state.spin_down.count_ones(), state.spin_up.count_ones()); + info!("Nsites: {}", state.n_sites); + + let opt_progress_bar = ProgressBar::new(NOPTITER as u64); + opt_progress_bar.set_prefix("Optimisation Progress: "); + opt_progress_bar.set_style(ProgressStyle::with_template("[{elapsed_precise}] {prefix} {bar:40.cyan/blue} {pos:>7}/{len:7} {msg}") + .unwrap() + .progress_chars("##-")); + + for opt_iter in 0..NOPTITER { + let (mean_energy, accumulated_states, deltae, correlation_time) = { + match COMPUTE_ENERGY_METHOD { + EnergyComputationMethod::MonteCarlo => compute_mean_energy(&mut rng, state, ¶meters, &system_params, &mut derivative), + EnergyComputationMethod::ExactSum => { + (compute_mean_energy_exact(¶meters, &system_params, &mut derivative), Vec::with_capacity(0), 0.0, 0.0) + }, + } + }; + if true { + let mut statesfp = File::create("states").unwrap(); + let mut out_str: String = String::new(); + for s in accumulated_states.iter() { + out_str.push_str(&format!("{}\n", s)); + } + statesfp.write(out_str.as_bytes()).unwrap(); + _save = false; + } + // Copy out the relevant terms. + work_derivative.mu = derivative.mu; + let mut i = 0; + for elem in derivative.visited.iter() { + work_derivative.visited[i] = *elem; + i += 1; + } + mapto_pairwf(&derivative, &mut work_derivative, &system_params); + + let mut x0 = vec![0.0; NFIJ + NVIJ + NGI]; + x0[(NGI + NVIJ)..(NGI + NVIJ + NFIJ)].copy_from_slice(&fij); + x0[NGI..(NGI + NVIJ)].copy_from_slice(parameters.vij); + x0[0..NGI].copy_from_slice(parameters.gi); + + if SET_EXPVALO_ZERO { + for i in 0..work_derivative.n as usize { + work_derivative.expval_o[i] = 0.0; + } + } + + // 68 misawa + let mut b: Vec = vec![0.0; work_derivative.n as usize]; + unsafe { + let incx = 1; + let incy = 1; + daxpy(work_derivative.n, -mean_energy, work_derivative.expval_o, incx, work_derivative.ho, incy); + dcopy(work_derivative.n, work_derivative.ho, incx, &mut b, incy); + } + //save_otilde(&mut der_fp, &derivative); + //save_otilde(&mut wder_fp, &work_derivative); + //println!(" = {:?}", work_derivative.ho); + //println!(" = {:?}", work_derivative.expval_o); + //spread_eigenvalues(&mut derivative); + //println!("x0 in = {:?}", x0); + let mut _flag: bool = true; + let ignored_columns = match OPTIMISE_ENERGY_METHOD { + EnergyOptimisationMethod::ExactInverse => { + exact_overlap_inverse(&work_derivative, &mut b, EPSILON_SHIFT, NPARAMS as i32, PARAM_THRESHOLD) + }, + EnergyOptimisationMethod::ConjugateGradiant => { + conjugate_gradiant(&work_derivative, &mut b, &mut x0, EPSILON_SHIFT, KMAX, NPARAMS as i32, PARAM_THRESHOLD, EPSILON_CG) + }, + }; + let mut delta_alpha = vec![0.0; NPARAMS]; + let mut j: usize = 0; + for i in 0..NPARAMS { + if ignored_columns[i] { + continue; + } + delta_alpha[i] = b[j]; + j += 1; + if !::is_finite(delta_alpha[i]) { + _flag = false; + } + } + let analytic_energy = mean_energy_analytic_2sites(¶meters, &system_params); + if OPTIMISE { + unsafe { + let incx = 1; + let incy = 1; + let alpha = - OPTIMISATION_TIME_STEP * ::exp(- (opt_iter as f64) * OPTIMISATION_DECAY); + if OPTIMISE_GUTZ { + daxpy(NGI as i32, alpha, &delta_alpha, incx, &mut parameters.gi, incy); + } + if OPTIMISE_JAST { + daxpy(NVIJ as i32, alpha, &delta_alpha[NGI..NPARAMS], incx, &mut parameters.vij, incy); + } + if OPTIMISE_ORB { + daxpy(NFIJ as i32, alpha, &delta_alpha[NGI + NVIJ..NPARAMS], incx, &mut fij, incy); + } + } + info!("Correctly finished optimisation iteration {}", opt_iter); + //info!("Rescaling the parameters."); + //let scale: f64 = unsafe { + // let incx = 1; + // let incy = 1; + // ddot(derivative.n, derivative.expval_o, incx, parameters.gi, incy) + //}; + //info!("Scale by : {}", scale); + //let ratio = 1.0 / (scale + 1.0); + //unsafe { + // let incx = 1; + // dscal(NPARAMS as i32, ratio, parameters.gi, incx) + //} + //info!("Scaled parameters by ratio = {}", ratio); + + // JastrowGutzwiller Shifting + let mut shift = 0.0; + for i in 0..NGI { + shift += parameters.gi[i]; + } + for i in 0..NVIJ { + shift += parameters.vij[i]; + } + shift = shift / (NGI + NVIJ) as f64; + for i in 0..NGI { + parameters.gi[i] -= shift; + } + for i in 0..NVIJ { + parameters.vij[i] -= shift; + } + } + // HARD CODE vij = vji + // Slater Rescaling + unsafe { + let incx = 1; + let max = fij[idamax(NFIJ as i32, &fij, incx) - 1]; + if ::abs(max) < 1e-16 { + error!("Pfaffian parameters are all close to 0.0. Rescaling might bring noise."); + panic!("Undefined behavior."); + } + info!("Max was: {}", max); + dscal(NFIJ as i32, 1.0 / max, &mut fij, incx); + } + unsafe { + dcopy( + NFIJ as i32, + &fij, + 1, + &mut parameters.fij[NFIJ..2*NFIJ], + 1 + ); + } + log_system_parameters(mean_energy, analytic_energy, deltae, correlation_time, &mut fp, ¶meters, &system_params, &delta_alpha); + zero_out_derivatives(&mut derivative); + //print_delta_alpha(&delta_alpha, NGI, NVIJ, NFIJ); + let opt_delta = unsafe { + let incx = 1; + dnrm2(work_derivative.n, &delta_alpha, incx) + }; + //error!("Changed parameters by norm {}", opt_delta); + opt_progress_bar.inc(1); + opt_progress_bar.set_message(format!("Changed parameters by norm: {:+>.05e} Current energy: {:+>.05e}", opt_delta, mean_energy)); + } + opt_progress_bar.finish() +} diff --git a/src/constants.rs b/src/constants.rs index b506ce5..08256ad 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -16,10 +16,12 @@ pub struct SysParams<'a> { pub hopping_bitmask: &'a [SpinState], pub nmcwarmup: usize, pub nmcsample: usize, + pub nbootstrap: usize, pub mcsample_interval: usize, pub clean_update_frequency: usize, pub tolerance_singularity: f64, pub tolerance_sherman_morrison: f64, + pub pair_wavefunction: bool, } pub type BitStruct = u8; diff --git a/src/density.rs b/src/density.rs index 9373ee7..3f5587e 100644 --- a/src/density.rs +++ b/src/density.rs @@ -5,17 +5,18 @@ use pyo3::{pyfunction, PyResult}; use crate::jastrow::{compute_jastrow_exp, fast_update_jastrow}; use crate::gutzwiller::{compute_gutzwiller_exp, fast_update_gutzwiller}; use crate::pfaffian::{construct_matrix_a_from_state, get_pfaffian_ratio, PfaffianState}; -use crate::{BitOps, FockState, Spin, SpinState, VarParams}; +use crate::{BitOps, FockState, Spin, SpinState, VarParams, SysParams}; pub fn compute_internal_product( state: FockState, params: &VarParams, + sys: &SysParams, ) -> f64 where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl { let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( - construct_matrix_a_from_state(params.fij, state), + construct_matrix_a_from_state(params.fij, state, sys), compute_jastrow_exp(state, params.vij, state.n_sites), compute_gutzwiller_exp(state, params.gi, state.n_sites) ) @@ -31,12 +32,13 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl pub fn compute_internal_product_parts( state: FockState, params: &VarParams, + sys: &SysParams, ) -> (PfaffianState, f64) where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl { let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( - construct_matrix_a_from_state(params.fij, state), + construct_matrix_a_from_state(params.fij, state, sys), compute_jastrow_exp(state, params.vij, state.n_sites), compute_gutzwiller_exp(state, params.gi, state.n_sites) ) diff --git a/src/fock_state.rs b/src/fock_state.rs index 96bfff0..4200ff9 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -556,11 +556,14 @@ impl Distribution> for Standard wher pub trait RandomStateGeneration { fn generate_from_nelec(rng: &mut R, nelec: usize, max_size: usize) -> Self; - fn generate_hopping(self: &Self, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin), sys: &SysParams) -> Self; + fn generate_hopping(self: &Self, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin)) -> Self; fn generate_exchange(self: &Self, rng: &mut R, out_idx: &mut (usize, usize)) -> Self; } -impl> RandomStateGeneration for FockState where Standard: Distribution { +impl RandomStateGeneration for FockState +where T: BitOps, + Standard: Distribution +{ fn generate_from_nelec(rng: &mut R, nelec: usize, max_size: usize) -> FockState { let mut state = FockState{spin_up: ::zeros(), spin_down: ::zeros(), n_sites: max_size}; let mut i = 0; @@ -584,30 +587,60 @@ impl> RandomStateGeneration for FockState where Standard } - fn generate_hopping(self: &FockState, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin), sys: &SysParams) -> FockState { + fn generate_hopping(self: &FockState, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin)) -> FockState { // This is cheap, don't sweat it - let all_hops = self.generate_all_hoppings(&sys.hopping_bitmask); + //let all_hops = self.generate_all_hoppings(&sys.hopping_bitmask); + + let spin = if rng.gen::() % 2 == 0 { + Spin::Up + } else { + Spin::Down + }; + let mut rand_from_idx; + let mut rand_to_idx; + match spin { + Spin::Up => { + // Choose spin from + rand_from_idx = (rng.gen::() % max_size) as usize; + while !self.spin_up.check(rand_from_idx) { + rand_from_idx = (rng.gen::() % max_size) as usize; + } + rand_to_idx = (rng.gen::() % max_size) as usize; + while self.spin_up.check(rand_to_idx) { + rand_to_idx = (rng.gen::() % max_size) as usize; + } + }, + Spin::Down => { + // Choose spin from + rand_from_idx = (rng.gen::() % max_size) as usize; + while !self.spin_down.check(rand_from_idx) { + rand_from_idx = (rng.gen::() % max_size) as usize; + } + rand_to_idx = (rng.gen::() % max_size) as usize; + while self.spin_down.check(rand_to_idx) { + rand_to_idx = (rng.gen::() % max_size) as usize; + } + } + } - let rand_index = rng.gen_range(0..all_hops.len()); - let hop = all_hops[rand_index]; let mut sup = self.spin_up; let mut sdown = self.spin_down; - match hop.2 { + match spin { Spin::Up => { // TODO: Combine these sets as a single bitmask. Maybe requires to modify the Trait BitOps - sup.set(hop.0); - sup.set(hop.1); + sup.set(rand_from_idx); + sup.set(rand_to_idx); }, Spin::Down => { - sdown.set(hop.0); - sdown.set(hop.1); + sdown.set(rand_from_idx); + sdown.set(rand_to_idx); } } - out_idx.0 = hop.0; - out_idx.1 = hop.1; - out_idx.2 = hop.2; + out_idx.0 = rand_from_idx; + out_idx.1 = rand_to_idx; + out_idx.2 = spin; FockState{ n_sites: max_size as usize, @@ -649,7 +682,6 @@ mod test { use super::*; use rand::SeedableRng; use rand::rngs::SmallRng; - type BitSize = u16; #[test] fn test_random_state_generator() { diff --git a/src/jastrow.rs b/src/jastrow.rs index 4e19b2c..87d966e 100644 --- a/src/jastrow.rs +++ b/src/jastrow.rs @@ -48,9 +48,17 @@ where let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); let k = indices[nk]; if n1.check(i) ^ n2.check(k) { - jastrow_out -= jastrow_params[i + k * n_sites]; + if i > k { + jastrow_out -= jastrow_params[i*(i-1)/2 + k]; + } else { + jastrow_out -= jastrow_params[k*(k-1)/2 + i]; + } } else { - jastrow_out += jastrow_params[i + k * n_sites]; + if i > k { + jastrow_out += jastrow_params[i*(i-1)/2 + k]; + } else { + jastrow_out += jastrow_params[k*(k-1)/2 + i]; + } } } i = regular_nor.leading_zeros() as usize; @@ -78,9 +86,17 @@ where let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); let k = indices[nk]; if n1.check(i) ^ n2.check(k) { - der.o_tilde[der.jas_off + i + k * n_sites + (der.n * der.mu) as usize] = -0.5; + if k > i { + der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = -0.5; + } else { + der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = -0.5; + } } else { - der.o_tilde[der.jas_off + i + k * n_sites + (der.n * der.mu) as usize] = 0.5; + if k > i { + der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = 0.5; + } else { + der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = 0.5; + } } } i = regular_nor.leading_zeros() as usize; @@ -110,9 +126,17 @@ fn jastrow_undo_update( } let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); if n1.check(i) ^ n2.check(index_j) { - *previous_jas += jastrow_params[i + index_j * n_sites]; + if i > index_j { + *previous_jas += jastrow_params[i*(i-1)/2 + index_j]; + } else { + *previous_jas += jastrow_params[index_j*(index_j-1)/2 + i]; + } } else { - *previous_jas -= jastrow_params[i + index_j * n_sites]; + if i > index_j { + *previous_jas -= jastrow_params[i*(i-1)/2 + index_j]; + } else { + *previous_jas -= jastrow_params[index_j*(index_j-1)/2 + i]; + } } i = spin_mask.leading_zeros() as usize; } @@ -141,9 +165,17 @@ fn jastrow_do_update( } let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); if n1.check(i) ^ n2.check(index_j) { - *previous_jas -= jastrow_params[i + index_j * n_sites]; + if i > index_j { + *previous_jas -= jastrow_params[i*(i-1)/2 + index_j]; + } else { + *previous_jas -= jastrow_params[index_j*(index_j-1)/2 + i]; + } } else { - *previous_jas += jastrow_params[i + index_j * n_sites]; + if i > index_j { + *previous_jas += jastrow_params[i*(i-1)/2 + index_j]; + } else { + *previous_jas += jastrow_params[index_j*(index_j-1)/2 + i]; + } } i = spin_mask.leading_zeros() as usize; } @@ -159,23 +191,38 @@ fn jastrow_single_update( index_j: usize, index_i: usize, sign: bool, - n_sites: usize, ) where T: BitOps, { if spin_mask.check(index_j) & spin_mask.check(index_i) { let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); - if n1.check(index_i) ^ n2.check(index_j) { - if sign { - *previous_jas -= jastrow_params[index_i + index_j * n_sites]; + if index_j > index_i { + if n1.check(index_i) ^ n2.check(index_j) { + if sign { + *previous_jas -= jastrow_params[index_j*(index_j-1)/2 + index_i]; + } else { + *previous_jas += jastrow_params[index_j*(index_j-1)/2 + index_i]; + } } else { - *previous_jas += jastrow_params[index_i + index_j * n_sites]; + if sign { + *previous_jas += jastrow_params[index_j*(index_j-1)/2 + index_i]; + } else { + *previous_jas -= jastrow_params[index_j*(index_j-1)/2 + index_i]; + } } } else { - if sign { - *previous_jas += jastrow_params[index_i + index_j * n_sites]; + if n1.check(index_i) ^ n2.check(index_j) { + if sign { + *previous_jas -= jastrow_params[index_i*(index_i-1)/2 + index_j]; + } else { + *previous_jas += jastrow_params[index_i*(index_i-1)/2 + index_j]; + } } else { - *previous_jas -= jastrow_params[index_i + index_j * n_sites]; + if sign { + *previous_jas += jastrow_params[index_i*(index_i-1)/2 + index_j]; + } else { + *previous_jas -= jastrow_params[index_i*(index_i-1)/2 + index_j]; + } } } } @@ -260,7 +307,6 @@ pub fn fast_update_jastrow( new_j, previous_j, false, - n_sites, ); // Now do @@ -292,7 +338,6 @@ pub fn fast_update_jastrow( new_j, previous_j, true, - n_sites, ); } @@ -467,6 +512,7 @@ mod test { const SIZE: usize = 8; let mut rng = SmallRng::seed_from_u64(42); const NSITES: usize = 8; + const NVIJ: usize = NSITES * (NSITES - 1) / 2; for _ in 0..100 { let up = rng.gen::(); let down = rng.gen::(); @@ -480,19 +526,22 @@ mod test { spin_down: down, n_sites: SIZE, }; - let mut jastrow_params: Vec = Vec::with_capacity(SIZE * SIZE); + let mut jastrow_params: Vec = Vec::with_capacity(NVIJ); + let mut jastrow_params_all: Vec = vec![0.0; NSITES * NSITES]; for _ in 0..SIZE * SIZE { jastrow_params.push(rng.gen::()); } - for i in 0..NSITES { - for j in 0..NSITES { - jastrow_params[j + i * NSITES] = jastrow_params[i + j * NSITES]; - if i == j {jastrow_params[i + j* NSITES] = 0.0;} + for j in 0..NSITES { + for i in 0..j { + jastrow_params_all[i + j * NSITES] = jastrow_params[(j - 1) * j / 2 + i]; + jastrow_params_all[j + i * NSITES] = jastrow_params[(j - 1) * j / 2 + i]; } } + println!("vij = {:?}", jastrow_params_all); + println!("closely packed vij = {:?}", jastrow_params); close( compute_jastrow_exp(fock_state1, &jastrow_params, NSITES), - compute_jastrow_easy_to_follow(fock_state2, &jastrow_params, NSITES), + compute_jastrow_easy_to_follow(fock_state2, &jastrow_params_all, NSITES), 1e-12, ) } @@ -501,6 +550,7 @@ mod test { #[test] fn test_fast_update_jastrow_u8() { const SIZE: usize = 8; + const NVIJ: usize = SIZE * (SIZE - 1) / 2; let mut rng = SmallRng::seed_from_u64(42); let mut up = rng.gen::(); let down = rng.gen::(); @@ -509,16 +559,11 @@ mod test { spin_down: down, n_sites: SIZE, }; - let mut jastrow_params: Vec = Vec::with_capacity(64); - for _ in 0..64 { + let mut jastrow_params: Vec = Vec::with_capacity(NVIJ); + for _ in 0..NVIJ { jastrow_params.push(rng.gen()); } - for i in 0..8 { - for j in 0..8 { - jastrow_params[j + i * 8] = jastrow_params[i + j * 8]; - } - } - let mut jastrow = compute_jastrow_exp(fock_state.clone(), &jastrow_params, 8); + let mut jastrow = compute_jastrow_exp(fock_state.clone(), &jastrow_params, SIZE); println!("previous_jas: {}", jastrow); for _ in 0..100 { let spin_update: f64 = rng.gen(); @@ -564,6 +609,7 @@ mod test { fn test_jastrow_u8_5sites() { let mut rng = SmallRng::seed_from_u64(42); const NSITES: usize = 5; + const NVIJ: usize = NSITES * (NSITES - 1) / 2; for _ in 0..100 { // Generate random state let mut fock_state1 = rng.gen::>(); @@ -573,20 +619,24 @@ mod test { // Copy the state let fock_state2 = fock_state1.clone(); - let mut jastrow_params: Vec = Vec::with_capacity(NSITES * NSITES); - for _ in 0..(NSITES * NSITES) { + let mut jastrow_params: Vec = Vec::with_capacity(NVIJ); + let mut jastrow_params_all: Vec = vec![0.0; NSITES * NSITES]; + for _ in 0..NVIJ { jastrow_params.push(rng.gen::()); } for i in 0..NSITES { - for j in 0..NSITES { - jastrow_params[j + i * NSITES] = jastrow_params[i + j * NSITES]; - if i == j {jastrow_params[i + j* NSITES] = 0.0;} + for j in i+1..NSITES { + if i == j { continue;} + jastrow_params_all[i + j * NSITES] = jastrow_params[(j - 1) * j / 2 + i]; + jastrow_params_all[j + i * NSITES] = jastrow_params[(j - 1) * j / 2 + i]; } } - println!("Params: {:?}", jastrow_params); + println!("vij = {:?}", jastrow_params_all); + println!("nvij = {}", jastrow_params_all.len()); + println!("closely packed vij = {:?}", jastrow_params); close( compute_jastrow_exp(fock_state1, &jastrow_params, NSITES), - compute_jastrow_easy_to_follow(fock_state2, &jastrow_params, NSITES), + compute_jastrow_easy_to_follow(fock_state2, &jastrow_params_all, NSITES), 1e-12, ) } diff --git a/src/lib.rs b/src/lib.rs index 711b0cc..949c4dc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,6 +1,8 @@ #[cfg(feature = "python-interface")] use pyo3::prelude::*; +use blas::{dcopy, dgemm}; + /// Input file parsing util. /// # Subfiles /// * __`orbitale.csv`__ - Variationnal parameters for the orbital. In csv @@ -78,10 +80,10 @@ impl<'a> std::fmt::Display for DerivativeOperator<'a> { let mut expval = " = ".to_owned(); let mut ho = " = ".to_owned(); let mut o_tilde = "O = ".to_owned(); - for mu in 0..(self.mu + 1) as usize { - expval.push_str(&format!("{:>width$.04e} ", self.expval_o[mu])); - ho.push_str(&format!("{:>width$.04e} ", self.ho[mu])); - for n in 0..self.n as usize { + for n in 0..self.n as usize { + expval.push_str(&format!("{:>width$.04e} ", self.expval_o[n])); + ho.push_str(&format!("{:>width$.04e} ", self.ho[n])); + for mu in 0..(self.mu + 1) as usize { o_tilde.push_str(&format!("{:>width$.04e}", self.o_tilde[n + n * mu])); } o_tilde.push_str("\n"); @@ -98,6 +100,118 @@ impl<'a> std::fmt::Display for DerivativeOperator<'a> { } } + +fn _save_otilde(der: &DerivativeOperator) { + let width = 16; + let mut c = vec![0.0; (der.n * der.n) as usize]; + println!("dim = {}", der.n * der.n); + unsafe { + dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, der.o_tilde, der.n, der.o_tilde, der.n, 0.0, &mut c, der.n); + } + let mut outstr = "".to_owned(); + outstr.push_str(&format!(" = ")); + for i in 0..der.n as usize { + outstr.push_str(&format!("\n ")); + for j in 0..der.n as usize { + outstr.push_str(&format!("{:>width$.04e}", c[i + der.n as usize * j])); + } + } + println!("{}", outstr); +} + +pub fn mapto_pairwf(input: &DerivativeOperator, output: &mut DerivativeOperator, sys: &SysParams) { + if input.mu != output.mu { + error!("Input dimension does not match output dimension. Make sure to match mu for both structures."); + panic!("Undefined Behavior."); + } + if input.mu < 0 { + error!("Dimension cannot be negative; is it empty?"); + panic!("Undefined Behavior"); + } + let nfij = sys.size*sys.size; + // Copy and scale fij from FIJ + for i in sys.ngi+sys.nvij+nfij..sys.ngi+sys.nvij+2*nfij { + unsafe { + //println!("{}", output.mu); + dcopy( + input.mu, + &input.o_tilde[i..input.n as usize + input.mu as usize * (i + 1)], + input.n, + &mut output.o_tilde[(i - nfij)..output.n as usize + output.mu as usize * (i - nfij + 1)], + output.n + ); + //for j in 0..nfij { + // dscal( + // output.mu, + // der_facteur, + // &mut output.o_tilde[j+sys.ngi+sys.nvij..output.n as usize + output.mu as usize * (j + 1 + sys.ngi + sys.nvij)], + // output.n, + // ) + //} + } + + } + // Copy Gutzwiller and jastrow + for i in 0..sys.ngi+sys.nvij { + unsafe { + dcopy( + input.mu, + &input.o_tilde[i..input.n as usize * (i + 1)], + input.n, + &mut output.o_tilde[i..output.n as usize * (i + 1)], + output.n + ); + } + + } + unsafe { + dcopy( + nfij as i32, + &input.expval_o[input.pfaff_off + nfij..input.pfaff_off + 2*nfij], + 1, + &mut output.expval_o[input.pfaff_off..output.pfaff_off + nfij], + 1 + ); + dcopy( + nfij as i32, + &input.ho[input.pfaff_off + nfij..input.pfaff_off + 2*nfij], + 1, + &mut output.ho[input.pfaff_off..input.pfaff_off + nfij], + 1 + ); + //dscal(nfij as i32, der_facteur, &mut output.expval_o[output.pfaff_off..output.pfaff_off + nfij], 1); + //dscal(nfij as i32, der_facteur, &mut output.ho[output.pfaff_off..output.pfaff_off + nfij], 1); + dcopy( + sys.nvij as i32, + &input.expval_o[input.jas_off..input.jas_off + sys.nvij], + 1, + &mut output.expval_o[output.jas_off..output.jas_off + sys.nvij], + 1 + ); + dcopy( + sys.nvij as i32, + &input.ho[input.jas_off..input.jas_off + sys.nvij], + 1, + &mut output.ho[output.jas_off..output.jas_off + sys.nvij], + 1 + ); + dcopy( + sys.ngi as i32, + &input.expval_o[0..sys.ngi], + 1, + &mut output.expval_o[0..sys.ngi], + 1 + ); + dcopy( + sys.ngi as i32, + &input.ho[0..sys.ngi], + 1, + &mut output.ho[0..sys.ngi], + 1 + ); + } +} + /// Module to calculate pfaffian /// # Usage /// Call the workspace query function on the $F^{\sigma\sigma'}_{ij}$ matrix and compute the diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs index 28972d7..c2b1260 100644 --- a/src/monte_carlo.rs +++ b/src/monte_carlo.rs @@ -2,8 +2,6 @@ use blas::daxpy; use log::{error, info, trace, warn}; use rand::distributions::{Distribution, Standard}; use rand::Rng; -use std::fs::File; -use std::io::Write; use crate::gutzwiller::compute_gutzwiller_der; use crate::jastrow::compute_jastrow_der; @@ -33,6 +31,7 @@ use crate::hamiltonian::{kinetic, potential}; // (ratio_ip, state2, updated_column, col_idx) //} +#[inline(always)] fn propose_hopping > @@ -47,13 +46,14 @@ T: BitOps + std::fmt::Display + std::fmt::Debug + From> ) -> (f64, FockState, Vec, usize) where Standard: Distribution { - let state2 = state.generate_hopping(rng, sys.size as u32, hop, sys); + let state2 = state.generate_hopping(rng, sys.size as u32, hop); let (ratio_ip, updated_column, col_idx) = { fast_internal_product(state, &state2, pfaff_state, &hop, previous_proj, params) }; (ratio_ip, state2, updated_column, col_idx) } +#[inline(always)] fn compute_hamiltonian(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams, sys: &SysParams) -> f64 { let kin = kinetic(state, pstate, proj, params, sys); let e = kin + potential(state, proj, pstate, sys); @@ -61,6 +61,7 @@ fn compute_hamiltonian(state: F e / (pstate.pfaff * ::exp(proj)) } +#[inline(always)] fn compute_derivative_operator> (state: FockState, pstate: &PfaffianState, der: &mut DerivativeOperator, sys: &SysParams) { @@ -95,7 +96,7 @@ fn make_update>( *ratio_prod *= ratio; } else { let tmp_pfaff = pstate.pfaff; - (*pstate, *proj) = compute_internal_product_parts(*state2, params); + (*pstate, *proj) = compute_internal_product_parts(*state2, params, sys); let inverse_proj = ::exp(*proj_copy_persistent - *proj); let err = ::abs(::abs(tmp_pfaff * *ratio * *ratio_prod * inverse_proj) - ::abs(pstate.pfaff)); if pstate.pfaff*pstate.pfaff < sys.tolerance_singularity { @@ -112,19 +113,187 @@ fn make_update>( } } +#[inline(always)] +fn warmup( + rng: &mut R, + state: &mut FockState, + hop: &mut (usize, usize, Spin), + proj: &mut f64, + proj_copy_persistent: &mut f64, + ratio_prod: &mut f64, + pstate: &mut PfaffianState, + n_accepted_updates: &mut usize, + params: &VarParams, + sys: &SysParams +) +where T: BitOps + From + std::fmt::Debug + std::fmt::Display, + R: Rng + ?Sized, + Standard: Distribution +{ + // Warmup + for _ in 0..sys.nmcwarmup { + let mut proj_copy = *proj; + let (mut ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, hop, rng, params, sys); + trace!("Current state: {}", state); + trace!("Proposed state: {}", state2); + trace!("Ratio: {}", ratio); + ratio *= ::exp(proj_copy - *proj); + let w = rng.gen::(); + if ::abs(ratio) * ::abs(ratio) >= w { + // We ACCEPT + trace!("Accept."); + *n_accepted_updates += 1; + make_update( + n_accepted_updates, + proj, + &mut proj_copy, + proj_copy_persistent, + &ratio, + ratio_prod, + state, + &state2, + &hop, + col, + colidx, + pstate, + params, + sys + ); + + } + } + +} + +#[inline(always)] +fn energy_error_estimation( + state_energy: f64, + previous_energies: &mut [f64], + energy_sums: &mut [f64], + energy_quad_sums: &mut [f64], + n_values: &mut [usize], + mc_it: usize, + error_estimation_level: usize, +) { + // Energy error estimation + let accumulation_level = ::log2((mc_it + 1) as f32) as usize; + for i in 0..accumulation_level { + if i >= error_estimation_level { break;} + if i == accumulation_level - 1{ + if i == 0 {break;} + previous_energies[i] = energy_sums[i - 1]; + } else { + n_values[i] += 1; + energy_sums[i] += 0.5 * (previous_energies[i] + state_energy); + energy_quad_sums[i] += 0.5 * (previous_energies[i]*previous_energies[i] + state_energy*state_energy); + } + } +} + +#[inline(always)] +fn accumulate_expvals(energy: &mut f64, state_energy: f64, der: &mut DerivativeOperator, rho: f64) { + // Accumulate Energy + *energy += state_energy; + let last_line_begin = (der.n * der.mu) as usize; + let last_line_end = (der.n * (der.mu + 1)) as usize; + // Accumulate + unsafe { + let incx = 1; + let incy = 1; + daxpy( + der.n, + state_energy, + &der.o_tilde[last_line_begin .. last_line_end], + incx, + &mut der.ho, + incy + ); + } + // Accumulate + for i in 0 .. der.n as usize { + der.expval_o[i] += der.o_tilde[i + last_line_begin] * rho; + } +} + +#[inline(always)] +fn normalize(energy: &mut f64, energy_bootstraped: &mut f64, expval_o: &mut [f64], ho: &mut [f64], nsample: f64, nbootstrap: f64, nparams: usize) { + *energy *= 1.0 / nsample; + *energy_bootstraped *= 1.0 / nsample; + *energy_bootstraped *= 1.0 / nbootstrap; + for i in 0..nparams { + expval_o[i] *= 1.0 / nsample; + ho[i] *= 1.0 / nsample; + } +} + +#[inline(always)] +fn sq(x: f64) -> f64 +{ + x * x +} + +pub fn compute_mean_energy_exact(params: &VarParams, sys: &SysParams, der: &mut DerivativeOperator) -> f64 +{ + if der.mu != -1 { + error!("The derivative operator current row was mu = {} on entry, is it reinitialized?", der.mu); + } + //let all_states = initial_state.generate_all_hoppings(sys.hopping_bitmask); + let all_states = vec![ + FockState{spin_up: 128u8, spin_down: 128u8, n_sites: 2}, + FockState{spin_up: 64u8, spin_down: 128u8, n_sites: 2}, + FockState{spin_up: 64u8, spin_down: 64u8, n_sites: 2}, + FockState{spin_up: 128u8, spin_down: 64u8, n_sites: 2}, + ]; + + let mut energy = 0.0; + let mut norm = 0.0; + der.mu = 0; + let mut state_energy; + + for state2 in all_states.iter() { + let (pstate, proj) = compute_internal_product_parts(*state2, params, sys); + norm += sq(::abs(pstate.pfaff) * ::exp(proj)); + let rho = ::abs(pstate.pfaff) * ::exp(proj); + println!("rho = {}", sq(rho)); + compute_derivative_operator(*state2, &pstate, der, sys); + state_energy = compute_hamiltonian(*state2, &pstate, proj, params, sys) * sq(rho); + accumulate_expvals(&mut energy, state_energy, der, sq(rho)); + der.visited[der.mu as usize] = 1; + let mut outstr = "".to_owned(); + outstr.push_str(&format!("H_mu O_k mu =")); + for i in 0..der.n as usize { + outstr.push_str(&format!("{} ", der.o_tilde[i + (der.n * der.mu) as usize] * state_energy)); + der.o_tilde[i + (der.n * der.mu) as usize] *= rho; + } + println!("{}", outstr); + der.mu += 1; + + } + println!("energy = {}", energy); + for i in 0.. der.n as usize { + der.expval_o[i] *= 1.0 / norm; + der.ho[i] *= 1.0 / norm; + } + for i in 0..der.mu as usize { + for j in 0..der.n as usize { + der.o_tilde[j + i * der.n as usize] *= 1.0 / ::sqrt(norm); + } + } + energy / norm +} + pub fn compute_mean_energy > -(rng: &mut R, initial_state: FockState, params: &VarParams, sys: &SysParams, derivatives: &mut DerivativeOperator) -> (f64, Vec>, f64) +(rng: &mut R, initial_state: FockState, params: &VarParams, sys: &SysParams, derivatives: &mut DerivativeOperator) -> (f64, Vec>, f64, f64) where Standard: Distribution { - let mut ratiofp = File::create("ratios").unwrap(); if derivatives.mu != -1 { warn!("The derivative operator current row was mu = {} on entry, is it reinitialized?", derivatives.mu); } let mut state = initial_state; let mut accumulated_states: Vec> = Vec::new(); - let (mut pstate, mut proj) = compute_internal_product_parts(state, params); + let (mut pstate, mut proj) = compute_internal_product_parts(state, params, sys); let mut hop: (usize, usize, Spin) = (0, 0, Spin::Up); let mut _lip = ::ln(::abs(::exp(proj) * pstate.pfaff)) * 2.0; let mut n_accepted_updates: usize = 0; @@ -133,49 +302,20 @@ where Standard: Distribution let mut proj_copy_persistent = proj; let mut ratio_prod = 1.0; let mut sample_counter: usize = 0; - if ::log2(sys.nmcsample as f64) <= 5.0 { + if ::log2(sys.nmcsample as f64) <= 6.0 { error!("Not enough monte-carlo sample for an accurate error estimation. NMCSAMPLE = {}", sys.nmcsample); - panic!("NMCSAMPLE is less than 32."); + panic!("NMCSAMPLE is less than 64."); } let error_estimation_level = ::log2(sys.nmcsample as f64) as usize - 5; let mut energy_sums = vec![0.0; error_estimation_level]; let mut energy_quad_sums = vec![0.0; error_estimation_level]; let mut previous_energies = vec![0.0; error_estimation_level + 1]; let mut n_values = vec![0; error_estimation_level]; + let mut energy_bootstraped = 0.0; info!("Starting the warmup phase."); - // Warmup - for _ in 0..sys.nmcwarmup { - let mut proj_copy = proj; - let (mut ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, &mut hop, rng, params, sys); - trace!("Current state: {}", state); - trace!("Proposed state: {}", state2); - trace!("Ratio: {}", ratio); - ratio *= ::exp(proj_copy - proj); - let w = rng.gen::(); - if ::abs(ratio) * ::abs(ratio) >= w { - // We ACCEPT - trace!("Accept."); - n_accepted_updates += 1; - make_update( - &mut n_accepted_updates, - &mut proj, - &mut proj_copy, - &mut proj_copy_persistent, - &ratio, - &mut ratio_prod, - &mut state, - &state2, - &hop, - col, - colidx, - &mut pstate, - params, - sys - ); - - } - } + warmup(rng, &mut state, &mut hop, &mut proj, &mut proj_copy_persistent, &mut ratio_prod, &mut + pstate, &mut n_accepted_updates, params, sys); info!("Starting the sampling phase."); // MC Sampling @@ -193,7 +333,6 @@ where Standard: Distribution trace!("Ratio: {}", ratio); ratio *= ::exp(proj_copy - proj); let w = rng.gen::(); - ratiofp.write(format!("De {}, à {} = {}\n", state, state2, ::abs(ratio) * ::abs(ratio)).as_bytes()).unwrap(); if ::abs(ratio) * ::abs(ratio) >= w { // We ACCEPT trace!("Accept."); @@ -222,59 +361,36 @@ where Standard: Distribution } if sample_counter >= sys.mcsample_interval { accumulated_states.push(state); - // Accumulate energy derivatives.visited[derivatives.mu as usize] += 1; let state_energy = compute_hamiltonian(state, &pstate, proj, params, sys); - energy += state_energy; - _energy_sq += state_energy*state_energy; - // Energy error estimation - let accumulation_level = ::log2((mc_it + 1) as f32) as usize; - for i in 0..accumulation_level { - if i >= error_estimation_level { break;} - if i == accumulation_level - 1{ - if i == 0 {break;} - previous_energies[i] = energy_sums[i - 1]; - } else { - n_values[i] += 1; - energy_sums[i] += 0.5 * (previous_energies[i] + state_energy); - energy_quad_sums[i] += 0.5 * (previous_energies[i]*previous_energies[i] + state_energy*state_energy); - } - } - // Accumulate - unsafe { - let incx = 1; - let incy = 1; - daxpy( - derivatives.n, - state_energy, - &derivatives.o_tilde[(derivatives.n * derivatives.mu) as usize .. (derivatives.n * (derivatives.mu + 1)) as usize], - incx, - &mut derivatives.ho, - incy - ); - } - // Accumulate the derivative operators - for i in 0 .. derivatives.n as usize { - derivatives.expval_o[i] += derivatives.o_tilde[i + (derivatives.n * derivatives.mu) as usize]; - } + accumulate_expvals(&mut energy, state_energy, derivatives, 1.0); + energy_error_estimation(state_energy, &mut previous_energies, &mut energy_sums, &mut + energy_quad_sums, &mut n_values, mc_it, error_estimation_level); sample_counter = 0; + if mc_it >= (sys.nmcsample - sys.nbootstrap) * sys.mcsample_interval { + energy_bootstraped += energy; + } } sample_counter += 1; } derivatives.mu += 1; info!("Final Energy: {:.2}", energy); - energy = energy / sys.nmcsample as f64; + normalize(&mut energy, &mut energy_bootstraped, &mut derivatives.expval_o, &mut derivatives.ho, + sys.nmcsample as f64, sys.nbootstrap as f64, derivatives.n as usize); info!("Final Energy normalized: {:.2}", energy); // Error estimation let mut error = vec![0.0; error_estimation_level]; for i in 0..error_estimation_level { error[i] = ::sqrt( - (energy_quad_sums[i] - energy_sums[i]*energy_sums[i] / (n_values[i] as f64)) / + (energy_quad_sums[i] - energy_sums[i]*energy_sums[i] / ((n_values[i]*n_values[i]) as f64)) / (n_values[i] * (n_values[i] - 1)) as f64 - ) + ); + } + if derivatives.mu == -1 { + warn!("Parameter mu was -1 on exit, was it updated?"); } let correlation_time = 0.5 * ((error[error_estimation_level-1]/error[0])*(error[error_estimation_level-1]/error[0]) - 1.0); - (energy, accumulated_states, correlation_time) + (energy_bootstraped, accumulated_states, error[error_estimation_level - 1], correlation_time) } diff --git a/src/optimisation.rs b/src/optimisation.rs index be91d84..ffc773c 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -1,15 +1,17 @@ -use blas::{daxpy, dcopy, ddot, dgemv, dscal}; -use log::trace; +use blas::{daxpy, dcopy, ddot, dgemm, dgemv, dger, dscal}; +use lapack::dsyev; +use log::{error, trace}; +use colored::Colorize; use crate::DerivativeOperator; -fn gradient(x: &[f64], a: &DerivativeOperator, b: &mut [f64], dim: i32) { +fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], dim: i32, mu: i32, nsamp: f64) { let alpha = -1.0; let incx = 1; let incy = 1; let mut work: Vec = vec![0.0; dim as usize]; // Compute Ax - compute_w(&mut work, a, x); + compute_w(&mut work, otilde, visited, expval_o, x, dim, mu, nsamp); unsafe { // Compute b - Ax daxpy(dim, alpha, &work, incx, b, incy); @@ -24,7 +26,7 @@ fn update_x(x: &mut [f64], pk: &[f64], alpha: f64, dim: i32) { }; } -fn compute_w(w: &mut [f64], a: &DerivativeOperator, p: &[f64]) { +fn compute_w(w: &mut [f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], p: &[f64], dim: i32, mu: i32, nsamp: f64) { // Computes Ap let incx = 1; let incy = 1; @@ -32,22 +34,24 @@ fn compute_w(w: &mut [f64], a: &DerivativeOperator, p: &[f64]) { let beta = 0.0; // Temp work vector - let mut work: Vec = vec![0.0; a.mu as usize]; + let mut work: Vec = vec![0.0; mu as usize]; unsafe { trace!("x_[n] = {:?}", p); - trace!("mu = {}, n = {}", a.mu, a.n); - dgemv(b"T"[0], a.n, a.mu, alpha, a.o_tilde, a.n, p, incx, beta, &mut work, incy); - for i in 0..a.mu as usize { - work[i] *= a.visited[i] as f64; + trace!("mu = {}, n = {}", mu, dim); + // 80 misawa + dgemv(b"T"[0], dim, mu, alpha, otilde, dim, p, incx, beta, &mut work, incy); + for i in 0..mu as usize { + work[i] *= visited[i] as f64; } trace!("O^[T]_[mu, n] x_[n] = {:?}", work); - trace!("Len(work) = {}, a.n = {}, a.mu = {}", work.len(), a.n, a.mu); - trace!("~O_[0, mu] = {:?}", a.o_tilde.iter().step_by(a.mu as usize).copied().collect::>()); + trace!("Len(work) = {}, a.n = {}, a.mu = {}", work.len(), dim, mu); + trace!("~O_[0, mu] = {:?}", otilde.iter().step_by(mu as usize).copied().collect::>()); // Sometimes segfaults - dgemv(b"N"[0], a.n, a.mu, 1.0 / a.nsamp, a.o_tilde, a.n, &work, incx, beta, w, incy); + dgemv(b"N"[0], dim, mu, 1.0 / nsamp, otilde, dim, &work, incx, beta, w, incy); trace!("O_[m, mu] O^[T]_[mu, n] x_[n] = {:?}", w); - let alpha = ddot(a.n, a.expval_o, incx, p, incy); - daxpy(a.n, - alpha, a.expval_o, incx, w, incy); + let alpha = ddot(dim, expval_o, incx, p, incy); + // 81 misawa + daxpy(dim, - alpha, expval_o, incx, w, incy); } } @@ -111,44 +115,340 @@ pub fn spread_eigenvalues(a: &mut DerivativeOperator) { } } +fn prefilter_overlap_matrix(a: &DerivativeOperator, _ignore_idx: &mut [bool], dim: i32, _diag_threshold: f64) -> usize { + // Loop over diagonal elements of S_km + // Reminder, S_{kk} = 1/N_{\rm MC.} \sum_\mu \tilde{O}^*_{k\mu}\tilde{O}^T_{\mu k} - + // \Re{\expval{O_k}}^2 + + let skip_param_count: usize = 0; + let mut diag_elem = vec![0.0; dim as usize]; + for k in 0..dim as usize { + // Start with 1/N_{\rm MC.} \sum_\mu \tilde{O}^*_{k\mu}\tilde{O}^T_{\mu k} + let z1: f64 = unsafe { + ddot( + a.mu, + &a.o_tilde[k..a.n as usize + a.mu as usize * (k + 1)], + a.n, + &a.o_tilde[k..a.n as usize + a.mu as usize * (k + 1)], + a.n + ) + }; + + // Now \Re{\expval{O_k}}^2 + let z2: f64 = a.expval_o[k] * a.expval_o[k]; + + diag_elem[k] = z1 - z2; + } + + let mut max_elem = ::MIN; + for k in 0..dim as usize { + if diag_elem[k] > max_elem { + max_elem = diag_elem[k]; + } + } + //for k in 0..dim as usize { + // if diag_elem[k] < threshold { + // skip_param_count += 1; + // ignore_idx[k] = true; + // } + //} + dim as usize - skip_param_count + +} + +fn cpy_segmented_matrix_to_dense(a: &DerivativeOperator, output_otilde: &mut [f64], output_expvalo: &mut [f64], ignore_idx: &[bool], dim: i32, nparams_opt: usize) { + let mut j: usize = 0; + for k in 0..dim as usize { + if ignore_idx[k] { + continue; + } + unsafe { + dcopy( + a.mu, + &a.o_tilde[k..a.n as usize + a.mu as usize * (k + 1)], + a.n, + &mut output_otilde[j..a.mu as usize * (j+1)], + nparams_opt as i32 + ); + } + output_expvalo[j] = a.expval_o[k]; + j += 1; + + } +} + +fn compute_s_explicit(otilde: &[f64], expval_o: &[f64], visited: &[usize], dim: i32, mu: i32, nsamp: f64, epsilon: f64) -> Vec { + // Computes Ap + let incx = 1; + let alpha = 1.0/nsamp; + let beta = -(1.0); + let mut otilde_w_visited = vec![0.0; (dim * mu) as usize]; + for i in 0..mu as usize { + for j in 0..dim as usize { + otilde_w_visited[j + i*dim as usize] = otilde[j + i*dim as usize] * visited[i] as f64; + } + } + + // Temp work vector + let mut work = vec![0.0; (dim * dim) as usize]; + unsafe { + // 80 misawa + dger(dim, dim, 1.0, &expval_o, incx, &expval_o, incx, &mut work, dim); + dgemm(b"N"[0], b"T"[0], dim, dim, mu, alpha, õ_w_visited, dim, õ_w_visited, dim, beta, &mut work, dim); + } + // Shift smallest eigenvalues + for i in 0..dim as usize { + work[i + dim as usize * i] += epsilon; + } + work +} + +fn _save_otilde(der: &[f64], mu: usize, n: usize) -> String { + let width = 16; + let mut o_tilde = "".to_owned(); + for m in 0..mu { + for i in 0..n { + if i == m { + o_tilde.push_str(&format!("{:>width$.04e}", der[i + m * n]).yellow()); + } else { + o_tilde.push_str(&format!("{:>width$.04e}", der[i + m * n])); + } + } + o_tilde.push_str("\n"); + } + o_tilde +} + +fn diagonalize_dense_matrix(s: &mut [f64], dim: i32) -> Vec { + let jobz = b"V"[0]; + let uplo = b"U"[0]; + let mut w = vec![0.0; dim as usize]; + let lwork = 3*(dim); + let mut work = vec![0.0; lwork as usize]; + let mut info = 0; + unsafe { + dsyev( + jobz, + uplo, + dim, + s, + dim, + &mut w, + &mut work, + lwork, + &mut info + ); + } + if info < 0{ + error!("Parameter {} had an illegal value in call to lapack::dsyev.", ::abs(info)); + } + else if info > 0{ + error!("Convergence was not achieved to diagonalize the overlap matrix. + {} off-diagonal elements did not converge to 0.", info); + } + + w +} + +fn compute_delta_from_eigenvalues(x0: &mut [f64], eigenvectors: &[f64], eigenvalues: &[f64], dim: i32) { + let trans = b"T"[0]; + let incx = 1; + let incy = 1; + let mut work = vec![0.0; dim as usize]; + unsafe { + dgemv(trans, dim, dim, 1.0, eigenvectors, dim, x0, incx, 0.0, &mut work, incy); + } + for i in 0..dim as usize{ + work[i] *= eigenvalues[i]; + } + let trans = b"N"[0]; + unsafe { + dgemv(trans, dim, dim, 1.0, eigenvectors, dim, &work, incx, 0.0, x0, incy); + } +} + +fn _compute_matrix_product(s: &mut [f64], eigenvectors: &[f64], eigenvalues: &[f64], dim: i32) { + let transb = b"T"[0]; + let incx = 1; + let incy = 1; + let mut work = vec![0.0; (dim * dim) as usize]; + let mut s_copy = vec![0.0; (dim * dim) as usize]; + let mut work_direct = vec![0.0; (dim * dim) as usize]; + let mut s_copy_direct = vec![0.0; (dim * dim) as usize]; + unsafe { + dcopy(dim*dim, eigenvectors, incx, &mut work, incy); + dcopy(dim*dim, eigenvectors, incx, &mut work_direct, incy); + //dgemm(transa, transb, dim, dim, dim, 1.0, eigenvectors, dim, s, dim, 0.0, &mut s_copy, dim); + } + for i in 0..dim as usize{ + for j in 0..dim as usize { + if eigenvalues[i] < 1e-1 { + work[j + dim as usize * i] *= 0.0; + work_direct[j + dim as usize * i] *= 0.0; + continue; + } + work[j + dim as usize * i] *= 1.0 / eigenvalues[i]; + work_direct[j + dim as usize * i] *= eigenvalues[i]; + } + } + let transa = b"N"[0]; + unsafe { + dgemm(transa, transb, dim, dim, dim, 1.0, eigenvectors, dim, &work, dim, 0.0, &mut s_copy, dim); + dgemm(transa, transb, dim, dim, dim, 1.0, eigenvectors, dim, &work_direct, dim, 0.0, &mut s_copy_direct, dim); + //dcopy(dim*dim, s, incx, &mut work, incy); + let transa = b"N"[0]; + let transb = b"T"[0]; + dgemm(transa, transb, dim, dim, dim, 1.0, &s_copy, dim, &s_copy_direct, dim, 0.0, s, dim); + } +} + +pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64, dim: i32, thresh: f64) -> Vec{ + // PRE FILTER + let mut ignore = vec![false; dim as usize]; + //println!("{}", save_otilde(a.o_tilde, a.mu as usize, a.n as usize)); + let mut unfiltered_s = compute_s_explicit(a.o_tilde, a.expval_o, a.visited, dim, a.mu, a.nsamp, epsilon); + //println!("dim = {}, Unfiltered S = ", dim); + //println!("{}", save_otilde(&unfiltered_s, dim as usize, dim as usize)); + let new_dim = prefilter_overlap_matrix(a, &mut ignore, dim, thresh); + //println!("ignore : {:?}", ignore); + let mut otilde = vec![0.0; new_dim * a.mu as usize]; + let mut expvalo = vec![0.0; new_dim]; + cpy_segmented_matrix_to_dense(a, &mut otilde, &mut expvalo, &ignore, dim, new_dim); + //let filtered_s = compute_s_explicit(õ, &expvalo, a.visited, new_dim as i32, a.mu, a.nsamp); + let mut eigenvalues = diagonalize_dense_matrix(&mut unfiltered_s, dim); + //let eigenvectors = &unfiltered_s; + //println!("S = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); + //compute_matrix_product(&mut s_copy, &eigenvectors, &eigenvalues, dim); + //println!("S^[-1] = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); + //panic!("stop"); + //let mut x0_raw = vec![0.0; dim as usize]; + //unsafe { + // let incx = 1; + // let incy = 1; + // dgemv(b"T"[0], dim, dim, 1.0, &s_copy, dim, x0, incx, 0.0, &mut x0_raw, incy); + //} + //println!("UD^[-1]U^[T] = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); + //println!("x0 = {:?}", x0_raw); + //println!("eigenvalues: {:?}", eigenvalues); + + // Remove problematic eigenvalue + let mut max = ::MIN; + for e in eigenvalues.iter() { + if *e > max { + max = *e; + } + } + let threshold = thresh * max; + let mut i = 0; + for e in eigenvalues.iter_mut() { + if *e < threshold { + *e = 0.0; + ignore[i] = true; + } + i+=1; + } + //Invert matrix + for e in eigenvalues.iter_mut() { + if *e == 0.0 { + continue; + } + *e = 1.0 / *e; + } + compute_delta_from_eigenvalues(b, &unfiltered_s, &eigenvalues, dim); + return ignore; +} + /// Computes the solution of $A\mathbf{x}-\mathbf{b}=0$ /// TODOC -pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], epsilon: f64, kmax: usize, dim: i32) { +pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], _epsilon: f64, kmax: usize, dim: i32, thresh: f64, epsilon_convergence: f64) -> Vec +{ + // PRE FILTER + let mut ignore = vec![false; dim as usize]; + //println!("{}", save_otilde(a.o_tilde, a.mu as usize, a.n as usize)); + //println!("dim = {}, Unfiltered S = ", dim); + //println!("{}", save_otilde(&unfiltered_s, dim as usize, dim as usize)); + let new_dim = prefilter_overlap_matrix(a, &mut ignore, dim, thresh); + //println!("ignore : {:?}", ignore); + let mut otilde = vec![0.0; new_dim * a.mu as usize]; + let mut expvalo = vec![0.0; new_dim]; + cpy_segmented_matrix_to_dense(a, &mut otilde, &mut expvalo, &ignore, dim, new_dim); + //let filtered_s = compute_s_explicit(õ, &expvalo, a.visited, new_dim as i32, a.mu, a.nsamp); + //let eigenvectors = &unfiltered_s; + //println!("S = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); + //compute_matrix_product(&mut s_copy, &eigenvectors, &eigenvalues, dim); + //println!("S^[-1] = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); + //panic!("stop"); + //let mut x0_raw = vec![0.0; dim as usize]; + //unsafe { + // let incx = 1; + // let incy = 1; + // dgemv(b"T"[0], dim, dim, 1.0, &s_copy, dim, x0, incx, 0.0, &mut x0_raw, incy); + //} + //println!("UD^[-1]U^[T] = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); + //println!("x0 = {:?}", x0_raw); + //println!("eigenvalues: {:?}", eigenvalues); + + // Remove problematic eigenvalue + //fp.write_all(save_otilde(&filtered_s, new_dim as usize, new_dim as usize).as_bytes()).unwrap(); + //println!("dim = {}, Filtered S = ", new_dim); + //println!("{}", save_otilde(&filtered_s, new_dim as usize, new_dim as usize)); + //println!("{}", save_otilde(õ, a.mu as usize, new_dim as usize)); + trace!("Initial guess x_0: {:?}", x0); trace!("Initial equation b: {:?}", b); - let mut w: Vec = vec![0.0; dim as usize]; + let mut w: Vec = vec![0.0; new_dim]; // Error threshold - let e = unsafe { - let incx = 1; - let incy = 1; - ddot(dim, b, incx, b, incy) * epsilon - }; + let mut e = 0.0; + for i in 0..dim as usize { + if ignore[i] { + continue; + } + e += b[i] * b[i]; + } + e *= epsilon_convergence; trace!("Error threshold e = {}", e); - gradient(x0, a, b, dim); - let mut p = Vec::with_capacity(dim as usize); + //println!("Error threshold e = {}", e); + gradient(x0, õ, a.visited, &expvalo, b, new_dim as i32, a.mu, a.nsamp); + let mut p = vec![0.0; new_dim]; + let mut j: usize = 0; + for i in 0..dim as usize { + if ignore[i] { + continue; + } + p[j] = b[i]; + j += 1; + } unsafe { - let incx = 1; - let incy = 1; - p.set_len(dim as usize); - dcopy(dim, b, incx, &mut p, incy) + dcopy(new_dim as i32, &p, 1, b, 1); } let mut alpha = 0.0; for k in 0..kmax { trace!("r_{} : {:?}", k, b); trace!("p_{} : {:?}", k, p); - compute_w(&mut w, a, &p); - let nrm2rk = alpha_k(b, &p, &w, &mut alpha, dim); + //println!("r_{} : {:?}", k, b); + //println!("p_{} : {:?}", k, p); + compute_w(&mut w, õ, a.visited, &expvalo, &p, new_dim as i32, a.mu, a.nsamp); + //println!("w_{} : {:?}", k, w); + let nrm2rk = alpha_k(b, &p, &w, &mut alpha, new_dim as i32); trace!("alpha_{} : {}", k, alpha); - update_x(x0, &p, alpha, dim); + //println!("alpha_{} : {}", k, alpha); + //if alpha < 0.0 { + // //error!("Input overlap matrix S was not positive-definite."); + // break; + // //panic!("p^T S p < 0.0"); + //} + update_x(x0, &p, alpha, new_dim as i32); trace!("x_{} : {:?}", k+1, x0); - update_r(b, &w, alpha, dim); - let beta = beta_k(b, nrm2rk, dim); + //println!("x_{} : {:?}", k+1, x0); + update_r(b, &w, alpha, new_dim as i32); + let beta = beta_k(b, nrm2rk, new_dim as i32); if beta * nrm2rk < e { trace!("Achieved convergence at {} iterations", k); break; } trace!("beta_{} : {}", k, beta); - update_p(b, &mut p, beta, dim); + update_p(b, &mut p, beta, new_dim as i32); } + ignore } diff --git a/src/pfaffian.rs b/src/pfaffian.rs index 53a8b44..ed769a6 100644 --- a/src/pfaffian.rs +++ b/src/pfaffian.rs @@ -14,7 +14,7 @@ use std::fmt; /// * __`matrix`__ - The matrix $A$. This is the matrix that we need the pfaffian /// to get the inner product $\braket{x}{\phi_{\text{PF}}}$ /// * __`curr_state`__ - The current state of that the matrix $A$ is written in. -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct PfaffianState { pub n_elec: usize, pub n_sites: usize, @@ -60,7 +60,7 @@ impl fmt::Display for PfaffianState { /// * __`a`__ - The matrix $A$ of the variationnal parameters. /// * __`n`__ - The dimension of the matrix, this correspond to the number of /// electrons. -fn invert_matrix(a: &mut [f64], n: i32) { +fn invert_matrix(a: &mut [f64], n: i32) -> f64 { // Info output of lapack let mut info1: i32 = 0; let mut info2: i32 = 0; @@ -69,11 +69,18 @@ fn invert_matrix(a: &mut [f64], n: i32) { let n_entry: i32 = n * n; // Workspaces let mut work: Vec = Vec::with_capacity(n_entry as usize); - let mut ipiv: Vec = Vec::with_capacity(n as usize); + let mut ipiv: Vec = vec![0; n as usize]; // Inverse matrix `a` inplace using L*U decomposition. + let mut determinant = 1.0; unsafe { dgetrf(n, n, a, n, &mut ipiv, &mut info1); + for i in 0..n as usize { + determinant *= a[i + i*n as usize]; + if ipiv[i] != (i + 1) as i32 { + determinant *= -1.0; + } + } dgetri(n, a, n, &ipiv, &mut work, n_entry, &mut info2); } @@ -88,6 +95,7 @@ fn invert_matrix(a: &mut [f64], n: i32) { ); panic!("Matrix invertion fail."); } + determinant } #[allow(dead_code)] @@ -117,7 +125,7 @@ fn transpose(a: &Vec, n: usize) -> Vec{ /// # Fields /// * __`fij`__ - All the variationnal parameters. /// * __`state`__ - The state of the system. -pub fn construct_matrix_a_from_state(fij: &[f64], state: FockState) -> PfaffianState +pub fn construct_matrix_a_from_state(fij: &[f64], state: FockState, sys: &SysParams) -> PfaffianState where T: BitOps + std::fmt::Display, { @@ -161,30 +169,32 @@ where -fij[indices2[jj] + size * indices[ii] + size*size]; } } - for jj in 0..indices.len() { - for ii in 0..indices.len() { - if indices[ii] == indices[jj] { - continue; + if !sys.pair_wavefunction { + for jj in 0..indices.len() { + for ii in 0..indices.len() { + if indices[ii] == indices[jj] { + continue; + } + a[ii + jj * n] = + fij[indices[ii] + size * indices[jj]] + -fij[indices[jj] + size * indices[ii]]; + a[jj + ii * n] = + fij[indices[jj] + size * indices[ii]] + -fij[indices[ii] + size * indices[jj]]; } - a[ii + jj * n] = - fij[indices[ii] + size * indices[jj]] - -fij[indices[jj] + size * indices[ii]]; - a[jj + ii * n] = - fij[indices[jj] + size * indices[ii]] - -fij[indices[ii] + size * indices[jj]]; } - } - for jj in 0..indices2.len() { - for ii in 0..indices2.len() { - if indices2[ii] == indices2[jj] { - continue; + for jj in 0..indices2.len() { + for ii in 0..indices2.len() { + if indices2[ii] == indices2[jj] { + continue; + } + a[ii + off + (jj + off) * n] = + fij[indices2[ii] + size * indices2[jj] + 3*size*size] + -fij[indices2[jj] + size * indices2[ii] + 3*size*size]; + a[jj + off + (ii + off) * n] = + fij[indices2[jj] + size * indices2[ii] + 3*size*size] + -fij[indices2[ii] + size * indices2[jj] + 3*size*size]; } - a[ii + off + (jj + off) * n] = - fij[indices2[ii] + size * indices2[jj] + 3*size*size] - -fij[indices2[jj] + size * indices2[ii] + 3*size*size]; - a[jj + off + (ii + off) * n] = - fij[indices2[jj] + size * indices2[ii] + 3*size*size] - -fij[indices2[ii] + size * indices2[jj] + 3*size*size]; } } @@ -224,22 +234,24 @@ pub fn compute_pfaffian_derivative(pstate: &PfaffianState, der: &mut DerivativeO trace!("~O_[{}, {}] = {}", der.pfaff_off + indices[ii] + size * indices2[jj] + 2 * size * size, der.mu, -a[(jj + off) * n + ii]); } } - for jj in 0..indices.len() { - for ii in 0..indices.len() { - if indices[ii] == indices[jj] { - continue; + if !sys.pair_wavefunction { + for jj in 0..indices.len() { + for ii in 0..indices.len() { + if indices[ii] == indices[jj] { + continue; + } + der.o_tilde[der.pfaff_off + indices[ii] + size * indices[jj] + (der.n * der.mu) as usize] = -a[ii + jj * n]; + der.o_tilde[der.pfaff_off + indices[jj] + size * indices[ii] + (der.n * der.mu) as usize] = -a[jj + ii * n]; } - der.o_tilde[der.pfaff_off + indices[ii] + size * indices[jj] + (der.n * der.mu) as usize] = -a[ii + jj * n]; - der.o_tilde[der.pfaff_off + indices[jj] + size * indices[ii] + (der.n * der.mu) as usize] = -a[jj + ii * n]; } - } - for jj in 0..indices2.len() { - for ii in 0..indices2.len() { - if indices2[ii] == indices2[jj] { - continue; + for jj in 0..indices2.len() { + for ii in 0..indices2.len() { + if indices2[ii] == indices2[jj] { + continue; + } + der.o_tilde[der.pfaff_off + indices2[ii] + size * indices2[jj] + 3*size*size + (der.n * der.mu) as usize] = -a[ii + off + (jj + off) * n]; + der.o_tilde[der.pfaff_off + indices2[jj] + size * indices2[ii] + 3*size*size + (der.n * der.mu) as usize] = -a[jj + off + (ii + off) * n]; } - der.o_tilde[der.pfaff_off + indices2[ii] + size * indices2[jj] + 3*size*size + (der.n * der.mu) as usize] = -a[ii + off + (jj + off) * n]; - der.o_tilde[der.pfaff_off + indices2[jj] + size * indices2[ii] + 3*size*size + (der.n * der.mu) as usize] = -a[jj + off + (ii + off) * n]; } } } @@ -262,7 +274,6 @@ pub fn get_pfaffian_ratio( spin: Spin, fij: &[f64], ) -> (f64, Vec, usize) { - // Rename let indx_up = &previous_pstate.indices.0; trace!("Up : {:?}", indx_up); @@ -631,6 +642,26 @@ mod tests { #[test] fn test_pfaffian_update_random_no_sign_correction() { const SIZE: usize = 8; + let sys = crate::SysParams { + size: SIZE, + nelec: 0, + array_size: (SIZE + 7) / 8, + cons_t: -1.0, + cons_u: 1.0, + nfij: 4*SIZE*SIZE, + nvij: SIZE*(SIZE-1)/2, + ngi: SIZE, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &[], + hopping_bitmask: &[], + clean_update_frequency: 0, + nmcwarmup: 0, + nmcsample: 0, + tolerance_sherman_morrison: 0.0, + tolerance_singularity: 0.0, + pair_wavefunction: false, + }; let mut rng = SmallRng::seed_from_u64(42); // Size of the system let mut params = vec![0.0; 4 * SIZE * SIZE]; @@ -655,7 +686,7 @@ mod tests { // Matrix needs to be even sized if n % 2 == 1 { continue;} println!("------------- Initial State ----------------"); - let mut pfstate = construct_matrix_a_from_state(¶ms, state); + let mut pfstate = construct_matrix_a_from_state(¶ms, state, &sys); println!("Inverse Matrix: {}", pfstate); let s: Spin; @@ -715,7 +746,7 @@ mod tests { let hop: (usize, usize, Spin) = (initial_index, final_index, s); println!("------------- Updated State Long way ----------------"); - let mut pfstate2 = construct_matrix_a_from_state(¶ms, state2); + let mut pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); println!("------------- Proposed Update ------------------"); println!("Jumps from: {}, Lands on: {}", initial_index, final_index); @@ -750,6 +781,26 @@ mod tests { #[test] fn test_pfaffian_8sites_u8() { const SIZE: usize = 8; + let sys = crate::SysParams { + size: SIZE, + nelec: 0, + array_size: (SIZE + 7) / 8, + cons_t: -1.0, + cons_u: 1.0, + nfij: 4*SIZE*SIZE, + nvij: SIZE*(SIZE-1)/2, + ngi: SIZE, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &[], + hopping_bitmask: &[], + clean_update_frequency: 0, + nmcwarmup: 0, + nmcsample: 0, + tolerance_sherman_morrison: 0.0, + tolerance_singularity: 0.0, + pair_wavefunction: false, + }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij params[7 + SIZE * 7] = 1.0; @@ -765,7 +816,7 @@ mod tests { spin_down: 3u8, n_sites: SIZE, }; - let pfstate = construct_matrix_a_from_state(¶ms, state); + let pfstate = construct_matrix_a_from_state(¶ms, state, &sys); println!("Inverse Matrix: {}", pfstate); close(pfstate.pfaff, 0.04, 1e-12); } @@ -773,6 +824,26 @@ mod tests { #[test] fn test_pfaffian_8sites_u8_update_spin_up() { const SIZE: usize = 8; + let sys = crate::SysParams { + size: SIZE, + nelec: 0, + array_size: (SIZE + 7) / 8, + cons_t: -1.0, + cons_u: 1.0, + nfij: 4*SIZE*SIZE, + nvij: SIZE*(SIZE-1)/2, + ngi: SIZE, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &[], + hopping_bitmask: &[], + clean_update_frequency: 0, + nmcwarmup: 0, + nmcsample: 0, + tolerance_sherman_morrison: 0.0, + tolerance_singularity: 0.0, + pair_wavefunction: false, + }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij params[7 + SIZE * 7] = 1.0; @@ -798,7 +869,7 @@ mod tests { spin_down: 3u8, n_sites: SIZE, }; - let pfstate = construct_matrix_a_from_state(¶ms, state); + let pfstate = construct_matrix_a_from_state(¶ms, state, &sys); println!("Inverse Matrix: {}", pfstate); close(pfstate.pfaff, 0.04, 1e-12); let state2 = crate::FockState { @@ -806,7 +877,7 @@ mod tests { spin_down: 3u8, n_sites: SIZE, }; - let pfstate2 = construct_matrix_a_from_state(¶ms, state2); + let pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); let pfaff_ratio = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Up, ¶ms).0; close(pfstate.pfaff * pfaff_ratio, pfstate2.pfaff, 1e-12); @@ -815,6 +886,26 @@ mod tests { #[test] fn test_pfaffian_8sites_u8_update_spin_down() { const SIZE: usize = 8; + let sys = crate::SysParams { + size: SIZE, + nelec: 0, + array_size: (SIZE + 7) / 8, + cons_t: -1.0, + cons_u: 1.0, + nfij: 4*SIZE*SIZE, + nvij: SIZE*(SIZE-1)/2, + ngi: SIZE, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &[], + hopping_bitmask: &[], + clean_update_frequency: 0, + nmcwarmup: 0, + nmcsample: 0, + tolerance_sherman_morrison: 0.0, + tolerance_singularity: 0.0, + pair_wavefunction: false, + }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij params[7 + SIZE * 7] = 1.0; @@ -856,7 +947,7 @@ mod tests { spin_down: 3u8, n_sites: SIZE, }; - let pfstate = construct_matrix_a_from_state(¶ms, state); + let pfstate = construct_matrix_a_from_state(¶ms, state, &sys); println!("Inverse Matrix: {}", pfstate); close(pfstate.pfaff, 0.84, 1e-12); let state2 = crate::FockState { @@ -864,7 +955,7 @@ mod tests { spin_down: 5u8, n_sites: SIZE, }; - let pfstate2 = construct_matrix_a_from_state(¶ms, state2); + let pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); let tmp = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Down, ¶ms); println!("B: {:?}", tmp.1); @@ -875,6 +966,26 @@ mod tests { #[test] fn test_pfaffian_8sites_u8_update_matrix() { const SIZE: usize = 8; + let sys = crate::SysParams { + size: SIZE, + nelec: 0, + array_size: (SIZE + 7) / 8, + cons_t: -1.0, + cons_u: 1.0, + nfij: 4*SIZE*SIZE, + nvij: SIZE*(SIZE-1)/2, + ngi: SIZE, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &[], + hopping_bitmask: &[], + clean_update_frequency: 0, + nmcwarmup: 0, + nmcsample: 0, + tolerance_sherman_morrison: 0.0, + tolerance_singularity: 0.0, + pair_wavefunction: false, + }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij params[7 + SIZE * 7] = 1.0; @@ -917,7 +1028,7 @@ mod tests { n_sites: SIZE, }; println!("------------- Initial State ----------------"); - let mut pfstate = construct_matrix_a_from_state(¶ms, state); + let mut pfstate = construct_matrix_a_from_state(¶ms, state, &sys); println!("Inverse Matrix: {}", pfstate); close(pfstate.pfaff, 0.84, 1e-12); let state2 = crate::FockState { @@ -927,7 +1038,7 @@ mod tests { }; let hop: (usize, usize, Spin) = (6, 5, Spin::Down); println!("------------- Updated State Long way ----------------"); - let pfstate2 = construct_matrix_a_from_state(¶ms, state2); + let pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); println!("------------- Proposed Update ------------------"); let tmp = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Down, ¶ms); @@ -942,4 +1053,5 @@ mod tests { close(*good, test, 1e-12); } } + } diff --git a/tests/vmc_2sites.rs b/tests/vmc_2sites.rs index 529e201..66e6be8 100644 --- a/tests/vmc_2sites.rs +++ b/tests/vmc_2sites.rs @@ -22,7 +22,7 @@ const TOL_SHERMAN: f64 = 1e-12; const TOL_SINGULARITY: f64 = 1e-12; const NFIJ: usize = 4*SIZE*SIZE; -const NVIJ: usize = SIZE*SIZE; +const NVIJ: usize = SIZE*(SIZE-1)/2; const NGI: usize = SIZE; pub const HOPPINGS: [f64; SIZE*SIZE] = [ @@ -31,7 +31,7 @@ pub const HOPPINGS: [f64; SIZE*SIZE] = [ ]; #[derive(Debug)] -enum State { +pub enum State { F3, F5, F6, @@ -55,11 +55,11 @@ fn norm(par: &VarParams) -> f64 { let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; let g0 = par.gi[0]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let a = ::exp(2.0 * g0 - 2.0 * v)*sq(::abs(f00ud - f00du)); let b = ::exp(2.0 * g1 - 2.0 * v)*sq(::abs(f11ud - f11du)); - let c = sq(::abs(f01uu - f10uu)); - let d = sq(::abs(f01dd - f10dd)); + let _c = sq(::abs(f01uu - f10uu)); + let _d = sq(::abs(f01dd - f10dd)); let e = sq(::abs(f10ud - f01du)); let f = sq(::abs(f01ud - f10du)); a + b + e + f @@ -76,7 +76,7 @@ fn print_ratios(par: &VarParams) { let f10du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; let g0 = par.gi[0]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let _psi5 = (f11ud - f11du) * ::exp(g1 - v); let _psi6 = f10ud - f01du; let _psi9 = f01ud - f10du; @@ -103,7 +103,7 @@ fn print_ip(par: &VarParams) { let f10du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; let g0 = par.gi[0]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let _psi5 = (f11ud - f11du) * ::exp(g1 - v); let _psi6 = f10ud - f01du; let _psi9 = f01ud - f10du; @@ -115,6 +115,45 @@ fn print_ip(par: &VarParams) { //statesipfp.write(&format!("<10|psi> = {}\n", psi10).as_bytes()).unwrap(); } +fn individual_state(state: &State, par: &VarParams) -> f64 { + match state { + State::F3 => { + let f01dd = par.fij[0 + 1 * SIZE + 3 * SIZE * SIZE]; + let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; + f01dd - f10dd + }, + State::F5 => { + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g1 = par.gi[1]; + let v = par.vij[0]; + (f11ud - f11du) * ::exp(g1 - v) + }, + State::F6 => { + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + f10ud - f01du + }, + State::F9 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + f01ud - f10du + }, + State::F10 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let v = par.vij[0]; + (f00ud - f00du) * ::exp(g0 - v) + }, + State::F12 => { + let f01uu = par.fij[0 + 1 * SIZE + 0 * SIZE * SIZE]; + let f10uu = par.fij[1 + 0 * SIZE + 0 * SIZE * SIZE]; + f01uu - f10uu + }, + } +} + fn energy_individual_state(state: &State, par: &VarParams) -> f64 { match state { State::F3 => { @@ -128,7 +167,7 @@ fn energy_individual_state(state: &State, par: &VarParams) -> f64 { let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let a = CONS_T * (f01ud + f10ud - f01du - f10du); let b = CONS_U * (f11ud - f11du) * ::exp(g1 - v); println!("|5> pot: {}, kin: {}", b, a); @@ -141,7 +180,7 @@ fn energy_individual_state(state: &State, par: &VarParams) -> f64 { let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; let g0 = par.gi[0]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let a = CONS_T * (f00ud - f00du) * ::exp(g0 - v); let b = CONS_T * (f11ud - f11du) * ::exp(g1 - v); println!("|6> pot: {}, kin: {}", 0, a + b); @@ -154,7 +193,7 @@ fn energy_individual_state(state: &State, par: &VarParams) -> f64 { let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; let g0 = par.gi[0]; let g1 = par.gi[1]; - let v = par.vij[1]; + let v = par.vij[0]; let a = CONS_T * (f00ud - f00du) * ::exp(g0 - v); let b = CONS_T * (f11ud - f11du) * ::exp(g1 - v); println!("|9> pot: {}, kin: {}", 0, a + b); @@ -168,7 +207,7 @@ fn energy_individual_state(state: &State, par: &VarParams) -> f64 { let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; let g0 = par.gi[0]; - let v = par.vij[1]; + let v = par.vij[0]; let a = CONS_T * (f01ud + f10ud - f01du - f10du); let b = CONS_U * (f00ud - f00du) * ::exp(g0 - v); println!("|10> pot: {}, kin: {}", b, a); @@ -191,20 +230,107 @@ fn analytic(par: &VarParams) -> f64 { - par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; let b = (par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) - * ::exp(par.gi[0]-par.vij[1]); + * ::exp(par.gi[0]-par.vij[0]); let c = (par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) - * ::exp(par.gi[1]-par.vij[1]); + * ::exp(par.gi[1]-par.vij[0]); let d = 2.0 * CONS_T * (b + c) * a; let e = sq(par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) - * ::exp(2.0*par.gi[1]-2.0*par.vij[1]) * CONS_U; + * ::exp(2.0*par.gi[1]-2.0*par.vij[0]) * CONS_U; let f = sq(par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) - * ::exp(2.0*par.gi[0]-2.0*par.vij[1]) * CONS_U; + * ::exp(2.0*par.gi[0]-2.0*par.vij[0]) * CONS_U; d + e + f } +fn analytic_derivatives_expval(par: &VarParams) -> Vec { + let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; + out_der[0] = sq(individual_state(&State::F10, par)); + out_der[1] = sq(individual_state(&State::F5, par)); + out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + //out_der[3] = + //out_der[4] = + //out_der[5] = + //out_der[6] = + out_der[7] = { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * sq(individual_state(&State::F10, par)) + }; + out_der[8] = individual_state(&State::F6, par); + out_der[9] = individual_state(&State::F9, par); + out_der[10] = { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * sq(individual_state(&State::F5, par)) + }; + out_der[11] = - { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * sq(individual_state(&State::F10, par)) + }; + out_der[12] = - individual_state(&State::F9, par); + out_der[13] = - individual_state(&State::F6, par); + out_der[14] = - { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * sq(individual_state(&State::F5, par)) + }; + //out_der[3] = + //out_der[4] = + //out_der[5] = + //out_der[6] = + + out_der +} + +fn analytic_ho_expval(par: &VarParams) -> Vec { + let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; + out_der[0] = individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par); + out_der[1] = individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par); + out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + //out_der[3] = + //out_der[4] = + //out_der[5] = + //out_der[6] = + out_der[7] = { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par) + }; + out_der[8] = energy_individual_state(&State::F6, par); + out_der[9] = energy_individual_state(&State::F9, par); + out_der[10] = { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par) + }; + out_der[11] = - { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par) + }; + out_der[12] = - energy_individual_state(&State::F9, par); + out_der[13] = - energy_individual_state(&State::F6, par); + out_der[14] = - { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par) + }; + //out_der[3] = + //out_der[4] = + //out_der[5] = + //out_der[6] = + + out_der +} + +fn print_der(der1: &[f64], der2: &[f64], npar: usize) { + println!("Monte-Carlo Analytic Ratio"); + for i in 0..npar { + println!("{:11.4e} {:10.4e} {:10.4e}", der1[i], der2[i], der2[i] / der1[i]); + } +} + #[test] fn comupte_energy_from_all_states() { @@ -225,6 +351,7 @@ fn comupte_energy_from_all_states() { nvij: NVIJ, ngi: NGI, mcsample_interval: 1, + nbootstrap: 1, transfert_matrix: &HOPPINGS, hopping_bitmask: &bitmask, clean_update_frequency: CLEAN_UPDATE_FREQUENCY, @@ -232,6 +359,7 @@ fn comupte_energy_from_all_states() { nmcsample: NMCSAMP, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, + pair_wavefunction: false, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -239,7 +367,7 @@ fn comupte_energy_from_all_states() { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ]; - let mut vij = [0.0, 0.3, 0.3, 0.0]; + let mut vij = [0.3]; let mut gi = [-0.7, -0.5]; println!("fij: {:?}", fij); println!("vij: {:?}", vij); @@ -276,7 +404,7 @@ fn comupte_energy_from_all_states() { let mut mean_energy = 0.0; for i in 0..4 { let state = all_states[i]; - let pstate = construct_matrix_a_from_state(¶meters.fij, state); + let pstate = construct_matrix_a_from_state(¶meters.fij, state, &sys); println!("X^[-1] = {}", pstate); println!("Pfaffian: {}", pstate.pfaff); let proj = compute_jastrow_exp(state, ¶meters.vij, SIZE) @@ -320,7 +448,7 @@ fn comupte_energy_from_all_states() { tmp }; - let (mc_mean_energy, accumulated_states, cor) = compute_mean_energy(&mut rng, initial_state, ¶meters, &sys, &mut der); + let (mc_mean_energy, accumulated_states, error, cor) = compute_mean_energy(&mut rng, initial_state, ¶meters, &sys, &mut der); let mut out_str: String = String::new(); for s in accumulated_states.iter() { out_str.push_str(&format!("{}\n", s)); @@ -328,7 +456,6 @@ fn comupte_energy_from_all_states() { //statesfp.write(out_str.as_bytes()).unwrap(); println!("Correlation time: {}", cor); - let error = ::sqrt((1.0 + 2.0 * cor) / (NMCSAMP as f64)); let mut energy_str: String = String::new(); energy_str.push_str(&format!("{} {} {}\n", mean_energy, mc_mean_energy, error)); //energyfp.write(energy_str.as_bytes()).unwrap(); @@ -336,4 +463,21 @@ fn comupte_energy_from_all_states() { mean_energy = mean_energy / norm(¶meters); println!("Monte-Carlo: {}, Analytic: {}", mc_mean_energy, mean_energy); close(mc_mean_energy, mean_energy, mean_energy * 1e-2); + + // Test derivatives + let exp_val = analytic_derivatives_expval(¶meters); + print_der(der.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + let psi = norm(¶meters); + println!("Norm: {:10.4e}", psi); + for i in 0..sys.ngi+sys.nvij+sys.nfij { + close(der.expval_o[i] * psi, exp_val[i], 1e-2); + } + + let exp_val_ho = analytic_ho_expval(¶meters); + print_der(der.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + let psi = norm(¶meters); + println!("Norm: {:10.4e}", psi); + for i in 0..sys.ngi+sys.nvij+sys.nfij { + close(der.ho[i] * psi, exp_val_ho[i], 2e-2); + } } diff --git a/tests/vmc_2sites_analytic_conv.rs b/tests/vmc_2sites_analytic_conv.rs new file mode 100644 index 0000000..7ea254f --- /dev/null +++ b/tests/vmc_2sites_analytic_conv.rs @@ -0,0 +1,485 @@ +use assert::close; +use blas::daxpy; +use blas::dcopy; +use impurity::monte_carlo::compute_mean_energy_exact; +use impurity::gutzwiller::compute_gutzwiller_exp; +use impurity::jastrow::compute_jastrow_exp; +use impurity::pfaffian::construct_matrix_a_from_state; +use impurity::{generate_bitmask, mapto_pairwf, DerivativeOperator, FockState, SysParams, VarParams}; +use impurity::hamiltonian::{kinetic, potential}; + +// Number of sites +const SIZE: usize = 2; +// Hubbard's model $U$ parameter +const CONS_U: f64 = 8.0; +// Hubbard's model $t$ parameter +const CONS_T: f64 = -1.0; +// Number of electrons +const NELEC: usize = 2; +const NMCSAMP: usize = 10000; +const NMCWARMUP: usize = 1000; +const CLEAN_UPDATE_FREQUENCY: usize = 32; +const TOL_SHERMAN: f64 = 1e-12; +const TOL_SINGULARITY: f64 = 1e-12; +const _NOPTITER: usize = 100; +const _EPSILON_SHIFT: f64 = 0.01; +const _DT: f64 = -0.01; + +const NFIJ: usize = 4*SIZE*SIZE; +const NVIJ: usize = SIZE*(SIZE-1)/2; +const NGI: usize = SIZE; + +pub const HOPPINGS: [f64; SIZE*SIZE] = [ + 0.0, 1.0, + 1.0, 0.0, +]; + +#[derive(Debug)] +pub enum State { + F3, + F5, + F6, + F9, + F10, + F12 +} + +fn individual_state(state: &State, par: &VarParams) -> f64 { + match state { + State::F3 => { + let f01dd = par.fij[0 + 1 * SIZE + 3 * SIZE * SIZE]; + let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; + f01dd - f10dd + }, + State::F5 => { + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g1 = par.gi[1]; + let v = par.vij[0]; + (f11ud - f11du) * ::exp(g1 - v) + }, + State::F6 => { + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + f10ud - f01du + }, + State::F9 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + f01ud - f10du + }, + State::F10 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let v = par.vij[0]; + (f00ud - f00du) * ::exp(g0 - v) + }, + State::F12 => { + let f01uu = par.fij[0 + 1 * SIZE + 0 * SIZE * SIZE]; + let f10uu = par.fij[1 + 0 * SIZE + 0 * SIZE * SIZE]; + f01uu - f10uu + }, + } +} + +fn norm(par: &VarParams) -> f64 { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f01uu = par.fij[0 + 1 * SIZE + 0 * SIZE * SIZE]; + let f10uu = par.fij[1 + 0 * SIZE + 0 * SIZE * SIZE]; + let f01dd = par.fij[0 + 1 * SIZE + 3 * SIZE * SIZE]; + let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = ::exp(2.0 * g0 - 2.0 * v)*sq(::abs(f00ud - f00du)); + let b = ::exp(2.0 * g1 - 2.0 * v)*sq(::abs(f11ud - f11du)); + let _c = sq(::abs(f01uu - f10uu)); + let _d = sq(::abs(f01dd - f10dd)); + let e = sq(::abs(f10ud - f01du)); + let f = sq(::abs(f01ud - f10du)); + a + b + e + f +} + +fn print_ratios(par: &VarParams) { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let f01ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let _psi5 = (f11ud - f11du) * ::exp(g1 - v); + let _psi6 = f10ud - f01du; + let _psi9 = f01ud - f10du; + let _psi10 = (f00ud - f00du) * ::exp(g0 - v); + //let mut statesfp = File::create("ratios_th").unwrap(); + //statesfp.write(&format!("<5|psi>/<6|psi> = {}\n", sq(psi5 / psi6)).as_bytes()).unwrap(); + //statesfp.write(&format!("<5|psi>/<9|psi> = {}\n", sq(psi5 / psi9)).as_bytes()).unwrap(); + //statesfp.write(&format!("<10|psi>/<6|psi> = {}\n", sq(psi10 / psi6)).as_bytes()).unwrap(); + //statesfp.write(&format!("<10|psi>/<9|psi> = {}\n", sq(psi10 / psi9)).as_bytes()).unwrap(); + //statesfp.write(&format!("<6|psi>/<5|psi> = {}\n", sq(psi6 / psi5)).as_bytes()).unwrap(); + //statesfp.write(&format!("<6|psi>/<10|psi> = {}\n", sq(psi6 / psi10)).as_bytes()).unwrap(); + //statesfp.write(&format!("<9|psi>/<5|psi> = {}\n", sq(psi9 / psi5)).as_bytes()).unwrap(); + //statesfp.write(&format!("<9|psi>/<10|psi> = {}\n", sq(psi9 / psi10)).as_bytes()).unwrap(); +} + +fn print_ip(par: &VarParams) { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let f01ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let _psi5 = (f11ud - f11du) * ::exp(g1 - v); + let _psi6 = f10ud - f01du; + let _psi9 = f01ud - f10du; + let _psi10 = (f00ud - f00du) * ::exp(g0 - v); + //let mut statesipfp = File::create("statesip").unwrap(); + //statesipfp.write(&format!("<5|psi> = {}\n", psi5).as_bytes()).unwrap(); + //statesipfp.write(&format!("<9|psi> = {}\n", psi9).as_bytes()).unwrap(); + //statesipfp.write(&format!("<6|psi> = {}\n", psi6).as_bytes()).unwrap(); + //statesipfp.write(&format!("<10|psi> = {}\n", psi10).as_bytes()).unwrap(); +} + +fn energy_individual_state(state: &State, par: &VarParams) -> f64 { + match state { + State::F3 => { + 0.0 + }, + State::F5 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = CONS_T * (f01ud + f10ud - f01du - f10du); + let b = CONS_U * (f11ud - f11du) * ::exp(g1 - v); + println!("|5> pot: {}, kin: {}", b, a); + a + b + }, + State::F6 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = CONS_T * (f00ud - f00du) * ::exp(g0 - v); + let b = CONS_T * (f11ud - f11du) * ::exp(g1 - v); + println!("|6> pot: {}, kin: {}", 0, a + b); + a + b + }, + State::F9 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = CONS_T * (f00ud - f00du) * ::exp(g0 - v); + let b = CONS_T * (f11ud - f11du) * ::exp(g1 - v); + println!("|9> pot: {}, kin: {}", 0, a + b); + a + b + }, + State::F10 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let v = par.vij[0]; + let a = CONS_T * (f01ud + f10ud - f01du - f10du); + let b = CONS_U * (f00ud - f00du) * ::exp(g0 - v); + println!("|10> pot: {}, kin: {}", b, a); + a + b + }, + State::F12 => { + 0.0 + }, + } +} + +fn sq(a: f64) -> f64 { + a*a +} + +fn analytic(par: &VarParams) -> f64 { + let a = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE] + - par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE] + + par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE] + - par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let b = (par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] + - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) + * ::exp(par.gi[0]-par.vij[0]); + let c = (par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] + - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) + * ::exp(par.gi[1]-par.vij[0]); + let d = 2.0 * CONS_T * (b + c) * a; + let e = sq(par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] + - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) + * ::exp(2.0*par.gi[1]-2.0*par.vij[0]) * CONS_U; + let f = sq(par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] + - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) + * ::exp(2.0*par.gi[0]-2.0*par.vij[0]) * CONS_U; + d + e + f +} + +fn analytic_ho_expval(par: &VarParams) -> Vec { + let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; + out_der[0] = individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par); + out_der[1] = individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par); + out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + //out_der[3] = + //out_der[4] = + //out_der[5] = + //out_der[6] = + out_der[3] = 1.0 * { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par) + }; + out_der[4] = 1.0 * energy_individual_state(&State::F6, par); + out_der[5] = 1.0 * energy_individual_state(&State::F9, par); + out_der[6] = 1.0 * { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par) + }; + + out_der +} + +fn analytic_derivatives_expval(par: &VarParams) -> Vec { + let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; + out_der[0] = sq(individual_state(&State::F10, par)); + out_der[1] = sq(individual_state(&State::F5, par)); + out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + // fij + out_der[3] = 1.0 * { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * sq(individual_state(&State::F10, par)) + }; + out_der[4] = 1.0 * individual_state(&State::F6, par); + out_der[5] = 1.0 * individual_state(&State::F9, par); + out_der[6] = 1.0 * { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * sq(individual_state(&State::F5, par)) + }; + out_der +} + +fn print_der(der1: &[f64], der2: &[f64], npar: usize) { + println!("Monte-Carlo Analytic Ratio"); + for i in 0..npar { + println!("{:11.4e} {:10.4e} {:10.4e}", der1[i], der2[i], der2[i] / der1[i]); + } +} + + +#[test] +fn comupte_energy_from_all_states() { + //let mut statesfp = File::create("states").unwrap(); + //let mut energyfp = File::create("energy").unwrap(); + env_logger::init(); + //let mut rng = thread_rng(); + let bitmask = generate_bitmask(&HOPPINGS, SIZE); + println!("bitmasks: {:?}", bitmask); + let sys = SysParams { + size: SIZE, + nelec: NELEC, + array_size: (SIZE + 7) / 8, + cons_t: CONS_T, + cons_u: CONS_U, + nfij: SIZE*SIZE, + nvij: NVIJ, + ngi: NGI, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &HOPPINGS, + hopping_bitmask: &bitmask, + clean_update_frequency: CLEAN_UPDATE_FREQUENCY, + nmcwarmup: NMCWARMUP, + nmcsample: NMCSAMP, + tolerance_sherman_morrison: TOL_SHERMAN, + tolerance_singularity: TOL_SINGULARITY, + pair_wavefunction: false, + }; + let mut fij = [ + 0.0, 0.0, 0.0, 0.0, + 1.093500753438337580e-01, 3.768419990611672210e-01, 3.769186909982900624e-01, 3.322533463612635796e-01, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + ]; + let mut vij = [5.079558854017672820e-01]; + let mut gi = [3.016937239100276336e-01, -8.096496093117950821e-01]; + println!("fij: {:?}", fij); + println!("vij: {:?}", vij); + println!("gi: {:?}", gi); + let parameters = VarParams { + size: SIZE, + fij: &mut fij, + gi: &mut gi, + vij: &mut vij + }; + print_ratios(¶meters); + print_ip(¶meters); + + // Generate all the states for 2 sites S=0 + let states_names: [State; 4] = [State::F10, State::F9, State::F6, State::F5]; + let all_states: [FockState; 4] = [ + // \ket{10} + FockState{spin_up: 128u8, spin_down: 128u8, n_sites: SIZE}, + // \ket{12} + //FockState{spin_up: 192u8, spin_down: 0u8, n_sites: SIZE}, + // \ket{3} + //FockState{spin_up: 0u8, spin_down: 192u8, n_sites: SIZE}, + // \ket{9} + FockState{spin_up: 128u8, spin_down: 64u8, n_sites: SIZE}, + // \ket{6} + FockState{spin_up: 64u8, spin_down: 128u8, n_sites: SIZE}, + // \ket{5} + FockState{spin_up: 64u8, spin_down: 64u8, n_sites: SIZE}, + ]; + println!("Basis states."); + for i in 0..4 { + println!("State {} : {}", i, all_states[i]); + } + let mut mean_energy = 0.0; + for i in 0..4 { + let state = all_states[i]; + let pstate = construct_matrix_a_from_state(¶meters.fij, state, &sys); + println!("X^[-1] = {}", pstate); + println!("Pfaffian: {}", pstate.pfaff); + let proj = compute_jastrow_exp(state, ¶meters.vij, SIZE) + +compute_gutzwiller_exp(state, ¶meters.gi, SIZE); + println!("Proj: {}", proj); + let ip = pstate.pfaff * ::exp(proj); + println!("<{:?}|psi> = {}", states_names[i], ip); + let pot = potential(state, proj, &pstate, &sys); + println!("Potential energy: {}", pot); + let kin = kinetic(state, &pstate, proj, ¶meters, &sys); + println!("Kinetic energy: {}", kin); + energy_individual_state(&states_names[i], ¶meters); + mean_energy += (kin + pot) * pstate.pfaff * ::exp(proj); + println!("Energy: {}", (kin + pot) / norm(¶meters)); + close( kin + pot, + energy_individual_state(&states_names[i], ¶meters), 1e-12); + } + close(mean_energy, analytic(¶meters), 1e-12); + mean_energy = mean_energy / norm(¶meters); + let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut visited: Vec = vec![0; NMCSAMP]; + let mut der = DerivativeOperator { + o_tilde: &mut otilde, + expval_o: &mut expvalo, + ho: &mut expval_ho, + n: (NFIJ + NVIJ + NGI) as i32, + nsamp: NMCSAMP as f64, + nsamp_int: 1, + mu: -1, + visited: &mut visited, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + epsilon: 0.0, + }; + let mut otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let mut expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut visited_pair: Vec = vec![0; NMCSAMP]; + let mut der_pair = DerivativeOperator { + o_tilde: &mut otilde_pair, + expval_o: &mut expvalo_pair, + ho: &mut expval_ho_pair, + n: (sys.size*sys.size + NVIJ + NGI) as i32, + nsamp: NMCSAMP as f64, + nsamp_int: 1, + mu: -1, + visited: &mut visited_pair, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + epsilon: 0.0, + }; + + let mean_energy_es = compute_mean_energy_exact(¶meters, &sys, &mut der); + der_pair.mu = der.mu; + //for s in accumulated_states.iter() { + // out_str.push_str(&format!("{}\n", s)); + //} + //statesfp.write(out_str.as_bytes()).unwrap(); + + //println!("Correlation time: {}", cor); + //let error = ::sqrt((1.0 + 2.0 * cor) / (NMCSAMP as f64)); + let mut energy_str: String = String::new(); + energy_str.push_str(&format!("{} {}\n", mean_energy, mean_energy_es)); + //energyfp.write(energy_str.as_bytes()).unwrap(); + //println!("Comparing monte-carlo energy, tol: {}", error); + println!("Monte-Carlo: {}, Analytic: {}", mean_energy_es, mean_energy); + close(mean_energy_es, mean_energy, 1e-12); + mapto_pairwf(&der, &mut der_pair, &sys); + + // Test derivatives + let exp_val = analytic_derivatives_expval(¶meters); + println!("Checking "); + print_der(der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + let psi = norm(¶meters); + println!("Norm: {:10.4e}", psi); + for i in 0..sys.ngi+sys.nvij+sys.nfij { + println!("{} == {}, tol = {}", der_pair.expval_o[i], exp_val[i] / psi, 2e-2); + close(der_pair.expval_o[i], exp_val[i] / psi, 1e-12); + } + + let exp_val_ho = analytic_ho_expval(¶meters); + println!("Checking "); + print_der(der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + let psi = norm(¶meters); + println!("Norm: {:10.4e}", psi); + for i in 0..sys.ngi+sys.nvij+sys.nfij { + println!("{} == {}, tol = {}", der_pair.ho[i] * psi, exp_val_ho[i], 2e-2); + close(der_pair.ho[i] * psi, exp_val_ho[i], 2e-2); + } + let mut x0 = vec![0.0; SIZE * SIZE + NVIJ + NGI]; + x0[0] = parameters.gi[0]; + x0[1] = parameters.gi[1]; + x0[2] = parameters.vij[0]; + for i in 0.. SIZE * SIZE { + x0[NGI + NVIJ + i] = parameters.fij[SIZE*SIZE + i]; + } + + // 68 misawa + let mut b: Vec = vec![0.0; der_pair.n as usize]; + unsafe { + let incx = 1; + let incy = 1; + daxpy(der_pair.n, -mean_energy, der_pair.expval_o, incx, der_pair.ho, incy); + dcopy(der_pair.n, der_pair.ho, incx, &mut b, incy); + } + println!("x0 = {:?}", x0); +} diff --git a/tests/vmc_2sites_exact_sum.rs b/tests/vmc_2sites_exact_sum.rs new file mode 100644 index 0000000..048d67a --- /dev/null +++ b/tests/vmc_2sites_exact_sum.rs @@ -0,0 +1,464 @@ +use assert::close; +use impurity::monte_carlo::compute_mean_energy_exact; +use impurity::gutzwiller::compute_gutzwiller_exp; +use impurity::jastrow::compute_jastrow_exp; +use impurity::pfaffian::construct_matrix_a_from_state; +use impurity::{generate_bitmask, mapto_pairwf, DerivativeOperator, FockState, SysParams, VarParams}; +use impurity::hamiltonian::{kinetic, potential}; + +// Number of sites +const SIZE: usize = 2; +// Hubbard's model $U$ parameter +const CONS_U: f64 = 8.0; +// Hubbard's model $t$ parameter +const CONS_T: f64 = -1.0; +// Number of electrons +const NELEC: usize = 2; +const NMCSAMP: usize = 10000; +const NMCWARMUP: usize = 1000; +const CLEAN_UPDATE_FREQUENCY: usize = 32; +const TOL_SHERMAN: f64 = 1e-12; +const TOL_SINGULARITY: f64 = 1e-12; + +const NFIJ: usize = 4*SIZE*SIZE; +const NVIJ: usize = SIZE*(SIZE-1)/2; +const NGI: usize = SIZE; + +pub const HOPPINGS: [f64; SIZE*SIZE] = [ + 0.0, 1.0, + 1.0, 0.0, +]; + +#[derive(Debug)] +pub enum State { + F3, + F5, + F6, + F9, + F10, + F12 +} + +fn individual_state(state: &State, par: &VarParams) -> f64 { + match state { + State::F3 => { + let f01dd = par.fij[0 + 1 * SIZE + 3 * SIZE * SIZE]; + let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; + f01dd - f10dd + }, + State::F5 => { + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g1 = par.gi[1]; + let v = par.vij[0]; + (f11ud - f11du) * ::exp(g1 - v) + }, + State::F6 => { + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + f10ud - f01du + }, + State::F9 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + f01ud - f10du + }, + State::F10 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let v = par.vij[0]; + (f00ud - f00du) * ::exp(g0 - v) + }, + State::F12 => { + let f01uu = par.fij[0 + 1 * SIZE + 0 * SIZE * SIZE]; + let f10uu = par.fij[1 + 0 * SIZE + 0 * SIZE * SIZE]; + f01uu - f10uu + }, + } +} + +fn norm(par: &VarParams) -> f64 { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f01uu = par.fij[0 + 1 * SIZE + 0 * SIZE * SIZE]; + let f10uu = par.fij[1 + 0 * SIZE + 0 * SIZE * SIZE]; + let f01dd = par.fij[0 + 1 * SIZE + 3 * SIZE * SIZE]; + let f10dd = par.fij[1 + 0 * SIZE + 3 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = ::exp(2.0 * g0 - 2.0 * v)*sq(::abs(f00ud - f00du)); + let b = ::exp(2.0 * g1 - 2.0 * v)*sq(::abs(f11ud - f11du)); + let _c = sq(::abs(f01uu - f10uu)); + let _d = sq(::abs(f01dd - f10dd)); + let e = sq(::abs(f10ud - f01du)); + let f = sq(::abs(f01ud - f10du)); + a + b + e + f +} + +fn print_ratios(par: &VarParams) { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let f01ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let _psi5 = (f11ud - f11du) * ::exp(g1 - v); + let _psi6 = f10ud - f01du; + let _psi9 = f01ud - f10du; + let _psi10 = (f00ud - f00du) * ::exp(g0 - v); + //let mut statesfp = File::create("ratios_th").unwrap(); + //statesfp.write(&format!("<5|psi>/<6|psi> = {}\n", sq(psi5 / psi6)).as_bytes()).unwrap(); + //statesfp.write(&format!("<5|psi>/<9|psi> = {}\n", sq(psi5 / psi9)).as_bytes()).unwrap(); + //statesfp.write(&format!("<10|psi>/<6|psi> = {}\n", sq(psi10 / psi6)).as_bytes()).unwrap(); + //statesfp.write(&format!("<10|psi>/<9|psi> = {}\n", sq(psi10 / psi9)).as_bytes()).unwrap(); + //statesfp.write(&format!("<6|psi>/<5|psi> = {}\n", sq(psi6 / psi5)).as_bytes()).unwrap(); + //statesfp.write(&format!("<6|psi>/<10|psi> = {}\n", sq(psi6 / psi10)).as_bytes()).unwrap(); + //statesfp.write(&format!("<9|psi>/<5|psi> = {}\n", sq(psi9 / psi5)).as_bytes()).unwrap(); + //statesfp.write(&format!("<9|psi>/<10|psi> = {}\n", sq(psi9 / psi10)).as_bytes()).unwrap(); +} + +fn print_ip(par: &VarParams) { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let f01ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let _psi5 = (f11ud - f11du) * ::exp(g1 - v); + let _psi6 = f10ud - f01du; + let _psi9 = f01ud - f10du; + let _psi10 = (f00ud - f00du) * ::exp(g0 - v); + //let mut statesipfp = File::create("statesip").unwrap(); + //statesipfp.write(&format!("<5|psi> = {}\n", psi5).as_bytes()).unwrap(); + //statesipfp.write(&format!("<9|psi> = {}\n", psi9).as_bytes()).unwrap(); + //statesipfp.write(&format!("<6|psi> = {}\n", psi6).as_bytes()).unwrap(); + //statesipfp.write(&format!("<10|psi> = {}\n", psi10).as_bytes()).unwrap(); +} + +fn energy_individual_state(state: &State, par: &VarParams) -> f64 { + match state { + State::F3 => { + 0.0 + }, + State::F5 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = CONS_T * (f01ud + f10ud - f01du - f10du); + let b = CONS_U * (f11ud - f11du) * ::exp(g1 - v); + println!("|5> pot: {}, kin: {}", b, a); + a + b + }, + State::F6 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = CONS_T * (f00ud - f00du) * ::exp(g0 - v); + let b = CONS_T * (f11ud - f11du) * ::exp(g1 - v); + println!("|6> pot: {}, kin: {}", 0, a + b); + a + b + }, + State::F9 => { + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let f11ud = par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE]; + let f11du = par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let g1 = par.gi[1]; + let v = par.vij[0]; + let a = CONS_T * (f00ud - f00du) * ::exp(g0 - v); + let b = CONS_T * (f11ud - f11du) * ::exp(g1 - v); + println!("|9> pot: {}, kin: {}", 0, a + b); + a + b + }, + State::F10 => { + let f01ud = par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE]; + let f10ud = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE]; + let f01du = par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE]; + let f10du = par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let f00ud = par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE]; + let f00du = par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]; + let g0 = par.gi[0]; + let v = par.vij[0]; + let a = CONS_T * (f01ud + f10ud - f01du - f10du); + let b = CONS_U * (f00ud - f00du) * ::exp(g0 - v); + println!("|10> pot: {}, kin: {}", b, a); + a + b + }, + State::F12 => { + 0.0 + }, + } +} + +fn sq(a: f64) -> f64 { + a*a +} + +fn analytic(par: &VarParams) -> f64 { + let a = par.fij[1 + 0 * SIZE + 1 * SIZE * SIZE] + - par.fij[0 + 1 * SIZE + 2 * SIZE * SIZE] + + par.fij[0 + 1 * SIZE + 1 * SIZE * SIZE] + - par.fij[1 + 0 * SIZE + 2 * SIZE * SIZE]; + let b = (par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] + - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) + * ::exp(par.gi[0]-par.vij[0]); + let c = (par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] + - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) + * ::exp(par.gi[1]-par.vij[0]); + let d = 2.0 * CONS_T * (b + c) * a; + let e = sq(par.fij[1 + 1 * SIZE + 1 * SIZE * SIZE] + - par.fij[1 + 1 * SIZE + 2 * SIZE * SIZE]) + * ::exp(2.0*par.gi[1]-2.0*par.vij[0]) * CONS_U; + let f = sq(par.fij[0 + 0 * SIZE + 1 * SIZE * SIZE] + - par.fij[0 + 0 * SIZE + 2 * SIZE * SIZE]) + * ::exp(2.0*par.gi[0]-2.0*par.vij[0]) * CONS_U; + d + e + f +} + +fn analytic_ho_expval(par: &VarParams) -> Vec { + let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; + out_der[0] = individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par); + out_der[1] = individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par); + out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + //out_der[3] = + //out_der[4] = + //out_der[5] = + //out_der[6] = + out_der[3] = 1.0 * { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * individual_state(&State::F10, par) + * energy_individual_state(&State::F10, par) + }; + out_der[4] = 1.0 * energy_individual_state(&State::F6, par); + out_der[5] = 1.0 * energy_individual_state(&State::F9, par); + out_der[6] = 1.0 * { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * individual_state(&State::F5, par) + * energy_individual_state(&State::F5, par) + }; + + out_der +} + +fn analytic_derivatives_expval(par: &VarParams) -> Vec { + let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; + out_der[0] = sq(individual_state(&State::F10, par)); + out_der[1] = sq(individual_state(&State::F5, par)); + out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + // fij + out_der[3] = 1.0 * { + ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) + * sq(individual_state(&State::F10, par)) + }; + out_der[4] = 1.0 * individual_state(&State::F6, par); + out_der[5] = 1.0 * individual_state(&State::F9, par); + out_der[6] = 1.0 * { + ::exp(par.gi[1] - par.vij[0]) / individual_state(&State::F5, par) + * sq(individual_state(&State::F5, par)) + }; + out_der +} + +fn print_der(der1: &[f64], der2: &[f64], npar: usize) { + println!("Monte-Carlo Analytic Ratio"); + for i in 0..npar { + println!("{:11.4e} {:10.4e} {:10.4e}", der1[i], der2[i], der2[i] / der1[i]); + } +} + + + +#[test] +fn comupte_energy_from_all_states() { + //let mut statesfp = File::create("states").unwrap(); + //let mut energyfp = File::create("energy").unwrap(); + env_logger::init(); + //let mut rng = thread_rng(); + let bitmask = generate_bitmask(&HOPPINGS, SIZE); + println!("bitmasks: {:?}", bitmask); + let sys = SysParams { + size: SIZE, + nelec: NELEC, + array_size: (SIZE + 7) / 8, + cons_t: CONS_T, + cons_u: CONS_U, + nfij: SIZE*SIZE, + nvij: NVIJ, + ngi: NGI, + mcsample_interval: 1, + nbootstrap: 1, + transfert_matrix: &HOPPINGS, + hopping_bitmask: &bitmask, + clean_update_frequency: CLEAN_UPDATE_FREQUENCY, + nmcwarmup: NMCWARMUP, + nmcsample: NMCSAMP, + tolerance_sherman_morrison: TOL_SHERMAN, + tolerance_singularity: TOL_SINGULARITY, + pair_wavefunction: false, + }; + let mut fij = [ + 0.0, 0.0, 0.0, 0.0, + 1.093500753438337580e-01, 3.768419990611672210e-01, 3.769186909982900624e-01, 3.322533463612635796e-01, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + ]; + let mut vij = [5.079558854017672820e-01]; + let mut gi = [3.016937239100276336e-01, -8.096496093117950821e-01]; + println!("fij: {:?}", fij); + println!("vij: {:?}", vij); + println!("gi: {:?}", gi); + let parameters = VarParams { + size: SIZE, + fij: &mut fij, + gi: &mut gi, + vij: &mut vij + }; + print_ratios(¶meters); + print_ip(¶meters); + + // Generate all the states for 2 sites S=0 + let states_names: [State; 4] = [State::F10, State::F9, State::F6, State::F5]; + let all_states: [FockState; 4] = [ + // \ket{10} + FockState{spin_up: 128u8, spin_down: 128u8, n_sites: SIZE}, + // \ket{12} + //FockState{spin_up: 192u8, spin_down: 0u8, n_sites: SIZE}, + // \ket{3} + //FockState{spin_up: 0u8, spin_down: 192u8, n_sites: SIZE}, + // \ket{9} + FockState{spin_up: 128u8, spin_down: 64u8, n_sites: SIZE}, + // \ket{6} + FockState{spin_up: 64u8, spin_down: 128u8, n_sites: SIZE}, + // \ket{5} + FockState{spin_up: 64u8, spin_down: 64u8, n_sites: SIZE}, + ]; + println!("Basis states."); + for i in 0..4 { + println!("State {} : {}", i, all_states[i]); + } + let mut mean_energy = 0.0; + for i in 0..4 { + let state = all_states[i]; + let pstate = construct_matrix_a_from_state(¶meters.fij, state, &sys); + println!("X^[-1] = {}", pstate); + println!("Pfaffian: {}", pstate.pfaff); + let proj = compute_jastrow_exp(state, ¶meters.vij, SIZE) + +compute_gutzwiller_exp(state, ¶meters.gi, SIZE); + println!("Proj: {}", proj); + let ip = pstate.pfaff * ::exp(proj); + println!("<{:?}|psi> = {}", states_names[i], ip); + let pot = potential(state, proj, &pstate, &sys); + println!("Potential energy: {}", pot); + let kin = kinetic(state, &pstate, proj, ¶meters, &sys); + println!("Kinetic energy: {}", kin); + energy_individual_state(&states_names[i], ¶meters); + mean_energy += (kin + pot) * pstate.pfaff * ::exp(proj); + println!("Energy: {}", (kin + pot) / norm(¶meters)); + close( kin + pot, + energy_individual_state(&states_names[i], ¶meters), 1e-12); + } + close(mean_energy, analytic(¶meters), 1e-12); + mean_energy = mean_energy / norm(¶meters); + let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut visited: Vec = vec![0; NMCSAMP]; + let mut der = DerivativeOperator { + o_tilde: &mut otilde, + expval_o: &mut expvalo, + ho: &mut expval_ho, + n: (NFIJ + NVIJ + NGI) as i32, + nsamp: NMCSAMP as f64, + nsamp_int: 1, + mu: -1, + visited: &mut visited, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + epsilon: 0.0, + }; + let mut otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let mut expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut visited_pair: Vec = vec![0; NMCSAMP]; + let mut der_pair = DerivativeOperator { + o_tilde: &mut otilde_pair, + expval_o: &mut expvalo_pair, + ho: &mut expval_ho_pair, + n: (sys.size*sys.size + NVIJ + NGI) as i32, + nsamp: NMCSAMP as f64, + nsamp_int: 1, + mu: -1, + visited: &mut visited_pair, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + epsilon: 0.0, + }; + + let mean_energy_es = compute_mean_energy_exact(¶meters, &sys, &mut der); + der_pair.mu = der.mu; + //for s in accumulated_states.iter() { + // out_str.push_str(&format!("{}\n", s)); + //} + //statesfp.write(out_str.as_bytes()).unwrap(); + + //println!("Correlation time: {}", cor); + //let error = ::sqrt((1.0 + 2.0 * cor) / (NMCSAMP as f64)); + let mut energy_str: String = String::new(); + energy_str.push_str(&format!("{} {}\n", mean_energy, mean_energy_es)); + //energyfp.write(energy_str.as_bytes()).unwrap(); + //println!("Comparing monte-carlo energy, tol: {}", error); + println!("Monte-Carlo: {}, Analytic: {}", mean_energy_es, mean_energy); + close(mean_energy_es, mean_energy, 1e-16); + mapto_pairwf(&der, &mut der_pair, &sys); + + // Test derivatives + let exp_val = analytic_derivatives_expval(¶meters); + println!("Checking "); + print_der(der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + let psi = norm(¶meters); + println!("Norm: {:10.4e}", psi); + for i in 0..sys.ngi+sys.nvij+sys.nfij { + println!("{} == {}, tol = {}", der_pair.expval_o[i], exp_val[i] / psi, 1e-12); + close(der_pair.expval_o[i], exp_val[i] / psi, 1e-12); + } + + let exp_val_ho = analytic_ho_expval(¶meters); + println!("Checking "); + print_der(der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + let psi = norm(¶meters); + println!("Norm: {:10.4e}", psi); + for i in 0..sys.ngi+sys.nvij+sys.nfij { + println!("{} == {}, tol = {}", der_pair.ho[i], exp_val_ho[i] / psi, 1e-12); + close(der_pair.ho[i], exp_val_ho[i] / psi, 1e-12); + } +} diff --git a/tests/vmc_clean_updates.rs b/tests/vmc_clean_updates.rs index 3b3d6fd..855c7b2 100644 --- a/tests/vmc_clean_updates.rs +++ b/tests/vmc_clean_updates.rs @@ -1,5 +1,3 @@ -use std::ptr::addr_of; -use std::slice::from_raw_parts as slice; use log::info; use assert::close; @@ -82,6 +80,7 @@ fn monte_carlo_first_iteration() { nvij: NVIJ, ngi: NGI, mcsample_interval: 1, + nbootstrap: 1, transfert_matrix: &HOPPINGS, hopping_bitmask: &bitmask, clean_update_frequency: CLEAN_UPDATE_FREQUENCY, @@ -89,6 +88,7 @@ fn monte_carlo_first_iteration() { nmcsample: NMCSAMP, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERANCE_SINGULARITY, + pair_wavefunction: false, }; log_system_parameters(&sys); @@ -173,7 +173,7 @@ fn monte_carlo_first_iteration() { info!("Initial Nelec: {}, {}", state.spin_down.count_ones(), state.spin_up.count_ones()); info!("Nsites: {}", state.n_sites); - let (energy, _, _) = compute_mean_energy(&mut rng, state, ¶meters, &sys, &mut der); + let (energy, _, _, _) = compute_mean_energy(&mut rng, state, ¶meters, &sys, &mut der); close(energy, -0.35, MONTE_CARLO_CONVERGENCE_TOLERANCE); close(energy, mean_energy_analytic_2sites(¶meters, &sys), MONTE_CARLO_CONVERGENCE_TOLERANCE); }