text:
## This file auto-generated by Yggdrasil/build_hashes.jl
shard_hash_table = Dict{CompilerShard,String}(
    CompilerShard("Rootfs", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz)
      => "09ca76f2134de87185b87ed20fe288a66b854fad199866806208809731c51873",
    CompilerShard("Rootfs", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs)
      => "225166fd5d70020e921187e121283b2dcd0dff586ecb6caecd2190277d66cb52",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "74eb805eb56600105514e0f9f404ddbcf8ecbfb7e1277ca3b3686efba7fa395c",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "dffb39a44c1f87fdb1354dbdb2fc02ef27024de3fa38543704a6316e62d39329",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "de9b6e3e07eed4f63aa6c335ff436e1e40f17c5c24df902ee72a18f17e059478",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "ae50e18ecf75e2ce64b3d53ca29618fd71ad33d7b03f98361432b9f8c177fa3c",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "f8f730cba90255135e85fa9c80d4be84f2a8e91ae9542a63cb4f6fd76d0f3adc",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "5251aa4135abebf16cba250ce23f3b4a89dfdacf9e1e8c8f8897e549f5015f30",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "e25b8d3f03d0488d03dc2173fc5ae1de49a1afd30b70b4f6513a3296afdbd78a",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "2a94622f6af2757c6f956614e090de15ef1f0175199e3fc4e8f7a103866a26d6",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "7c9dc0124830aa33440677afc8d6e82ffb8a2049d9bdee4ed243d46420150e4b",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "7cc4f4a56b311aa13e7952c4b3735f70a2ff173786ffb5aedb0edd618abcbb83",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "d1e6e9b9c2e614eee59bbe4764391c5cdda54f58d51f0cfef2348afc367ea6aa",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "fc77a555d6c8b4fb17229d205c98567f4bd3c03a0c8774f791c3b947caf47b1e",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "f2a7f622b9ca4d8957eba19b75a5a970a78e2edbb87615cadb30f0ed11452568",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "b09731fd8ca8a8341a9eb4301f05dcef1a5d742a84d946139faa61abd67fb372",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "64582217195062887eca489604286c19e606d4e73749682fa3e0fa7b9ec0c077",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "30b8631a7c2475e0d30eca921cc1af235716d058355cd6582017a55274eeb33b",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "cf0422cac34c8015128917749c3d79148ce156d142f1cdd58d47c2619c3fe86d",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "314fd8a4054697cffa3263f18bb9524a58699bde429170abee88c8e219e101f1",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "7cc8e1f997ab684740835e3e92175dce3874975f10018cf4968e6fee968a939b",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "b78e550d1ae17ea6a5c2e9b26bc7690e43f56ab02c474c3315bc36da470ae084",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "13f5c4a836f0c0f82f2b811b77cf8379f93768fe4399cca85511c0f48f7a595a",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "700871928ed87ebc1827f66ab4a42004950fd352c315ed4f8495716d5a3829ca",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "4427c61f832aca7f23074d9e1fa95aa6d1a73847f90a3e27fca9cca68185b4da",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "8e33240653afd3e26a4d59b07a42c38d7f4babce4eb1c49f8ec45866acc5befd",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "df7a46e825a84966584f23fb5439d3c9d405a6479bd2d28e315ed7c91f842ab8",
    CompilerShard("BaseCompilerShard", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "700abed47c371033a5f745405a472aa86207ed241d23a982a0b87c59db9fa216",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "805c2ae8562a51945b8b6a6c263f20ee52d110891e7799b21093750730cb74c7",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "727be698216c43199b730c84de01b160ed4d83c2d176b07a9949bf17a9c9a9be",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "00d19deb7a9ac79505d7ffd48b6bf86a82b82ec77436751872d37b7acf53b725",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "70fda2388e9c74c0e05ff1810faca532bec4ecd49c7550f9f1b72b9ac81f797f",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "d8931c2b598562dddf5cb60fec1f39a5c1ca187d42dad744d1bdbfc7b6d9f9ba",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "9d74a42210e42b607c91ca16d94b0530825eae416e287601eea7ee312caf02ac",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "e1230183f44b3b2aa9e7462322e0fa06fc5584fec3233758a548a88fa47df074",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "49cab4b5fdb8a8eaf2d4551ab6e3f9f501eeda5909d02179e3502d9d3871cf92",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "ab93bb13f68cbd733f3196e58e0b99838201aa15f81fe110ef3ad3fbb7de4cf8",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "a65bb81dde47704221960c69e3cede330de52db88dd50ce5bad6cc7e5d16ae83",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "269e2b68de0b5ed0ae9a612b4885f8d69acc8eed0d7e44527da5a8c546ab8f3d",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "e337851acabab44a9c66fa075569a71480b6ea6d7e8e17081985e97cfa06e044",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "a775a17c89575938345a0a33eaa86fcbb58ed38f7c13072e3ad2eb166463749b",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "ce31f00c378627746f6c404064991d680f357db07cf1ced9fdf47a42f0b5fd0f",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "3c06e0dd811c3fa86a21d1286beeb47ef7c144873bf8b9d6f2dfcf130fb5e1cc",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "d926f4aa8fec3cb4c2468c3970ec9c815c3b6fcf29c46d5abe653770dce31ba6",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "90469f4863c425e5e4e488f445ecdd252db428ce67c527107e0ed4c1c9c6d4d1",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "074e108776b699b77eef886c8308d5e4b17a83549068dd0ad1f0410c8c05cf33",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "070f3bca044cc6268fbad762fade04ac84ef31620fe08812c9efaf4dc1aced2a",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "a955805f16b621021607b0d22f6d892757ab1d076120fcf718a87712d279459b",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "9b83746ecb5ac3de1b71ad4679a23cfbafb314340d1e173a3b90c5c6a55dc84e",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "40649207e3e95021c324ceca02d30f64c01c5f8f8a63e79e26061274bd36e2ab",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "dfae9a7ecdce8bbeca99c85f769074cd8da018a45e5a1dffe7aab6c6d01de22d",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "8f1ca5b0ae46555001d3097be7882295d41b8b2b641c33eeb14a9972625479d4",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "9ffb011e5ab0a95bdf73ef5c0fac6b6e8f1b943bb1bc5118838f668a80563923",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "cc3f48dd1ed7571ee0339a7f4e6d24a9a3f5804c834904668560fffbc9b1b1ae",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "dc89cf8867b62443bb4688897d946b4c1787f5d2792a442ba2e70985c4df93eb",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "393de81c2f62261573cc4045ff1b732b3b7e7d2bfa9e71303edbd2d5bc461dd3",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "7dfd087ad1586aa14d76169c4b16927effe7ff9cf793f0cb34991db83edf950d",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "abfe05895cf2518e341fc92c60a59e1eda9397db6675ff87262ebe452a603a06",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "8a1a671a755a5b800f88a6869dba09b29f9f687051303dbbe4ed7ea4a38cc615",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "8269129795cafdac88c0214d79ae99da354eaf0581958d6121297ae7e824fbfe",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "66546aac3a9a8c8a71a88c0f8e73460d957096f1887277f96efe857b8dd8f5fa",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "9167296c93526121529e05d95401b9c0a955d5b1b2663c4911af3e506db44ecb",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "8c56c53a133456381926ba69a521bf3c7fa7c33db04041700f114c80161226e7",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "d04941f924ae59636d6d86311554fd8f79d54d666b11afc0f877d310c4f7c02a",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "e77919a8592341bdd6bff320b795a2109e21800185a6be29fd738aedb38bfc84",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "bfdfc2dd2f07d68b7e5b1f60420ca953bfa1d156548a8536fb8773ce4f1388e6",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "d80c2f0826c6148a9a2dfb6d090039a0cd7daf02798014dcf8d5f379a139612f",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "12a0e2fcab2f8264f789a58e9ca9b58e491f58c134684a32c5c6af1db506baaa",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "2eb7602ef8b5c1f7d31f48786b55e23788ad26a3cda0af876a79ac53e1a126e5",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "7c1bca1273d68b4d4405d0b9423288c0f4c6748bce0d5485101d9668f0e90f04",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "2f03a420e6db1a9b1d6c1251166bd542b561331e13e9276ae02632d6dc4cfa0c",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "72620dbc6d00bca2c43334c8a7c7910e1545451dbf97bbf022d63aeda446d876",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "7846ccb2e1162f8b7f575bb202ef8824cd543750ba98074f0c3d3c20a2a2b1c1",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "ad15632f1ced9b844884ba6289c08bc30e725ca9b2656016d04a988109ee1d05",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "c5d3df2fb165d5ee7b168f2912f924cb06ae0df033e300651a8ee8a4e2792594",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "0c259429e649441b91e7784ed8ff2d2d634af4a79fc16bf4f6a2546caa616485",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "988f5b33cbc6ee9560fcf4dfedc8e3130f98f586ee24c2d7f3fc888beb3aa7df",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "2b0a722688afeaf9ed6b7e05e8cbea7368d7921deff2a742a7146774d23c6800",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "1b6b47071ba3cb9b7732fd9c8cdb6ff02d1cf7ff8200da74310815b0eeb0b30c",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "dc66edfd3963eaf08d58e55d72652ea09c7548faae314d7ba80c9c05ec38ef47",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "1df741772e611fbc6336958929bc37393f1f5588cac428ac832c26ca944c0a14",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "9fd1f8aa9835351397452188ad727334873e5fda9679291a2ba14e3748a988f8",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "6f7f1b253becdafd52dd3fa4ce2d6c1bb3ea443b804716d8c96250a5e7aaf1d9",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "6fe10529104c8e21b372c4b9be035add7ca5de17b2b17adfd79df2f0b58b87e8",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "352b0cc82630e68aad43a96fe6ca924637f1ec069d225560c2cdfe0de9dd1a32",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "7389cc02a31edb8edd4e7d117ca72186d95894549383450df9e33f04eef26625",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "20760f2150fae6d9e1902a0b4245dae0ce5c16afb3ba7c0db3ff1e829ecce19c",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "c3dd3293303f8c69dca65362ea2c579f868c8135fa5c689052287e65a012bc9d",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "7b6f0ad6c414afd3552850debf272dd6411131b26e621d60ff0047649925e578",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "df5f4d790c40e5ba2c44cb7bb27325092c5eb3f21fb21188e31611cba8fd6b3a",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "c17dc5eeebb307ae8da514ac9adbfb4ee0e3952f7262dcd53d3c8624e109b6bf",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "823be423cea868df4eb23a2426e26112f9381f1388091aea410e8325a3f80471",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "32370b97fc164be7624dce792fd9c08a9ee30ce214f899015eedfbbc0e30b3bb",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "874dfe5d4150018bfcfb2360c8c2c3420ccdcb838860b723468170aafef5407a",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "84f05dcdc5ba04346c52d2fed4e0c044e0fcc315cc2e5ba1026f98cfd0d2b72d",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "035ce50908500b546ccce7ae42ab1427a7b46f69b2cebf375ad60f21b1aab0a1",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "fe8f3f6cd5ec91583f96b8e077d63cd0c7de4e4c2afdb1f302ac9c4bc8088811",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "3ddf4860b8954f7cd9f2d33156a4825563a9564b1a145870ea865363d4920083",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "6daf98e374ec30317338330735ae1be4151a21f6264e8e390890aa7901e8e588",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "1b23c273a47da39f8a674bd65af565862f1cea2651549b2dc5c28b8ca6add3b7",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "e7c53ca794f8ac674e2d9c18879575887d595cbb38c5e7357e9eae17d23791c4",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "fb262e83559d799cab12963bc603058d8eec163a246dda512e4bd3a108b8e75c",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "483b2586fbe7c1939e068b7c2fe7237db45c797bce3be932bf86fd741a2b0a3f",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "cdd907ee2c8a89650a0b5d041b923c9ef91a8ef91d8bdbc407dae8a1d0f528bc",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "e84ddd12f798531d02048c10628402947d5587da978fc1d3d89dbdfee45b907a",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "2af59fbd5668a570f59aeedfb83fdb74be7b7d0c51408553cacd9808088bb6cf",
    CompilerShard("LLVM", v"6.0.1", Linux(:x86_64, libc=:glibc), :targz)
      => "72c383dc42eaee7d6ee64b1cd34540dcea0217ce17b7eb121b706138c515e58a",
    CompilerShard("LLVM", v"6.0.1", Linux(:x86_64, libc=:glibc), :squashfs)
      => "f246f422e7b9f921de3c01b77c838a85e6a537eb103d16d5fef8c7149a7b0942",
)


text:
using NewtonsMethod
using Test

@testset "NewtonsMethod.jl" begin
    # Write your own tests here.
end

ext: jl


text:
using Plots
using Random
Random.seed!(1)

include_implicit_dynamics()
include_ddp()
include(joinpath(pwd(),
    "examples/implicit_dynamics/models/planar_push/model_v2.jl"))

# visualize
include(joinpath(pwd(), "models/visualize.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/models/planar_push/visuals.jl"))

vis = Visualizer()
render(vis)

# Implicit dynamics
T = 26
h = 0.1
data = dynamics_data(model, h, T,
        r_func, rz_func, rθ_func, rz_array, rθ_array;
        idx_ineq = idx_ineq,
		idx_soc = idx_soc,
		z_subset_init = z_subset_init,
        diff_idx = 2,
        opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = false))

model_implicit = ImplicitDynamics{Midpoint, FixedTime}(2 * model.dim.q, model.dim.u, 0, data)

n = model_implicit.n
m = model_implicit.m

# Problem setup

# Initial and final states
q1 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, -0.01]
x1 = [q1; q1]

x_goal = 0.5
y_goal = 0.5
# x_goal = 0.25
# y_goal = 0.25
θ_goal = 0.5 * π
# qT = q1
qT = [x_goal, y_goal, θ_goal, x_goal-r_dim, y_goal-r_dim]
xT = [qT; qT]

# Objective
V = 1.0 * Diagonal([1.0, 1.0, 1.0, 0.1, 0.1])
Q_velocity = [V -V; -V V] ./ h^2.0
_Q_track = [1.0, 1.0, 1.0, 0.1, 0.1]
Q_track = 1.0 * Diagonal([_Q_track; _Q_track])

Q = [t < T ? Q_velocity + Q_track : Q_velocity + 1.0 * Q_track for t = 1:T]
q = [-2.0 * (t == T ? 1.0 : 1.0) * Q_track * xT for t = 1:T]
R = [Diagonal(1.0e-2 * ones(m)) for t = 1:T-1]
r = [zeros(m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
ul = -5.0 * ones(m)
uu = 5.0 * ones(m)
p = [t < T ? 2 * m : n - 2 * 2 for t = 1:T]
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		# c .= [ul - u; u - uu]
	else
		c .= (x - cons.con[T].info[:xT])[collect([(1:3)..., (6:8)...])]
	end
end

q0 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, -0.01]
q1 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, -0.01]
x1 = [q0; q1]
ū = [t < 5 ? [1.0; 0.0] : t < 10 ? [0.5; 0.0] : [0.0; 0.0] for t = 1:T-1]
w = [zeros(model_implicit.d) for t = 1:T-1]

x̄ = rollout(model_implicit, x1, ū, w, h, T)
q̄ = state_to_configuration(x̄)

# visualize!(vis, model, q̄, ū, Δt = h, r = r_dim)

prob = problem_data(model_implicit, obj, con_set, copy(x̄), copy(ū), w, h, T,
	analytical_dynamics_derivatives = true)

# Solve
@time stats = constrained_ddp_solve!(prob,
	max_iter = 1000,
    grad_tol = 1.0e-3,
    max_al_iter = 5,
	ρ_init = 1.0,
    ρ_scale = 10.0,
	con_tol = 0.001)

@show ilqr_iterations(stats)
@show objective(obj, prob.m_data.x̄, prob.m_data.ū)
@show prob.s_data.c_max

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/planar_push_rotate.jld2") x u
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/planar_push_rotate.jld2") x u

# compute comparable objective
J = 0.0
x = x̄
u = ū
for t = 1:T
    if t < T
        J += x[t]' * Q[t] * x[t] + q[t]' * x[t] + u[t][1:2]' * R[t][1:2, 1:2] * u[t][1:2] + r[t][1:2]' * u[t][1:2]
    elseif t == T
        J += x[t]' * Q[t] * x[t] + q[t]' * x[t]
    else
        J += 0.0
    end
end
@show J

q̄ = state_to_configuration(x)
v̄ = [(q̄[t+1] - q̄[t]) ./ h for t = 1:length(q̄)-1]

vis = Visualizer()
render(vis)
open(vis)
visualize!(vis, model,
    # q̄,
    [[q̄[1] for i = 1:15]..., q̄..., [q̄[end] for i = 1:15]...],
    Δt = h, r = r_dim, r_pusher = 0.25 * r_dim)
default_background!(vis)
settransform!(vis["/Cameras/default"],
    compose(Translation(0.0, 0.0, 50.0), LinearMap(RotZ(0.5 * pi) * RotY(-pi/2.5))))
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 50)

_create_planar_push!(vis, model,
        i = -1,
        r = r_dim, r_pusher = 0.25 * r_dim,
        box_color = RGBA(0.0, 1.0, 0.0, 0.5),
        pusher_color = RGBA(0.5, 0.5, 0.5, 0.0))

_set_planar_push!(vis, model, qT, i = -1)

# t = 1
# id = t
# tl = 0.5
# _create_planar_push!(vis, model,
#         r = r_dim,
#         r_pusher = 0.25 * r_dim,
#         tl = tl,
#         i = id)
# _set_planar_push!(vis, model, q̄[t], i = id)
#
# t = 5
# id = t
# tl = 0.6
# _create_planar_push!(vis, model,
#         r = r_dim,
#         r_pusher = 0.25 * r_dim,
#         tl = tl,
#         i = id)
# _set_planar_push!(vis, model, q̄[t], i = id)
#
# t = 10
# id = t
# tl = 0.7
# _create_planar_push!(vis, model,
#         r = r_dim,
#         r_pusher = 0.25 * r_dim,
#         tl = tl,
#         i = id)
# _set_planar_push!(vis, model, q̄[t], i = id)
#
# t = 15
# id = t
# tl = 0.8
# _create_planar_push!(vis, model,
#         r = r_dim,
#         r_pusher = 0.25 * r_dim,
#         tl = tl,
#         i = id)
# _set_planar_push!(vis, model, q̄[t], i = id)
#
# t = 20
# id = t
# tl = 0.9
# _create_planar_push!(vis, model,
#         r = r_dim,
#         r_pusher = 0.25 * r_dim,
#         tl = tl,
#         i = id)
# _set_planar_push!(vis, model, q̄[t], i = id)
#
# t = 26
# id = t
# tl = 1.0
# _create_planar_push!(vis, model,
#         r = r_dim,
#         r_pusher = 0.25 * r_dim,
#         tl = tl,
#         i = id)
# _set_planar_push!(vis, model, q̄[t], i = id)
#
# box_line_mat = LineBasicMaterial(color=color=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0), linewidth=10.0)
# pusher_line_mat = LineBasicMaterial(color=color=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0), linewidth=10.0)
#
# points_box = Vector{Point{3,Float64}}()
# points_pusher = Vector{Point{3,Float64}}()
#
# for (i, xt) in enumerate(x̄)
# 	push!(points_box, Point(xt[1], xt[2], 0.0))
#     push!(points_pusher, Point(xt[4], xt[5], 0.0))
# end
# setobject!(vis[:box_traj], MeshCat.Line(points_box, box_line_mat))
# setobject!(vis[:pusher_traj], MeshCat.Line(points_pusher, pusher_line_mat))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Test, Logging

include("../gpu_batch.jl")

"""
For an array of length `N` * 1024 with value type `T`, test that each sequential
group of 1024 elements is stably sorted.
"""
function test_multi_batch(N, T)
    dest = CUDA.zeros(Int32, N * 1024)
    host_vals = rand(T, N * 1024)
    vals = CuArray(host_vals)
    odd = collect(filter(x -> mod(x, 2) == 1, host_vals[1 + 1024 * (i - 1):1024 * i]) for i in 1:N)
    even = collect(filter(x -> mod(x, 2) == 0, host_vals[1 + 1024 * (i - 1):1024 * i]) for i in 1:N)
    @assert Set(host_vals) == union(reduce(union, map(Set, odd)), reduce(union, map(Set, even)))
    flags = map(x -> x % 2 == 0, vals)
    batch_sums = CUDA.zeros(Int32, N)

    @cuda blocks=N threads=1024 index_sort_kernel(dest, flags, batch_sums)
    @cuda blocks=N threads=1024 index_move_kernel(vals, dest)
    synchronize()
    @test Array(vals) == reduce(vcat, collect(vcat(odd[i], even[i]) for i in 1:N))
end

test_multi_batch(1, Int16)
test_multi_batch(2, Int16)
test_multi_batch(1, Int64)
test_multi_batch(2, Int64)
test_multi_batch(64, Int64)
@info "passed"


text:
import BenchPerf
import Random
using BenchmarkTools

include("utils.jl")

function increment_at!(xs, indices)
    for i in indices
        @inbounds xs[i] += 1
    end
end

function threaded_increment_at!(xss, iss, nrepeats)
    Threads.@threads for j in eachindex(xss, iss)
        for _ in 1:nrepeats
            increment_at!(xss[j], iss[j])
        end
    end
end

group = BenchmarkGroup()
for e in 5:22
    n = 2^e
    nrepeats = nrepeats_from_n(n)
    group["n=$n"] = bench = @benchmarkable(
        threaded_increment_at!(xss, iss, $nrepeats),
        setup = begin
            Random.seed!(43)
            # Iteration over nthreads can be done in a single process.  But it's
            # done via BenchmarkConfigSweeps just for the demo.
            nt = Threads.nthreads()
            n = $n
            xss = [zeros(Int, n) for _ in 1:nt]
            iss = [rand(1:n, n) for _ in 1:nt]
        end,
        # Larger `samples` (than the default: 10_000) to bound the benchmark by
        # the time limit.
        samples = 1_000_000,
    )
    # Not tuning it like ../random_increments/benchmarks.jl since `nrepeats`
    # already takes care of it. To verify this, uncomment the following line,
    # run the benchmarks, and then check `unique(df.evals)`.
    #=
    if e < 15
        tune!(bench)
    end
    =#
end

include("../quick.jl")
maybe_quick!(group)

SUITE = BenchPerf.wrap(group; detailed = 1)


text:
@symbol_func function gray_cost(cur_reactor::AbstractReactor)
  cur_cost = W_M(cur_reactor)

  cur_cost /= P_F(cur_reactor)

  cur_cost
end


text:
# Define PDF & CDF for Student's t distribution:

using SpecialFunctions
using GSL
using Optim
using DelimitedFiles

include("../../../src/loglinspace.jl")

function students_t_pdf(x,nx,sig,nu)
# PDF of the Student's t distribution.
pdf = zeros(nx)
alp = 0.5*(nu+1)
fac = gamma(alp)/sqrt(pi*nu)/gamma(0.5*nu)
for i=1:nx
  pdf[i] = fac/(1+(x[i]/sig)^2/nu)^alp
end
return pdf/sig
end

function students_t_cdf(x,nx,sig,nu)
# Computes the cumulative distribution function for Student's t distribution:
cdf = zeros(nx)
alp = 0.5*(nu+1)
fac = gamma(alp)/sqrt(pi*nu)/gamma(0.5*nu)
for i=1:nx
  y = (x[i]/sig)^2/nu
  if y < 1.0
    cdf[i] = 0.5+x[i]/sig*fac*hypergeom([0.5,alp],1.5,-y)
  else
    cdf[i] = 0.5+x[i]/sig*fac*hypergeom([0.5,1.5-alp],1.5,y/(1+y))/sqrt(1+y)
  end
end
return cdf
end

function fit_students_t(resid,s1,ndof)
# Fits a discrete distribution:
  # Fits a students-t distribution to the
  x0 = [s1,ndof]
  nresid = size(resid)[1]
  function student_model(x)
    s1,ndof = x
    println("s1: ",s1," ndof: ",ndof)
    if ndof >= 0.5 && s1 > 0.0
      pdf = students_t_pdf(resid,nresid,s1,ndof)
    else
      return 0.0
    end
    # Return a "chi-square":
    return -sum(log.(pdf))
  end
  xbest = Optim.optimize(student_model, x0)
return xbest.minimizer
end

function fit_students_t_hist(resid,binhist,s1,ndof)
# Fits a binned distribution:
  norm_hist = sum(binhist)
  # Fits a students-t distribution to the
  x0 = [s1,ndof,norm_hist]
  sigbin = copy(binhist)
  sigbin[sigbin .== 0.0] .= 1.0
  function student_model(x)
    s1,ndof,norm_hist = x
    println("s1: ",s1," ndof: ",ndof," norm: ",norm_hist)
    if ndof >= 0.0
      pdf = students_t_pdf(resid,size(resid)[1],s1,ndof)
    else
      return 0.0
    end
    # Return a "chi-square":
#    return sum((binhist[binhist .!= 0.0] .- pdf[binhist .!= 0.0].*norm_hist).^2 ./binhist[binhist .!= 0.0])
    return sum((binhist .- pdf.*norm_hist).^2 ./sigbin)
  end
  xbest = Optim.optimize(student_model, x0)
return xbest.minimizer
end

# Load in Julia wrapper of matplotlib:
using PyPlot
# Load Special Functions package for computing erf:
using SpecialFunctions

fig,axes = subplots(1,2,figsize=(10,5))

#clf()
# Read in the residuals:
dev = readdlm("../../../data/residuals_nomixture_rescale.txt")
# Plot the cumulative distribution function of the residuals:
ax = axes[1]
nobs = 447
devsort = sort(dev[:,1])
iobs = linearspace(0,1,nobs)
ax.plot(devsort,iobs,linewidth=3,label="CDF")
ax.plot(devsort[devsort .< -3],iobs[devsort .< -3],linewidth=3,color="C5")
ax.plot(devsort[devsort .>  3],iobs[devsort .>  3],linewidth=3,color="C5")
# Create an array for plotting normalized residuals:
nx = 10000
x = linearspace(-7.5,7.5,nx)
# Plot CDF of a Gaussian with \sigma = 0.9:
#p = 0.8
p = 1.0 - 0.3150298440418438
s1 = 1.0
d1 = 0.5*(1 .+erf.(x./sqrt(2)./s1))
ax.plot(x,d1,label="Single Gaussian",linewidth=1)
# Plot CDF of a Gaussian with \sigma = 2.0:
#s1 = 0.87
s1 = 0.7312227481615532
d1 = 0.5*(1 .+erf.(x./sqrt(2)./s1))
#s2 = 2.0
s2 = 1.9192905742029427
d2 = 0.5*(1 .+erf.(x./sqrt(2)./s2))
#plot(x,d2)
# Sum of CDFs:
#ax.plot(x,p*d1 .+ (1-p)*d2,label="Double Gaussian",linewidth=2,linestyle="--")
#ax.set_title("Cumulative distribution of normalized residuals of best TTV model")
# Student's t distribution - first optimize the fit to the histogram:
sig_est = 0.835
dof = 3.27
ax.plot(x,students_t_cdf(x,nx,sig_est,dof),label="Student's-t",linewidth=2,linestyle=":")
ax.set_title("Cumulative distribution")
# Label the plot:
ax.set_xlabel(L"$z = (t_{\mathrm{obs},ij}-t_{ij}(\mathbf{x}))/\sigma_{ij}$")
ax.set_ylabel(L"$P(>z)$")
ax.legend(loc = "lower right",fontsize=10)
ax.axis([-8,8,-0.05,1.05])

# Now plot a histogram:
nbin = 50; x10 = -7.5; x20 = 7.5; dx = (x20-x10)/nbin
dev_hist = zeros(nbin)
for i=1:length(dev)
  ibin = ceil(Int64,(dev[i]-x10)/dx)
  if ibin > 0 && ibin <= nbin
    dev_hist[ibin] += 1.0
  end
end
res = []
resc = []
hbin = []
for ibin=1:nbin
  push!(res,x10 + (ibin-1)*dx)
  push!(resc,x10 + (ibin-0.5)*dx)
  push!(hbin,dev_hist[ibin])
  push!(res,x10 + ibin*dx)
  push!(hbin,dev_hist[ibin])
end
ax = axes[2]
ax.semilogy(res,hbin,linestyle="None")
errs = sqrt.(dev_hist)
#errs[errs .== 0.0] .= 1.0
ax.errorbar(resc,dev_hist,yerr=errs,fmt=".",color="#1f77b4")
# Plot in a diffent color the bins with outliers:
indx = (resc .< -3.0) .| (resc .> 3.0)
ax.errorbar(resc[indx],dev_hist[indx],yerr=errs[indx],fmt=".",color="C5")

# Optimize Student's-t distribution:
ndof = 2.8
s1 = 0.79
#xbest = fit_students_t_hist(resc,dev_hist,s1,ndof)
xbest = fit_students_t(dev,s1,ndof)
sig_est = xbest[1]
dof = xbest[2]
norm_hist = size(dev)[1]

# Now overplot single and double Gaussian models:
nx = 1000
devs = linearspace(-8,8,nx)
s1 = 1.0
g1 = dx*nobs/sqrt(2*pi*s1^2).*exp.(-0.5.*(devs/s1).^2)
ax.plot(devs,g1,linewidth=2)
ax.axis([-8,8,1e-1,1e2])
#ax.set_xlabel("x = (data-model)/sigma")
ax.set_xlabel(L"$z = (t_{\mathrm{obs},ij}-t_{ij}(\mathbf{x}))/\sigma_{ij}$")
ax.set_ylabel(L"$dP/dz$")
#s1 = 0.87
s1 = 0.7312227481615532
g2 = p*dx*nobs/sqrt(2*pi*s1^2).*exp.(-0.5.*(devs/s1).^2).+
     (1-p)*dx*nobs/sqrt(2*pi*s2^2).*exp.(-0.5.*(devs/s2).^2)
#ax.plot(devs,g2,linestyle="--",linewidth=2)
ax.set_title("Histogram")

# Overplot Student's t distribution:
stpdf = students_t_pdf(devs,nx,sig_est,dof)
ax.plot(devs,stpdf*norm_hist*0.3,linestyle=":",linewidth=3)
#ax.legend(ncol=3,fontsize=8)

savefig("../Students_t_optimized.pdf",bbox_inches="tight")


text:
using Oceananigans.Solvers: calculate_pressure_source_term_fft_based_solver!

import Oceananigans.Solvers: solve_for_pressure!, source_term_storage, source_term_kernel, solution_storage

source_term_storage(solver::DistributedFFTBasedPoissonSolver) = first(solver.storage)

source_term_kernel(::DistributedFFTBasedPoissonSolver) = calculate_pressure_source_term_fft_based_solver!

solution_storage(solver::DistributedFFTBasedPoissonSolver) = first(solver.storage)


text:
module GZK_sim

using CSV
using DataFrames
using Random
using ThreadsX

include("./pions.jl")
include("./proton.jl")

using .PionFunctions
using .ProtonFunctions

export propagate

"""
    'propagate(z_dist, E_dist)'

Propagates protons to the Earth (z=0), calculating the energy they lose and the neutrinos they create along the way. This is the main loop of the simulation.

Inputs:
- z_dist = array of the initial redshift (z) of each proton
- E_dist = array of the initial energy (log10(E/eV)) of each proton

Output:
Dataframe with columns for the event #, injection energy, injection redshift, arriving proton energy, # of arriving neutrinos, and an array of the arriving neutrino energies.
"""
function propagate(z_dist::AbstractArray{T1}, E_dist::AbstractArray{T2}, seed = false) where {T1<:Real, T2<:Real}

    if seed != false
        rng = MersenneTwister(seed)
    else
        rng = MersenneTwister()
    end

    # load in the pγ cross-section data
    data = CSV.read("./data/p-gamma_cross_section.csv", DataFrame) #Source: M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
    ϵ = data[!,1]*1e3 # photon energy in MeV
    σ = data[!,2]*1e-13 # cross-section in nanometer^2

    primaries = [E_dist z_dist]

    # preallocate the output vectors. Fill it with "missing" so we can tell if something went wrong by looking at the output CSV.
    event = Vector{Union{Missing, Int}}(missing, size(primaries,1))
    injection_energy = Vector{Union{Missing, Float64}}(missing, size(primaries,1))
    redshift = Vector{Union{Missing, Float64}}(missing, size(primaries,1))
    proton_energies = Vector{Union{Missing, Float64}}(missing, size(primaries,1))
    num_of_neutrinos = Vector{Union{Missing, Float64}}(missing, size(primaries,1))
    neutrino_energies = Vector{Union{Missing, Array}}(missing, size(primaries,1))

    # Create a separate RNG for each thread
    rng_thread = map(i->MersenneTwister(rand(rng,UInt32)),1:Threads.nthreads())

    # loop over each injected primary.
    #Threads.@threads for i in (1:size(primaries,1))
    ThreadsX.foreach(1:size(primaries,1), basesize=ceil(Int, size(primaries,1)/Threads.nthreads())) do i

        mp = 938.272e6 # proton mass, eV
        Γ_initial = (big(10)^primaries[i,:][1])/mp # the initial Lorentz factor of the proton
        z = LinRange(primaries[i,:][2],0,100) # z ranges from the injection redshift to Earth (z=0)
        initial_rate = interaction_rate(Γ_initial, z[1], ϵ, σ)
        u = rand(rng_thread[Threads.threadid()]) # this random number will be used to decide when the interaction takes place
        p_initial = 1 # initial interaction probability is 0

        # define arrays to store values
        proton = [Γ_initial,initial_rate,p_initial,u]
        neutrinos = []

        # loop over each step in redshift
        for j in (2:size(z,1))

            proton[1] = adiabatic_energy_loss(proton[1], z[j-1], z[j])
            rate = interaction_rate(proton[1], z[j], ϵ, σ)
            p = survival_probability(proton[3], proton[2], z[j-1], rate, z[j])
            proton[3] = p
            if p < proton[4] # if p < u then the interaction takes place
                if rand(rng_thread[Threads.threadid()], 1:3) == 1 # neutral pions have a 1/3 chance of being produced
                    # neutral pion production creates a new proton and no neutrinos
                    proton = neutral_pion_production(proton[1], z[j], rate, rng_thread[Threads.threadid()])
                else # charged pions have a 2/3 chance of being produced
                    # charged pion production creates a new proton and 4 neutrinos
                    proton, neutrino1, neutrino2, neutrino3, neutrino4 = charged_pion_production(proton[1], z[j], rate, rng_thread[Threads.threadid()])
                    append!(neutrinos, [neutrino1, neutrino2, neutrino3, neutrino4])
                end
            else
                proton[2] = rate # if the proton doesn't interact simply update the interaction rate and loop again
            end
        end

        # loop over each neutrino that was created
        for k in (1:size(neutrinos,1))
            # neutrinos do not undergo any energy loss as they propagate except from the adiabatic expansion of the universe.
            # the energy of a neutrino arriving to Earth is therefore E/(1+z), where E and z are the energy and redshift at which the neutrino was created.
            neutrinos[k][1] = neutrinos[k][1]/(1+neutrinos[k][2])
        end

        # fill in the output vectors
        event[i] = i
        injection_energy[i] = primaries[i,:][1]
        redshift[i] = primaries[i,:][2]
        proton_energies[i] = log10(proton[1]*mp)
        num_of_neutrinos[i] = size(neutrinos,1)
        neutrino_energies[i] = log10.([convert(Float64,x[1]) for x in neutrinos].*1e6)

    end

    # construct the dataframe
    df = DataFrame(Event=event, Injection_Energy=injection_energy, Injection_Z=redshift, Arriving_Proton_Energy=proton_energies, Arriving_Neutrinos=num_of_neutrinos,
    Neutrino_Energies=neutrino_energies)

    return df

end

end


text:
module Repressilator

using ..MultiFidelityABC
using LinearAlgebra
using Dierckx
using Random

# Repressilator model specification contains:
# - Model in repressilator_model_generator - what to simulate
# - mf-abc specification (what is low-fidelity, high-fidelity, etc) in repressilator_mfabc_problem
# - Prior parameter distribution in repressilator_prior

const z =  zeros(3,3)
const nu = vcat(hcat(I,-I,z,z),hcat(z,z,I,-I))
const x0 = [0.0; 0.0; 0.0; 40.0; 20.0; 60.0]
const k_nominal = (1.0, 2.0, 5.0, 1000.0, 20.0)
const T = 10.0

function propensity!(v::Array{Float64,1}, x::Array{Float64,1}, k::NTuple{5,Float64})

    v[1:3] = k[1] .+ (k[4]*(k[5]^k[2]))./((k[5]^k[2]) .+ (x[[6,4,5]].^k[2]))
    v[4:6] = x[1:3]
    v[7:9] = k[3].*x[1:3]
    v[10:12] = k[3].*x[4:6]

    return nothing
end

function diff_propensity!(dv::Array{Float64,2}, x::Array{Float64,1}, k::NTuple{5,Float64})

    dv[1,6] = -k[4]*k[2].*(k[5].^k[2]).*(x[6].^(k[2]-1))*(k[5].^k[2] + x[6].^k[2]).^(-2);
    dv[2,4] = -k[4]*k[2].*(k[5].^k[2]).*(x[4].^(k[2]-1))*(k[5].^k[2] + x[4].^k[2]).^(-2);
    dv[3,5] = -k[4]*k[2].*(k[5].^k[2]).*(x[5].^(k[2]-1))*(k[5].^k[2] + x[5].^k[2]).^(-2);
    dv[4:6,1:3] = Matrix{Float64}(I,3,3);
    dv[7:9,1:3] = k[3]*Matrix{Float64}(I,3,3);
    dv[10:12,4:6] = k[3]*Matrix{Float64}(I,3,3);

    return nothing
end

gm = GillespieModel(nu, propensity!, x0, k_nominal, T)
tlm = TauLeapModel(nu, propensity!, x0, k_nominal, T, diff_propensity!, 0.01, 3.0, 0.01)

function summary_statistics(t::Times, x::States)
    return vcat([Spline1D(t,x[j,:];k=1)(0:T) for j in 1:size(x,1)]...)
end

Random.seed!(123)
t, x = simulate(gm) # Simulate the nominal model with a fixed seed
Random.seed!()
yo = summary_statistics(t, x)

function distance(y::Array{Float64,1})
    return norm(y-yo)/T
end

function lofi(k::Parameters)
    t, x, c_pp = simulate(tlm, k)
    return distance(summary_statistics(t,x)), c_pp
end

# The following hifi *couples* low fidelity and high fidelity simulations:
function hifi(k::Parameters, c_pp::Coarse_PP)
    t, x = complete(tlm, k, c_pp)
    return distance(summary_statistics(t,x))
end

# # This hifi version would produce independent (uncoupled) high fidelity simulations:
# function hifi(k,pass)
#     gm = GillespieModel(rep_mg, k)
#     t, x = gillespieDM(gm)
#     return summary_statistics(t,x)
# end

using Distributions
function prior()
    k2 = rand(Uniform(1,4))
    k5 = rand(Uniform(10,30))
    return k_nominal[1], k2, k_nominal[3], k_nominal[4], k5
end

mf_prob = MFABC(prior, lofi, hifi)

# Parameter functions we're interested in estimating
F = [k->Float64(1.9 < k[2] < 2.1),
     k->Float64(1.2 < k[2] < 1.4),
     k->k[2]]

end

max_stars_repo_name: tpprescott/mf-abc

max_issues_repo_name: tpprescott/mf-abc

max_forks_repo_name: tpprescott/mf-abc


text:
using ITensors, Printf, Profile, ProfileView

function main()

N = 100
sites = spinOneSites(N)
H = Heisenberg(sites)
psi0 = randomMPS(sites)
sw = Sweeps(5)
cutoff!(sw,1E-12)

#G.S. energy for N=10,S=1/2 Heisenberg is -4.258035206805

maxdim!(sw,10,20,100,100,200)
energy,psi = @time dmrg(H,psi0,sw,maxiter=3)
@printf "Final energy = %.12f\n" energy

Profile.clear()  # in case we have any previous profiling data
@profile dmrg(H,psi0,sw,maxiter=3)
ProfileView.view()

return
end


text:
abstract type ReserveDirection end
abstract type ReserveUp <: ReserveDirection end
abstract type ReserveDown <: ReserveDirection end
abstract type ReserveSymmetric <: ReserveDirection end
abstract type AbstractReserve <: Service end
abstract type Reserve{T <: ReserveDirection} <: AbstractReserve end
abstract type ReserveNonSpinning <: AbstractReserve end


text:
import Gtk

function gtkwindow(name, w, h, closecb=nothing)
    c = Gtk.Canvas()
    win = Gtk.Window(c, name, w, h)
    if closecb !== nothing
        Gtk.on_signal_destroy(closecb, win)
    end
    c
end

function display(c::Gtk.Canvas, pc::PlotContainer)
    c.draw = function(_)
        ctx = getgc(c)
        set_source_rgb(ctx, 1, 1, 1)
        paint(ctx)
        Winston.page_compose(pc, Gtk.cairo_surface(c))
    end
    Gtk.draw(c)
end


text:
type CallSpecifiedPathwiseMultiProduct{M <: MarketModelPathwiseMultiProduct, E <: ExerciseStrategy, M2 <: MarketModelPathwiseMultiProduct} <: MarketModelPathwiseMultiProduct
  underlying::M
  strategy::E
  rebate::M2
  evolution::EvolutionDescription
  isPresent::Vector{BitArray{1}}
  cashFlowTimes::Vector{Float64}
  rebateOffset::Int
  wasCalled::Bool
  dummyCashFlowsThisStep::Vector{Int}
  dummyCashFlowsGenerated::Vector{Vector{MarketModelPathWiseCashFlow}}
  currentIndex::Int
  callable::Bool
end

function CallSpecifiedPathwiseMultiProduct(underlying::MarketModelPathwiseMultiProduct,
                                  strategy::ExerciseStrategy)
  underlying = clone(underlying)
  strategy = clone(strategy)

  products = number_of_products(underlying)
  d1 = get_evolution(underlying)
  rateTimes1 = d1.rateTimes
  evolutionTimes1 = d1.evolutionTimes
  exerciseTimes = strategy.exerciseTimes

  # TODO check if rebate is there
  description = EvolutionDescription(rateTimes1, exerciseTimes)
  amounts = zeros(products, length(exerciseTimes))
  rebate = MarketModelPathwiseCashRebate(description, exerciseTimes, amounts, products)

  allEvolutionTimes = Vector{Vector{Float64}}(4)
  allEvolutionTimes[1] = evolutionTimes1
  allEvolutionTimes[2] = exerciseTimes
  allEvolutionTimes[3] = rebate.evolution.evolutionTimes
  allEvolutionTimes[4] = relevant_times(strategy)

  mergedEvolutionTimes, isPresent = merge_times(allEvolutionTimes)

  # TODO Add relevant rates
  evolution = EvolutionDescription(rateTimes1, mergedEvolutionTimes)
  cashFlowTimes = possible_cash_flow_times(underlying)
  rebateOffset = length(cashFlowTimes)
  rebateTimes = possible_cash_flow_times(rebate)
  cashFlowTimes = vcat(cashFlowTimes, rebateTimes)

  dummyCashFlowsThisStep = zeros(Int, products)
  n = max_number_of_cashflows_per_product_per_step(rebate)

  dummyCashFlowsGenerated = fill(MarketModelPathWiseCashFlow[MarketModelPathWiseCashFlow(d1.numberOfRates + 1) for i = 1:n], products)

  return CallSpecifiedPathwiseMultiProduct(underlying, strategy, rebate, evolution, isPresent, cashFlowTimes, rebateOffset, false, dummyCashFlowsThisStep,
                                  dummyCashFlowsGenerated, 1, true)
end

function next_time_step!(cs::CallSpecifiedPathwiseMultiProduct, currentState::CurveState, numberCashFlowsThisStep::Vector{Int}, genCashFlows::Vector{Vector{MarketModelPathWiseCashFlow}})
  isUnderlyingTime = cs.isPresent[1][cs.currentIndex]
  isExerciseTime = cs.isPresent[2][cs.currentIndex]
  isRebateTime = cs.isPresent[3][cs.currentIndex]
  isStrategyRelevantTime = cs.isPresent[4][cs.currentIndex]

  isDone = false

  if ~cs.wasCalled && isStrategyRelevantTime
    next_step!(cs.strategy, currentState)
  end

  if ~cs.wasCalled && isExerciseTime && cs.callable
    cs.wasCalled = get_exercise!(cs.strategy, currentState)
  end

  if cs.wasCalled
    if isRebateTime
      isDone = next_time_step!(cs.rebate, currentState, numberCashFlowsThisStep, genCashFlows)
      for i in eachindex(numberCashFlowsThisStep)
        for j = 1:numberCashFlowsThisStep[i]
          genCashFlows[i][j].timeIndex += cs.rebateOffset
        end
      end
    end
  else
    if isRebateTime
      next_time_step!(cs.rebate, currentState, cs.dummyCashFlowsThisStep, cs.dummyCashFlowsGenerated)
    end

    if isUnderlyingTime
      isDone = next_time_step!(cs.underlying, currentState, numberCashFlowsThisStep, genCashFlows)
    end
  end

  cs.currentIndex += 1
  return isDone || cs.currentIndex == length(cs.evolution.evolutionTimes) + 1
end

number_of_products(cs::CallSpecifiedPathwiseMultiProduct) = number_of_products(cs.underlying)
already_deflated(cs::CallSpecifiedPathwiseMultiProduct) = already_deflated(cs.underlying)
possible_cash_flow_times(cs::CallSpecifiedPathwiseMultiProduct) = cs.cashFlowTimes
max_number_of_cashflows_per_product_per_step(cs::CallSpecifiedPathwiseMultiProduct) =
            max(max_number_of_cashflows_per_product_per_step(cs.underlying), max_number_of_cashflows_per_product_per_step(cs.rebate))

function reset!(cs::CallSpecifiedPathwiseMultiProduct)
  reset!(cs.underlying)
  reset!(cs.rebate)
  reset!(cs.strategy)
  cs.currentIndex = 1
  cs.wasCalled = false

  return cs
end

clone(cs::CallSpecifiedPathwiseMultiProduct) = CallSpecifiedPathwiseMultiProduct(clone(cs.underlying), clone(cs.strategy), clone(cs.rebate), clone(cs.evolution), deepcopy(cs.isPresent),
                                        copy(cs.cashFlowTimes), cs.rebateOffset, cs.wasCalled, copy(cs.dummyCashFlowsThisStep), deepcopy(cs.dummyCashFlowsGenerated),
                                        cs.currentIndex, cs.callable)


text:
using TaskGraphs

## set up the environment
using GraphUtils

indicator_grid = IndicatorGrid(zeros(Bool,7,5))                           # initialize a 10x10 environment with no obstacles
indicator_grid[4,4] = true                                                                # obstacle at position 4,4
indicator_grid[1:6,2] .= true
indicator_grid[1:6,4] .= true                                                          # obstacle covering positions 6:7,6:7
env = construct_factory_env_from_indicator_grid(indicator_grid) # initialize the factory environment
@show env.vtx_map                                                                     # visualize the mapping from grid position to vertex indices

#   1   0   2   0   3
#   4   0   5   0   6
#   7   0   8   0   9
#  10   0  11   0  12
#  13   0  14   0  15
#  16   0  17   0  18
#  19  20  21  22  23

## set up the environment
## Define the initial conditions of the robots
robot_ics = [
    ROBOT_AT(1,23), # robot 1 starts at vertex 2
    ROBOT_AT(2,20), # robot 2 starts at vertex 9
]

## Define the manufacturing project
spec = ProjectSpec()
# set initial conditions of "raw materials"
set_initial_condition!(spec,OBJECT_AT(1,23)) # object 1 starts at vertex 4
set_initial_condition!(spec,OBJECT_AT(2,20))  # object 2 starts at vertex

# define the operations that need to take place

op1 = Operation(Δt=0)
op2 = Operation(Δt=0)

# inputs
set_precondition!(op1,OBJECT_AT(1,17))
set_precondition!(op2,OBJECT_AT(2,11))

# outputs
set_postcondition!(op1,OBJECT_AT(3,14))
set_postcondition!(op2,OBJECT_AT(4,8))
add_operation!(spec,op1)
add_operation!(spec,op2)

# add a terminal operation with no outputs
op3 = Operation(Δt=0)
op4 = Operation(Δt=0)

# outputs
set_precondition!(op3,OBJECT_AT(3,2))
set_precondition!(op4,OBJECT_AT(4,5))
add_operation!(spec,op3)
add_operation!(spec,op4)

# define assignment_dict
assignment_dict = Dict(
  1=>[1,3],2=>[2,4],
)
# export
#     apply_assignment_dict!

## define solver
solver = NBSSolver()
# finalize problem construction (the solver is passed as an argument here
# because it determines what cost model is used for the problem)
prob = pctapf_problem(solver,spec,env,robot_ics)
TaskGraphs.apply_assignment_dict!(get_schedule(get_env(prob)),assignment_dict,get_problem_spec(get_env(prob)))
prob = PC_MAPF(get_env(prob))
route_planner = solver.path_planner
# solve the problem
solution, cost = solve!(route_planner,prob)
println(env.vtx_map)
println(solution)
println(cost)
# check if the problem was solved to optimality
@show feasible_status(solver)
@show optimal_status(solver)


text:
using MCTS
using POMDPToolbox
using POMDPs
using Multilane
using POMCPOW
using ParallelPOMCPOW

@everywhere begin
    using POMDPs
    using ParallelPOMCPOW
    using Multilane
    POMDPs.actions(p::MLPOMDP, fsrb::ParallelPOMCPOW.FilteredSRBelief) = actions(p, fsrb.srb)
end

tii = true
wup = WeightUpdateParams(smoothing=0.0, wrong_lane_factor=0.5)

solver = POMCPOWSolver(tree_queries=1_000_000,
                               criterion=MaxUCB(10.0),
                               max_depth=20,
                               max_time=10.0,
                               enable_action_pw=false,
                               k_observation=4.0,
                               alpha_observation=1/8,
                               estimate_value=FORollout(SimpleSolver()),
                               check_repeat_obs=false,
                               node_sr_belief_updater=AggressivenessPOWFilter(wup),
                               tree_in_info=tii
                              )

psolver = ParallelPOMCPOWSolver(solver, 100)

behaviors = standard_uniform(correlation=0.75)
pp = PhysicalParam(4, lane_length=100.0)
dmodel = NoCrashIDMMOBILModel(10, pp,
                              behaviors=behaviors,
                              p_appear=1.0,
                              lane_terminate=true,
                              max_dist=500.0
                             )
rmodel = SuccessReward(lambda=2,
                       target_lane=4,
                       brake_penalty_thresh=4.0
                      )
pomdp = NoCrashPOMDP{typeof(rmodel), typeof(behaviors)}(dmodel, rmodel, 0.95, false)
mdp = NoCrashMDP{typeof(rmodel), typeof(behaviors)}(dmodel, rmodel, 0.95, false)
is = initial_state(pomdp, Base.GLOBAL_RNG)
ips = MLPhysicalState(is)

agg_up = AggressivenessUpdater(pomdp, 500, 0.1, 0.1, wup, MersenneTwister(50000))
ib = initialize_belief(agg_up, ips)

planner = solve(solver, pomdp)
pplanner = solve(psolver, pomdp)

println("Single Thread")
action(planner, ib)
@time a, info = action_info(planner, ib)
@show info[:search_time_us]/1e6
@show info[:tree_queries]

println("Parallel")
action(pplanner, ib)
@time a, info = action_info(pplanner, ib)
@show info[:search_time_us]/1e6
@show info


text:
export create_superset
function convert_to_float(value::String)

    # determining whether set is in reps or duration
    colon_loc = findfirst(":", value)
    if typeof(colon_loc) != Nothing
        return parse(Float64,value[1:colon_loc[1]-1])+
               parse(Float64,value[colon_loc[1]+1:end])/60
    else
        return parse(Float64,value)
    end
end

function convert_to_time(value::Float64)
    return string(floor(Int, value))*":"*lpad(string(floor(Int,(value%1)*60)), 2, "0")
end

function select_reps(exercise)
    difficulty = exercise[:difficulty]
    reps_0 = convert_to_float(exercise[:reps_0])
    reps_10 = convert_to_float(exercise[:reps_10])

    value = reps_0 + (reps_10 - reps_0)*difficulty/10

    # determining whether set is in reps or duration
    colon_loc = findfirst(":", exercise[:reps_0])

    if typeof(colon_loc) != Nothing
        return convert_to_time(value)
    else
        return string(floor(Int, value))
    end

end

function select_set!(superset)
    # select set
    set = copy(superset[1][1][1,:])

    # clear weights of chosen set
    superset[1][1][1, :priority] = 0
    superset[1][2] = 0

    # updating weights
    update_weights!(superset)

    # sorting
    sort_sets!(superset)

    return set
end

function update_weights!(superset)
    for i = 1:size(superset)[1]
        superset[i][2] += rand()
        for j = 1:size(superset[i][1])[1]
            superset[i][1][j,:priority] += rand()
        end
    end
end

function sort_sets!(superset)
    for i = 1:size(superset)[1]
        sort!(superset[i][1], :priority, rev=true)
    end

    sort!(superset, by=x->x[2], rev=true)
end

# Ideally, this is a heap of heaps, using BinaryHeap(Base.By(x->x[2]), data)
function create_superset(core_file, mobility_file, arms_file, legs_file)

    core_set = DataFrame(CSV.File("exercises/core.dat"))
    mobility_set = DataFrame(CSV.File("exercises/mobility.dat"))
    arms_set = DataFrame(CSV.File("exercises/arms.dat"))
    legs_set = DataFrame(CSV.File("exercises/legs.dat"))

    # Add a priority column to the data frames
    for set in [core_set, mobility_set, arms_set, legs_set]
        set.priority = rand(size(set)[1])
    end

    superset = [[core_set, rand()],
                [mobility_set, rand()],
                [arms_set, rand()],
                [legs_set, rand()]]

    return sort_sets!(superset)

end

max_stars_repo_name: leios/SwoLeios

max_issues_repo_name: leios/SwoLeios

max_forks_repo_name: leios/SwoLeios


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0202 - Miscellaneous topics

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
## Pkg.resolve()
## Pkg.instantiate() # run this if you didn't in Segment 02.01
using Random
Random.seed!(123)

# ## Plotting

# Within the many possible packages to plot in Julia we use here the `Plots` package that allows at run time to choose the plot's backend.
# The defauld backend is `gr` (that, if you want to go back to it after you choosen another backend, you can activate with `gr()`). Other common backends are the Python MatplotLib (`pyplot()`, requiring package `PyPlot`) and Plotly (`plotlyjs()` from package `PlotlyJS`).
# Actually we use `StatsPlots` that just adds a set of convenient functionalities (in `plots` terminology called "_recipes_") on top of `Plots`.

using StatsPlots # no need to `using Plots` as `Plots` elements are reexported by StatsPlots

# The basic idea is that we can draw different "graphical elements" in a Plot figure. The `plot(.)` function will first create a new figure, while `plot!(.)` will modify an existing plot (by drawing new elements on it, such as a new serie) taking the "current" plot as default if no plot object is passed as first argument.

# ### Plotting functions

# Let's start plotting a function. The FIRST TIME you invoke `plot` will take a while. This is the famous "time to first plot" problem due to the JIT compilation, but it refer only to the first plotting in a working session
plot(cos) # default [-5,+5] range
savefig("currentPlot1.svg"); #src
# ![](currentPlot.svg)
plot!(x->x^2, -2,2 ) # more explicit, with ranges
savefig("currentPlot1.svg"); #src
# ![](currentPlot.svg)
plot!(x->max(x,2), label="max function", linestyle=:dot, color=:black, title="Chart title", xlabel= "X axis", ylabel="Y axis", legend=:topleft) # a bit of design
savefig("currentPlot2.svg"); #src
# ![](currentPlot2.svg)
plot!(twinx(),x->20x,colour=RGB(20/255,120/255,13/255)) # secondary axis
savefig("currentPlot3.svg"); #src
# ![](currentPlot3.svg)

# ### Plotting data

using DataFrames
x = 11:15
data = DataFrame(a=[4,8,6,6,3], b=[2,4,5,8,6], c=[20,40,15,5,30])
plot(x, Matrix(data)) # x, series (in column)
savefig("currentPlot4.svg"); #src
# ![](currentPlot4.svg)
@df data plot(x, :a, seriestype=:bar, legend=:topleft)
savefig("currentPlot5.svg"); #src
# ![](currentPlot5.svg)
plot!(x, data.b, seriestype=:line)
savefig("currentPlot6.svg"); #src
# ![](currentPlot6.svg)
scatter!(twinx(), x, data.c) # alias for `plot!(..., seriestype=:scatter)`
savefig("currentPlot7.svg"); #src
# ![](currentPlot7.svg)

# ### Layouts with multiple plots

l = @layout [a ; b c] # a,b,c, here are just placeholders, not related with the df column names..
p1 = plot(x, data.a)
p2 = scatter(x, data.b)
p3 = plot(x, data.c)
plot(p1, p2, p3, layout = l)
savefig("currentPlot8.svg"); #src
# ![](currentPlot8.svg)

# ### Saving the plot..
savefig("myplot.png")
savefig("myplot.pdf")
savefig("myplot.svg")

# ## Probabilistic analysis

# Julia has a very elegant package to deal with probability analysis, the `Distributions` package.

using Distributions

# The idea is that we first create a `Distribution` object, of the distribution family and with the parameter required, and then we operate on it, for example to sample or to retrieve a quantile.
# The following table gives the distribution constructors for some of the most common distributions:

# | Discrete distr.   | Constructor                       |  Continuous distr. | Constructor              |
# | ----------------- | --------------------------------- | ------------------ | ------------------------ |
# | Discrete uniform  | `DiscreteUniform(lRange,uRange)`  | Uniform            | `Uniform(lRange,uRange)` |
# | Bernoulli         | `Bernoulli(p)`                    | Exponential        | `Exponential(rate)`      |
# | Binomial          | `Binomial(n,p)`                   | Laplace            | `Laplace(loc, scale)`    |
# | Categorical       | `Categorical(ps)`                 | Normal             | `Normal(μ,σ)`            |
# | Multinomial       | `Multinomial(n, ps)`              | Erlang             | `Erlang(n,rate)`         |
# | Geometric         | `Geometric(p)`                    | Cauchy             | `Cauchy(μ, σ)`           |
# | Hypergeometric    | `Hypergeometric(nS, nF, nTrials)` | Chisq              | `Chisq(df)`              |
# | Poisson           | `Poisson(rate)`                   | T Dist             | `TDist(df)`              |
# | Negative Binomial | `NegativeBinomial(nSucc,p)`       | F Dist             | `FDist(df1, df2)`        |
# |                   |                                   | Beta Dist          | `Beta(shapeα,shapeβ)`    |
# |                   |                                   | Gamma Dist         | `Gamma(shapeα,1/rateβ)`  |

d = Normal(10,3) # note that the parameter is the standard deviation, not the variance

mean(d)
var(d)
median(d)
quantile(d,0.9)
cdf(d,13.844)
pdf(d,0)
rand(d,3,4,2)

sample = rand(d,1000);
density(sample)
savefig("currentPlot9.svg"); #src
# ![](currentPlot9.svg)
plot!(d)
savefig("currentPlot10.svg"); #src
# ![](currentPlot10.svg)
fit(Normal, sample) # using MLE

# ## Curve fitting

# `LsqFit` is a very flexible data fitting package for linear and nonlinear arbitrary functions (using least squares).
# Let's use it to fit a logistic growth curve (Verhulst model) of volumes for a forest stand.
using LsqFit
data = DataFrame(
    age = 20:5:90,
    vol = [64,112,170,231,293,352,408,459,505,546,582,613,640,663,683]
) # Scots Pine, data from the UK Forestry Commission https://web.archive.org/web/20170119072737/http://forestry.gov.uk/pdf/FCBK048.pdf/$FILE/FCBK048.pdf
plot(data.vol)
savefig("currentPlot11.svg"); #src
# ![](currentPlot11.svg)

logisticModel(age,parameters) = parameters[1]/(1+exp(-parameters[2] * (age-parameters[3]) ))
logisticModelVec(age,parameters) = logisticModel.(age,Ref(parameters))

initialParameters = [1000,0.02,50] #max growth; growth rate, mid age
fitobject         = curve_fit(logisticModelVec, data.age, data.vol, initialParameters)
fitparams         = fitobject.param

fitobject.resid
residuals = logisticModelVec(data.age,fitparams) .- data.vol
fitobject.resid == residuals

sigma            = stderror(fitobject)
confidence_inter = confidence_interval(fitobject, 0.05) # 5% significance level

x = 0:maximum(data.age)*1.5
plot(x->logisticModel(x,fitparams),0,maximum(x), label= "Fitted vols", legend=:topleft)
plot!(data.age, data.vol, seriestype=:scatter, label = "Obs vols")
plot!(data.age, residuals, seriestype=:bar, label = "Residuals")
savefig("currentPlot12.svg"); #src
# ![](currentPlot12.svg)

# ## Constrained optimisation

# `JuMP` is the leading library to express complex optimisation problems in a clear, mathematical friendly syntax, compute the information required by the solver engines to solve the optimisation problem, pass the problem to the aforementioned solver engines and retrieve the solutions.

# JuMP has the same flexibility and expressivity of dedicated _algeabric modelling languages_ as `GAMS` or `AMPL` but with the advantage of being a library within a much more general programming language, with larger community, development tools, language constructs and possibility to interface the specific "optimisation component" of a model with the rest of the model.

# We will see how to specify the variables, the constraints and the objective function of an optimisation model, how to "solve" it and how to retrieve the optimal values

# ### A linear example: the classical "transport" problem

# Obj: minimise transport costs $c$ from several plants $p$ to several markets $m$ under the contraint to satisfy the demand $d_m$ at each market while respecting the production capacity $c_p$ of each plant:
# $min_{x_{p,m}} \sum_p \sum_m c_{p,m} * x_{p,m}$
# subject to:
# $\sum_m x_{p,m} \leq d_m$
# $\sum_p x_{p,m} \geq c_m$

using  JuMP, GLPK, DataFrames, CSV

# #### "Sets" and exogenous parameters definition

# index sets
orig = ["Epinal", "Bordeaux", "Grenoble"] # plant/sawmills origin of timber
dest = ["Paris", "Lyon", "Nantes", "Tulouse", "Lille", "Marseille", "Strasbourg"] # markets
prod = ["Fuelwood", "Sawnwood", "Pannels"]

# Read input data into DataFrames and then into Dictionaries with 2D or 3D tuples of index sets as keys

# supply(prod, orig) amounts available at origins
supplytable = CSV.read(IOBuffer("""
prod     Epinal Bordeaux Grenoble
Fuelwood 400    700      800
Sawnwood 800    1600     1800
Pannels  200    300      300
"""), DataFrame, delim=" ", ignorerepeated=true,copycols=true)
supply = Dict( (r[:prod],o) => r[Symbol(o)] for r in eachrow(supplytable), o in orig)

# demand(prod, dest) amounts required at destinations
demandtable = CSV.read(IOBuffer("""
prod      Paris Lyon Nantes Tulouse Lille Marseille Strasbourg
Fuelwood  300   300  100    75      650   225       250
Sawnwood  500   750  400    250     950   850       500
Pannels   100   100  0      50      200   100       250
"""), DataFrame, delim=" ", ignorerepeated=true,copycols=true)
demand = Dict( (r[:prod],d) => r[Symbol(d)] for r in eachrow(demandtable), d in dest)

# limit(orig, dest) of total units from any origin to destination
defaultlimit = 625.0
limit = Dict((o,d) => defaultlimit for o in orig, d in dest)

# cost(prod, orig, dest) Shipment cost per unit
costtable = CSV.read(IOBuffer("""
prod     orig     Paris Lyon Nantes Tulouse Lille Marseille Strasbourg
Fuelwood Epinal   30    10   8      10      11    71        6
Fuelwood Bordeaux 22    7    10     7       21    82        13
Fuelwood Grenoble 19    11   12     10      25    83        15

Sawnwood Epinal   39    14   11     14      16    82        8
Sawnwood Bordeaux 27    9    12     9       26    95        17
Sawnwood Grenoble 24    14   17     13      28    99        20

Pannels  Epinal   41    15   12     16      17    86        8
Pannels  Bordeaux 29    9    13     9       28    99        18
Pannels  Grenoble 26    14   17     13      31    104       20
"""), DataFrame, delim=" ", ignorerepeated=true,copycols=true)
cost = Dict( (r[:prod],r[:orig],d) => r[Symbol(d)] for r in eachrow(costtable), d in dest)

# #### Optimisation model definition

trmodel = Model(GLPK.Optimizer)
set_optimizer_attribute(trmodel, "msg_lev", GLPK.GLP_MSG_ON)

# #### Model's endogenous variables definition

@variables trmodel begin
    x[p in prod, o in orig, d in dest] >= 0
end;

# #### Constraints definition

@constraints trmodel begin
    supply[p in prod, o in orig], # observe supply limit at plant/sawmill origin o
        sum(x[p,o,d] for d in dest) <= supply[p,o]
    demand[p in prod, d in dest], # satisfy demand at market dest d
        sum(x[p,o,d] for o in orig) >= demand[p,d]
    c_total_shipment[o in orig, d in dest],
        sum(x[p,o,d] for p in prod) <= limit[o,d]
end;

# #### Objective function definition
@objective trmodel Min begin
    sum(cost[p,o,d] * x[p,o,d] for p in prod, o in orig, d in dest)
end

# #### Human-readable visualisation of the model

print(trmodel)

# #### Model resolution
optimize!(trmodel)
status = termination_status(trmodel)

# #### Post-resolution information retrieval

# Here, after the model has been "solved", we can retrieve information as the optimal level of the endogenous variables, the value of the opjective function at these optimal levels and the shadow costs of the contraints.

if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT) && has_values(trmodel)
    println("#################################################################")
    if (status == MOI.OPTIMAL)
        println("** Problem solved correctly **")
    else
        println("** Problem returned a (possibly suboptimal) solution **")
    end
    println("- Objective value (total costs): ", objective_value(trmodel))
    println("- Optimal routes:\n")
    optRoutes = value.(x)
    for p in prod
      println("\n* $(p):")
      [println("$o --> $d: $(optRoutes[p,o,d])") for o in orig, d in dest]
      println("- Shadow prices of supply:")
      [println("$o = $(dual(supply[p,o]))") for o in orig]
      println("- Shadow prices of demand:")
      [println("$d = $(dual(demand[p,d]))") for d in dest]
    end
else
    println("The model was not solved correctly.")
    println(status)
end

# ### A nonlinear example: portfolio optimisation

# The problem objective is to choose the shares of different assets in the portfolio (here forest species, but the example is exactly the same considering other assets, for example financial investments) that maximise the portfolio expected returns while minimising its expected variance under the portfolio owner risk aversion risk.
# Here the "returns" are based on the timber production and the covariance between individual species of the portfolio is based on the observed volume growth covariances.
# The idea is that within the infinite possible allocations, the locus of those allocations for which is not possible to increase the portfolio profitability without increasing also its variance and the converse whose variance can not be lowered without at the same time lower its expected profitability are efficient in the Pareto meaning and form an "efficient frontier". Within this frontier the problem is to find the unique point that maximise the utility of the portfolio's owner given its risk aversion characteristic.
# Graphically the problem is depicted i nthe following picture:

# ![The efficient frontier and the owner utility curves](https://raw.githubusercontent.com/sylvaticus/IntroSPMLJuliaCourse/main/lessonsSources/02_-_JULIA2_-_Scientific_programming_with_Julia/graph_eff_frontier_v2.png)

# Data originally from the Institut national de l'information géographique et forestière (IGN) of France. See the paper [A. Dragicevic, A. Lobianco, A. Leblois (2016), ”Forest planning and productivity-risk trade-off through the Markowitz mean-variance model“, Forest Policy and Economics, Volume 64](http://dx.doi.org/10.1016/j.forpol.2015.12.010) for a thorough discussion of this model.

# Declare the packages we are going to use:
using JuMP, Ipopt, StatsPlots

# Forest species names
species   = ["Chêne pédonculé", "Chêne sessile", "Hêtre", "Pin sylvestre"]
nSpecies  = length(species)
# Average productiities by specie
# This is implemented in a dictionary: key->value
y   = Dict( "Chêne pédonculé" => 1.83933333333333,
            "Chêne sessile"   => 2.198,
            "Hêtre"           => 3.286,
            "Pin sylvestre"   => 3.3695)

# Covariance matrix between species
σtable = [[0.037502535947712	0.016082745098039	0.027797176470588	-0.025589882352942]
          [0.016082745098039	0.015177019607843	0.018791960784314	-0.102880470588234]
          [0.027797176470588	0.018791960784314	0.031732078431373	-0.166391058823529]
          [-0.025589882352942	-0.102880470588234	-0.166391058823529	2.02950454411765]]
# We reshape the covariance matrix in a dictionary (sp1,sp2) -> value
# The function (ix,x) = enumerate(X) returns a tuple of index position and element
# for each element of an array
σ = Dict((i,j) => σtable[i_ix,j_ix] for (i_ix,i) in enumerate(species), (j_ix,j) in enumerate(species))

################################################################################
###### Showing the possible mean/variance of the portfolio by simulation #######
################################################################################

nSamples = 1000
shares   = rand(nSamples,nSpecies);
# Converting to probabilities:
import BetaML.Utils:softmax
[shares[i,:]  = softmax(shares[i,:], β=one.(shares[i,:]) .* 5) for i in 1:nSamples]

pScores = Array{Float64,2}(undef,0,2)
for i in 1:nSamples
    global pScores
    pVar    = sum(shares[i,j1] * shares[i,j2] * σ[species[j1],species[j2]] for j1 in 1:nSpecies, j2 in 1:nSpecies)
    pY      = sum(shares[i,j]*y[species[j]] for j in 1:nSpecies)
    pScores = vcat(pScores,[pVar pY])
end

scatter(pScores[:,1],pScores[:,2],colour=:blue)
savefig("currentPlot13.svg"); #src
# ![](currentPlot13.svg)

################################################################################
### Finding (one) optimal portfolio ############################################
################################################################################

# Risk aversion coefficient
α = 0.1

# We declare an optimisation problem, we name it "m" and we let JuMP associate it with the most
# suitable solver within the one installed:
port = Model(Ipopt.Optimizer)

# We declare a set of variables, indicized by the species name:
@variables port begin
    x[i in species] >= 0
end

# We declare the constraint shat the sum of shares must be equal to 1
@constraint(port, c_share, sum(x[i] for i in species) == 1)

#=
@objective port Min begin
  α *  sum(x[i] * x[j] * σ[i,j] for i in species for j in species) - sum(x[i] * y[i] for i in species)
end
=#

@NLobjective port Min α *  sum(x[i] * x[j] * σ[i,j] for i in species for j in species) - sum(x[i] * y[i] for i in species)

# Print the optimisation model in nice human-readable format:
print(port)

# Solve the model and return the solving status:
optimize!(port)
status = termination_status(port)

# Return the objective:
println("Objective value: ", objective_value(port))
# Return the value of the decision variable:
optShares = value.(x)
[println("$sp = $(optShares[sp])") for sp in species];
pOptVar = sum(optShares[species[j1]] * optShares[species[j2]] * σ[species[j1],species[j2]] for j1 in 1:nSpecies, j2 in 1:nSpecies)
pOptY   = sum(optShares[species[j]]*y[species[j]] for j in 1:nSpecies)

function computeOptimalPortfolio(species,y,σ,α)
    port = Model(Ipopt.Optimizer)
    set_optimizer_attribute(port, "print_level", 0)
    @variables port begin
        x[i in species] >= 0
    end
    @constraint(port, c_share, sum(x[i] for i in species) == 1)
    @NLobjective port Min α *  sum(x[i] * x[j] * σ[i,j] for i in species for j in species) - sum(x[i] * y[i] for i in species)
    optimize!(port)
    status = termination_status(port)
    optShares = value.(x)
    pOptVar = sum(optShares[species[j1]] * optShares[species[j2]] * σ[species[j1],species[j2]] for j1 in 1:nSpecies, j2 in 1:nSpecies)
    pOptY   = sum(optShares[species[j]]*y[species[j]] for j in 1:nSpecies)
    return (pOptVar,pOptY)
end

αs = [1000,100,10,1,0.1,0.05,0.02,0.01]
pOptScores = Array{Float64,2}(undef,0,2)
for α in αs
    global pOptScores
    pVar,pY =computeOptimalPortfolio(species,y,σ,α)
    pOptScores = vcat(pOptScores,[pVar pY])
end
scatter!(pOptScores[:,1],pOptScores[:,2],colour=:red)
savefig("currentPlot14.svg"); #src
# ![](currentPlot14.svg)

αs = [82.45,50,30,20,15,12,10,9,8,7,6,5]
pOptScores = Array{Float64,2}(undef,0,2)
for α in αs
    global pOptScores
    pVar,pY =computeOptimalPortfolio(species,y,σ,α)
    pOptScores = vcat(pOptScores,[pVar pY])
end

scatter(pOptScores[:,1],pOptScores[:,2],colour=:red)
savefig("currentPlot15.svg"); #src
# ![](currentPlot15.svg)

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
using ChainRulesCore
using Documenter
using DocThemeIndigo
using Markdown

DocMeta.setdocmeta!(
    ChainRulesCore,
    :DocTestSetup,
    quote
        using Random
        Random.seed!(0)  # frule doctest shows output

        using ChainRulesCore, LinearAlgebra
        # These rules are all actually defined in ChainRules.jl, but we redefine them here to
        # avoid the dependency.
        @scalar_rule(sin(x), cos(x))  # frule and rrule doctest
        @scalar_rule(sincos(x), @setup((sinx, cosx) = Ω), cosx, -sinx)  # frule doctest
        @scalar_rule(hypot(x::Real, y::Real), (x / Ω, y / Ω))  # rrule doctest
    end,
)

indigo = DocThemeIndigo.install(ChainRulesCore)

makedocs(;
    modules=[ChainRulesCore],
    format=Documenter.HTML(;
        prettyurls=false,
        assets=[indigo],
        mathengine=MathJax3(
            Dict(
                :tex => Dict(
                    "inlineMath" => [["\$", "\$"], ["\\(", "\\)"]],
                    "tags" => "ams",
                    # TODO: remove when using physics package
                    "macros" => Dict(
                        "ip" => ["{\\left\\langle #1, #2 \\right\\rangle}", 2],
                        "Re" => "{\\operatorname{Re}}",
                        "Im" => "{\\operatorname{Im}}",
                        "tr" => "{\\operatorname{tr}}",
                    ),
                ),
            ),
        ),
    ),
    sitename="ChainRules",
    authors="Jarrett Revels and other contributors",
    pages=[
        "Introduction" => "index.md",
        "FAQ" => "FAQ.md",
        "Rule configurations and calling back into AD" => "config.md",
        "Opting out of rules" => "opting_out_of_rules.md",
        "Writing Good Rules" => "writing_good_rules.md",
        "Complex Numbers" => "complex.md",
        "Deriving Array Rules" => "arrays.md",
        "Debug Mode" => "debug_mode.md",
        "Gradient Accumulation" => "gradient_accumulation.md",
        "Usage in AD" => "use_in_ad_system.md",
        "Converting ZygoteRules" => "converting_zygoterules.md",
        "Tips for making packages work with AD" => "tips_for_packages.md",
        "Design" => [
            "Changing the Primal" => "design/changing_the_primal.md",
            "Many Differential Types" => "design/many_differentials.md",
        ],
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/JuliaDiff/ChainRulesCore.jl.git", devbranch="main", push_preview=true
)


text:
using GaussQuadrature
using PyPlot

x = collect(linspace(0, 1, 201))
const n = 5

a, b = shifted_legendre_coefs(Float64, n)
p = orthonormal_poly(x, a, b)
grid(true)

figure(1)
plot(x, p)
const r = 0
α, β = logweight_coefs(Float64, n+1, r)
q = orthonormal_poly(x, α[1:n], β)
title("Orthonormal shifted Legendre polynomials")

figure(2)
plot(x, q)
grid(true)
s = "Orthonormal polymials with weight \$\\log x^{-1}\$"
title(latexstring(s))


text:
export readSrcRcvLocationFile,generateSrcRcvProjOperators,readDataFileToDataMat,writeSrcRcvLocFile,writeDataFile,limitDataToOffset

function readSrcRcvLocationFile(filename,M::RegularMesh)
# read files for sources and receivers location in UTM and translate them into indices locations.
# file contains a table: idx xUTM yUTM zUTM
A = readdlm(filename);
domainBoundaryUTM = M.domain;
n = M.n;
# domainBoundaryUTM: [x_start, x_end, y_start, y_end, z_start, z_end];
# n: number of cells. (n-1: number of cells)

for k=1:M.dim
	A[:,k+1] = A[:,k+1] - domainBoundaryUTM[2*k-1];
	A[:,k+1] = A[:,k+1] ./ (domainBoundaryUTM[2*k] - domainBoundaryUTM[2*k-1]);
	A[:,k+1] = round.(A[:,k+1] .* (n[k]))+1.0;
end
A = round.(Int,A);
return A;
end

function generateSrcRcvProjOperators(n_nodes::Vector,SrcRcvNodeLoc::Array{Int64})
# gets a sources/receivers location table (in indices) and generates operators.
nSrcRcv = size(SrcRcvNodeLoc,1);
I = zeros(Int64,nSrcRcv);
J = zeros(Int64,nSrcRcv);
V = zeros(Float64,nSrcRcv);
for k=1:nSrcRcv
	J[k] = k;
	if length(n_nodes)==2
		I[k] = sub2ind(tuple(n_nodes...),SrcRcvNodeLoc[k,2],SrcRcvNodeLoc[k,3]);
	else
		I[k] = sub2ind(tuple(n_nodes...),SrcRcvNodeLoc[k,2],SrcRcvNodeLoc[k,3],SrcRcvNodeLoc[k,4]);
	end
	V[k] = 1.0;
end
P = sparse(I,J,V,prod(n_nodes),nSrcRcv);
end

function readDataFileToDataMat(filename,SrcNodeLoc::Array{Int64},RcvNodeLoc::Array{Int64})
# reads a data file into a data matrix D for the code to work with
Data = readdlm(filename);
# Data is a table of #src #rcv tau(src,rcv) Wd(src,rcv) and the rest is zero.
# CellLocMap is a table of #src i j k
numSrc = size(SrcNodeLoc,1);
numRcv = size(RcvNodeLoc,1);

if size(Data,2) == 4
	D = zeros(Float64,numRcv,numSrc);
	Wd = zeros(Float64,numRcv,numSrc);
	SrcIdxs = SrcNodeLoc[:,1];
	RcvIdxs = RcvNodeLoc[:,1];
	for k=1:size(Data,1)
		srcIdx = convert(Int64,round(Data[k,1]));
		rcvIdx = convert(Int64,round(Data[k,2]));
		isrc = binarySearch(SrcIdxs,srcIdx);
		ircv = binarySearch(RcvIdxs,rcvIdx);
		D[ircv,isrc] = Data[k,3] ;#+ 1im*Data[k,4];
		Wd[ircv,isrc] = Data[k,4] ;#+ 1im*Data[k,5];
	end
elseif size(Data,2) == 6
	D = zeros(Complex128,numRcv,numSrc);
	Wd = zeros(Complex128,numRcv,numSrc);
	SrcIdxs = SrcNodeLoc[:,1];
	RcvIdxs = RcvNodeLoc[:,1];
	for k=1:size(Data,1)
		srcIdx = convert(Int64,round(Data[k,1]));
		rcvIdx = convert(Int64,round(Data[k,2]));
		isrc = binarySearch(SrcIdxs,srcIdx);
		ircv = binarySearch(RcvIdxs,rcvIdx);
		D[ircv,isrc] = Data[k,3] + 1im*Data[k,5];
		Wd[ircv,isrc] = Data[k,4] + 1im*Data[k,6];
	end

end
return D,Wd;
end

function binarySearch(arr, value)
    low = 1
    high = length(arr)
    while low <= high
        mid = round(Int,(low+high)/2)
        if arr[mid] > value
          high = mid-1
        elseif arr[mid] < value
          low = mid+1
        else
          return mid
        end
    end
    return -1
end

function writeSrcRcvLocFile(filename,Msh::RegularMesh,pad::Int64,jump::Int64)
# writes an equally distanced sources/receivers file.
n_tup = tuple((Msh.n+1)...); # n_tup is nodes...
Q = zeros(n_tup);

if Msh.dim==2
	Q[pad+1:jump:end-pad,1] = 1.0
else
	Q[pad+1:jump:end-pad,pad+1:jump:end-pad,1] = 1.0;
end
numSrcRcv = round(Int,sum(Q));
Q = find(Q);
SrcRcvTable = zeros(numSrcRcv,1+Msh.dim);
SrcRcvIDs = randperm(numSrcRcv);
SrcRcvIDs = sort(SrcRcvIDs[1:numSrcRcv]);
SrcRcvTable[:,1] = SrcRcvIDs;
startLocations = Msh.domain[1:2:end];
for k=1:numSrcRcv
	SrcRcvTable[k,2:end] = (collect(ind2sub(n_tup,Q[k]))-1).*Msh.h + startLocations;
end
writedlm(filename,SrcRcvTable);
end

function writeDataFile(filename,D::Union{Array{Float64,2},Array{Complex128,2}},Wd::Union{Array{Float64,2},Array{Complex128,2}},SrcNodeLoc::Array{Int64},RcvNodeLoc::Array{Int64})
srcIDs = SrcNodeLoc[:,1];
rcvIDs = RcvNodeLoc[:,1];

if eltype(D)==Float64
	println("Writing real data file");
	Data = zeros(Float32,countnz(Wd),4);
	k = 1;
	for j = 1:size(D,2)
		for i = 1:size(D,1)
			if Wd[i,j] != 0.0
				Data[k,1] = srcIDs[j];
				Data[k,2] = rcvIDs[i];
				Data[k,3] = D[i,j];
				Data[k,4] = Wd[i,j];
				k+=1;
			end
		end
	end
else
	println("Writing ``complex'' data file");
	Data = zeros(Float32,countnz(Wd),6);
	k = 1;
	for j = 1:size(D,2)
		for i = 1:size(D,1)
			if Wd[i,j] != 0.0
				Data[k,1] = srcIDs[j];
				Data[k,2] = rcvIDs[i];
				Data[k,3] = real(D[i,j]);
				Data[k,4] = real(Wd[i,j]);
				Data[k,5] = imag(D[i,j]);
				Data[k,6] = imag(Wd[i,j]);
				k+=1;
			end
		end
	end
end
writedlm(filename,Data);
return;
end

function limitDataToOffset(Wd::Union{Array{Float64,2},Array{Complex128,2}},SrcNodeLoc::Array{Int64},RcvNodeLoc::Array{Int64},offsetInNodes::Int64)
for j = 1:size(Wd,2)
	srcLoc = SrcNodeLoc[j,2:end];
	for i = 1:size(Wd,1)
		rcvLoc = RcvNodeLoc[i,2:end];
		dist = norm(srcLoc - rcvLoc);
		# dist = rcvLoc[1] - srcLoc[1];
		if dist > offsetInNodes || dist < 2
			Wd[i,j] = 0.0;
		end
	end
end
return Wd;
end


text:
# functions to undo incycle, containRoot, gammaz
# originally in functions.jl
# Claudia March 2015

# ---------------------------------------- undo update of new hybridization --------------------------------

# function to undo updateInCycle which returns an array
# of edges/nodes changed
function undoInCycle!(edges::Array{Edge,1},nodes::Array{Node,1})
    for e in edges
        e.inCycle = -1;
    end
    for n in nodes
        n.inCycle = -1;
    end
end

# function to undo updateContainRoot (which returns an array
# of edges changed) and gives value bool, by default true
function undoContainRoot!(edges::Array{Edge,1}, bool::Bool)
    for e in edges
        e.containRoot = bool
    end
end

undoContainRoot!(edges::Vector{Edge}) = undoContainRoot!(edges,true)

# function to undo updateGammaz which returns an array
# of edges changed
# it only changes the status of istIdentifiable to true
function undoistIdentifiable!(edges::Array{Edge,1})
    for e in edges
        !e.istIdentifiable ? e.istIdentifiable = true : e.istIdentifiable = false;
    end
end

"""
    undoGammaz!(node, network)

Undo `updateGammaz!` for the 2 cases: bad diamond I,II.
`node` should be a hybrid node.
Set length to edges that were not identifiable and
change edges' `gammaz` attribute to -1.0.
Recalculate branch lengths in terms of `gammaz`.
*warning*: needs to know `incycle` attributes
"""
function undoGammaz!(node::Node, net::HybridNetwork)
    node.hybrid || error("cannot undo gammaz if starting node is not hybrid")
    if(node.isBadDiamondI)
        edge_maj, edge_min, tree_edge2 = hybridEdges(node);
        other_maj = getOtherNode(edge_maj,node);
        other_min = getOtherNode(edge_min,node);
        edgebla,tree_edge_incycle1,tree_edge = hybridEdges(other_min);
        edgebla,tree_edge_incycle2,tree_edge = hybridEdges(other_maj);
        other_min.gammaz != -1 || error("bad diamond I in node $(node.number) but no gammaz updated correctly")
        setLength!(tree_edge_incycle1,-log(1-other_min.gammaz))
        other_maj.gammaz != -1 || error("bad diamond I in node $(node.number) but no gammaz updated correctly")
        setLength!(tree_edge_incycle2,-log(1-other_maj.gammaz))
        if approxEq(other_maj.gammaz,0.0) && approxEq(other_min.gammaz,0.0)
            setGamma!(edge_maj,0.0, true) # gamma could be anything if both gammaz are 0.0, but will set to 0.0
            setLength!(edge_maj,0.0)
            setLength!(edge_min,0.0)
        else
            setGamma!(edge_maj,other_maj.gammaz / (other_maj.gammaz+other_min.gammaz), true)
        end
        other_min.gammaz = -1.0
        other_maj.gammaz = -1.0
        tree_edge_incycle1.istIdentifiable = true;
        tree_edge_incycle2.istIdentifiable = true;
        edge_maj.istIdentifiable = true;
        edge_min.istIdentifiable = true;
        node.isBadDiamondI = false
        net.numBad -= 1
    elseif(node.isBadDiamondII)
        edge_maj, edge_min, tree_edge2 = hybridEdges(node);
        tree_edge2.istIdentifiable = true
        node.isBadDiamondII = false
    elseif(node.isBadTriangle)
        edge_maj, edge_min, tree_edge2 = hybridEdges(node);
        tree_edge2.istIdentifiable = true
        node.isBadTriangle = false
    elseif(node.isVeryBadTriangle || node.isExtBadTriangle)
        node.isVeryBadTriangle = false
        node.isExtBadTriangle = false
        net.hasVeryBadTriangle = false
    else
        edge_maj, edge_min, tree_edge2 = hybridEdges(node);
        edge_maj.istIdentifiable = isEdgeIdentifiable(edge_maj)
        edge_min.istIdentifiable = isEdgeIdentifiable(edge_min)
    end
end


text:
"""
    Microphysics_0M

Zero-moment bulk microphysics scheme that instantly removes
moisture above certain threshold.
This is equivalent to instanteneous conversion of cloud condensate
into precipitation and precipitation fallout with infinite
terminal velocity.

"""
module Microphysics_0M

using ..Thermodynamics

using CLIMAParameters
using CLIMAParameters.Atmos.Microphysics_0M

const APS = AbstractParameterSet

export remove_precipitation

"""
    remove_precipitation(param_set::APS, q; q_vap_sat)

 - `param_set` - abstract parameter set
 - `q` - current PhasePartition
 - `q_vap_sat` - water vapor specific humidity at saturation

Returns the `q_tot` tendency due to the removal of precipitation.
The tendency is obtained assuming a relaxation with a constant timescale
to a state with precipitable water removed.
The threshold for when to remove `q_tot` is defined either by the
condensate specific humidity or supersaturation.
The thresholds and the relaxation timescale are defined in
CLIMAParameters.
"""
function remove_precipitation(
    param_set::APS,
    q::PhasePartition{FT},
) where {FT <: Real}

    _τ_precip::FT = τ_precip(param_set)
    _qc_0::FT = qc_0(param_set)

    return -max(0, (q.liq + q.ice - _qc_0)) / _τ_precip
end
function remove_precipitation(
    param_set::APS,
    q::PhasePartition{FT},
    q_vap_sat::FT,
) where {FT <: Real}

    _τ_precip::FT = τ_precip(param_set)
    _S_0::FT = S_0(param_set)

    return -max(0, (q.liq + q.ice - _S_0 * q_vap_sat)) / _τ_precip
end

end #module Microphysics_0M.jl


text:
#Test code
#-------------------------------------------------------------------------------

using SpiceData

@warn("No real test code yet... just ensuring that \"using\" works.")

@info("See sample/demo*.jl for sample usage.")

:Test_Complete


text:
# Interactive Expectation Maximization demo for the 'old faithful' data set
# This script should be executed from the commandline

using DataFrames, CSV
include("gmm_plot.jl")

function continueOrExit()
    # Helper function: wait for keypress, quit on q
    if (readline() == "q") exit() end
end

old_faithful = CSV.read("../datasets/old_faithful.csv")
X = convert(Matrix{Float64}, [old_faithful[1] old_faithful[2]]') # Every column is a data point
N = size(X, 2)

# Initialization
#clusters = MvNormal[MvNormal([3.;50.].*rand(2)+[2.;40.], diagm([.5;10^2])); MvNormal([3.;50.].*rand(2)+[2.;40.], diagm([.5;10^2]))]
clusters = [MvNormal([4.;60.], diagm([.5;10^2]));
            MvNormal([2.;80.], diagm([.5;10^2]))]
π_hat = [0.5;0.5]
γ = fill!(Matrix{Float64}(2,N), NaN)

ion()
plotGMM(X, clusters, γ)
print("Press [enter] for E-step or [q] to quit: "); continueOrExit()

while(true)
    # E-step: update γ
    norm = [pdf(clusters[1], X) pdf(clusters[2], X)] * π_hat
    γ[1,:] = (π_hat[1] * pdf(clusters[1],X) ./ norm)'
    γ[2,:] = 1 - γ[1,:]
    clf(); plotGMM(X, clusters, γ)
    print("Press [enter] for M-step or [q] to quit: "); continueOrExit()

    # M-step: update cluster parameters and π_hat
    m = sum(γ, 2)
    π_hat = m / N
    μ_hat = (X * γ') ./ m'
    for k=1:2
        μ_k = μ_hat[:,k]
        Z = (X .- μ_k)
        Σ_k = Hermitian(((Z .* (γ[k,:])') * Z') / m[k])
        clusters[k] = MvNormal(μ_k, convert(Matrix, Σ_k))
    end
    clf(); plotGMM(X, clusters, γ)
    print("Press [enter] for E-step or [q] to quit: "); continueOrExit()
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
module LangrethsRules

using LinearAlgebra

export TimeOrderedGreenFunction, TimeOrderedConvolution
export greater, lesser

include("common.jl")
include("rules.jl")

end # module


text:
## This wrapper for json allows us to override certain types
ec_json(x::Any) = json(x)

#function ec_json(x::DateTime) =  "new Date($(year(x)), $(month(x)), $(day(x)), $(hour(x)), $(minute(x)), $(second(x)), 0)"
#    ## Date(year, month, day, hours, minutes, seconds, milliseconds)

## Open a url using our heuristic
function open_url(url::String)
    @osx_only     run(`open $url`)
    @linux_only   run(`xdg-open $url`)
#    @windows_only run(`start $url`)
end

#

echarts_tpl = """
<meta charset="utf-8">

"""

## Render charts
## io -- render to io stream
## fname -- render to file
## none -- create html file, show in browser

function render(io::IO, opt::Dict; online::Bool = false)

    js_url = online ? js_online : js_local

    d = Dict(
        :id => ECharts.set_id(),
        :chart_options => json(opt),
        :js_url => js_url
        )

    out = Mustache.render(ECharts.echarts_tpl, d)

    print(io, out)
end

function render(fname::String, opt::Dict; online::Bool = false)
    io = open(fname, "w")
    render(io, opt; online = online)
    close(io)
    open_url(fname)
end

function render(opt::Dict; online::Bool = false)
    @windows_only   fname = "EC_" * randstring(10) * ".html"
    @unix_only      fname = tempname() * ".html"
    render(fname, opt; online = online)
end