From db98ecd1c21150a1a76d9ac555b70299de69d9bc Mon Sep 17 00:00:00 2001 From: pjentsch <pjentsch@uwaterloo.ca> Date: Thu, 8 Apr 2021 01:22:34 -0400 Subject: [PATCH] finally contact network done --- CovidAlertVaccinationModel/Manifest.toml | 183 +++++++++--------- CovidAlertVaccinationModel/Project.toml | 1 + .../scratch_code/model_cuda.jl | 103 ++++++++++ .../other_model_implementation.jl | 83 ++++++++ .../src/CovidAlertVaccinationModel.jl | 26 +-- CovidAlertVaccinationModel/src/agents.jl | 63 +++--- .../src/contact_vectors.jl | 40 ++++ CovidAlertVaccinationModel/src/graphs.jl | 177 ++++++----------- .../src/mixing_distributions.jl | 15 -- CovidAlertVaccinationModel/src/model.jl | 90 +++++---- CovidAlertVaccinationModel/test/runtests.jl | 25 ++- IntervalsModel/src/IntervalsModel.jl | 3 - .../src/ZeroWeightedDistributions.jl | 4 + 13 files changed, 486 insertions(+), 327 deletions(-) create mode 100644 CovidAlertVaccinationModel/scratch_code/model_cuda.jl create mode 100644 CovidAlertVaccinationModel/scratch_code/other_model_implementation.jl create mode 100644 CovidAlertVaccinationModel/src/contact_vectors.jl diff --git a/CovidAlertVaccinationModel/Manifest.toml b/CovidAlertVaccinationModel/Manifest.toml index def003c..dd67016 100644 --- a/CovidAlertVaccinationModel/Manifest.toml +++ b/CovidAlertVaccinationModel/Manifest.toml @@ -21,10 +21,10 @@ uuid = "ec485272-7323-5ecc-a04f-4719b315124d" version = "0.1.0" [[ArrayInterface]] -deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays"] -git-tree-sha1 = "ee07ae00e3cc277dcfa5507ce25be522313ecc3e" +deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] +git-tree-sha1 = "2fbfa5f372352f92191b63976d070dc7195f47a4" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "3.1.1" +version = "3.1.7" [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -40,9 +40,9 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" [[BenchmarkTools]] deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"] -git-tree-sha1 = "9e62e66db34540a0c919d72172cc2f642ac71260" +git-tree-sha1 = "068fda9b756e41e6c75da7b771e6f89fa8a43d15" uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -version = "0.5.0" +version = "0.7.0" [[Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -52,9 +52,9 @@ version = "1.0.6+5" [[CSV]] deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"] -git-tree-sha1 = "1f79803452adf73e2d3fc84785adb7aaca14db36" +git-tree-sha1 = "6d4242ef4cb1539e7ede8e01a47a32365e0a34cd" uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" -version = "0.8.3" +version = "0.8.4" [[Cairo_jll]] deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] @@ -64,39 +64,39 @@ version = "1.16.0+6" [[CategoricalArrays]] deps = ["DataAPI", "Future", "JSON", "Missings", "Printf", "Statistics", "StructTypes", "Unicode"] -git-tree-sha1 = "99809999c8ee01fa89498480b147f7394ea5450f" +git-tree-sha1 = "f713d583d10fc036252fd826feebc6c173c522a8" uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597" -version = "0.9.2" +version = "0.9.5" [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "de4f08843c332d355852721adb1592bce7924da3" +git-tree-sha1 = "44e9f638aa9ed1ad58885defc568c133010140aa" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.9.29" +version = "0.9.37" [[ColorSchemes]] deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"] -git-tree-sha1 = "3141757b5832ee7a0386db87997ee5a23ff20f4d" +git-tree-sha1 = "d3cf83862f70d430d4b34e43ed65e74bd50ae0e0" uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" -version = "3.10.2" +version = "3.11.0" [[ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "4bffea7ed1a9f0f3d1a131bbcd4b925548d75288" +git-tree-sha1 = "32a2b8af383f11cbb65803883837a149d10dfe8a" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.10.9" +version = "0.10.12" [[Colors]] -deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"] -git-tree-sha1 = "ac5f2213e56ed8a34a3dd2f681f4df1166b34929" +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "82f4e6ff9f847eca3e5ebc666ea2cd7b48e8b47e" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.6" +version = "0.12.7" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "919c7f3151e79ff196add81d7f4e45d91bbf420b" +git-tree-sha1 = "4fecfd5485d3c5de4003e19f00c6898cccd40667" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.25.0" +version = "3.26.0" [[CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -131,9 +131,9 @@ version = "1.6.0" [[DataFrames]] deps = ["CategoricalArrays", "Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] -git-tree-sha1 = "b0db5579803eabb33f1274ca7ca2f472fdfb7f2a" +git-tree-sha1 = "d50972453ef464ddcebdf489d11885468b7b83a3" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -version = "0.22.5" +version = "0.22.7" [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] @@ -165,9 +165,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[Distributions]] deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"] -git-tree-sha1 = "5a9a742ae30f13d6172c7ea245988d932134e25b" +git-tree-sha1 = "e64debe8cd174cc52d7dd617ebc5492c6f8b698c" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.24.13" +version = "0.24.15" [[Downloads]] deps = ["ArgTools", "LibCURL", "NetworkOptions"] @@ -210,9 +210,9 @@ version = "4.3.1+4" [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "e384d3cff80ac79c7a541a817192841836e46331" +git-tree-sha1 = "31939159aeb8ffad1d4d8ee44d07f8558273120a" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.11.2" +version = "0.11.7" [[FixedPointNumbers]] deps = ["Statistics"] @@ -250,27 +250,27 @@ uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" [[GLFW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"] -git-tree-sha1 = "a1bbf700b5388bffc3d882f4f4d625cf1c714fd7" +git-tree-sha1 = "bd1dbf065d7a4a0bdf7e74dd26cf932dda22b929" uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89" -version = "3.3.2+1" +version = "3.3.3+0" [[GR]] deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"] -git-tree-sha1 = "aaebdf5588281c2902f499b49e67953f2b409c9c" +git-tree-sha1 = "f74b42150042d11a1d94badf27c6125867c7e9bd" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -version = "0.54.0" +version = "0.57.1" [[GR_jll]] -deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "8aee6fa096b0cbdb05e71750c978b96a08c78951" +deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "578527027f2d6c29a8de3e2eb6887d8850ef755c" uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" -version = "0.53.0+0" +version = "0.57.1+0" [[GeometryBasics]] deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"] -git-tree-sha1 = "4d4f72691933d5b6ee1ff20e27a102c3ae99d123" +git-tree-sha1 = "4136b8a5668341e58398bb472754bff4ba0456ff" uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326" -version = "0.3.9" +version = "0.3.12" [[Gettext_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] @@ -290,10 +290,10 @@ uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe" version = "1.0.0" [[HTTP]] -deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets", "URIs"] -git-tree-sha1 = "942c1a9c750bbe79912b7bd060a420932afd35b8" +deps = ["Base64", "Dates", "IniFile", "MbedTLS", "NetworkOptions", "Sockets", "URIs"] +git-tree-sha1 = "c9f380c76d8aaa1fa7ea9cf97bddbc0d5b15adc2" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.9.3" +version = "0.9.5" [[IfElse]] git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef" @@ -378,21 +378,21 @@ uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" version = "2.10.0+3" [[LaTeXStrings]] -git-tree-sha1 = "c7aebfecb1a60d59c0fe023a68ec947a208b1e6b" +git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" -version = "1.2.0" +version = "1.2.1" [[LabelledArrays]] deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"] -git-tree-sha1 = "5e288800819c323de5897fa6d5a002bdad54baf7" +git-tree-sha1 = "df09e970c816637191ef8794ef5c5c9f8950db88" uuid = "2ee39098-c373-598a-b85f-a56591580800" -version = "1.5.0" +version = "1.6.0" [[Latexify]] deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"] -git-tree-sha1 = "3a0084cec7bf157edcb45a67fac0647f88fe5eaf" +git-tree-sha1 = "7c72983c6daf61393ee8a0b29a419c709a06cede" uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" -version = "0.14.7" +version = "0.14.12" [[LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] @@ -573,9 +573,9 @@ uuid = "91d4177d-7536-5919-b921-800302f37372" version = "1.3.1+3" [[OrderedCollections]] -git-tree-sha1 = "d45739abcfc03b51f6a42712894a593f74c80a23" +git-tree-sha1 = "4fa2ba51070ec13fcc7517db714445b4ab986bdf" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.3.3" +version = "1.4.0" [[PCRE_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -584,16 +584,16 @@ uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc" version = "8.42.0+4" [[PDMats]] -deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"] -git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc" +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "f82a0e71f222199de8e9eb9a09977bd0767d52a0" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.10.1" +version = "0.11.0" [[Parsers]] deps = ["Dates"] -git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714" +git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.15" +version = "1.1.0" [[Pipe]] git-tree-sha1 = "6842804e7867b115ca9de748a0cf6b364523c16d" @@ -624,15 +624,15 @@ version = "1.0.10" [[Plots]] deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"] -git-tree-sha1 = "cab13323a50caf17432793269677b289234f02ca" +git-tree-sha1 = "cc4eb1be2576984d7a0f7f51478827dee816138b" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.10.4" +version = "1.11.2" [[PooledArrays]] -deps = ["DataAPI"] -git-tree-sha1 = "0e8f5c428a41a81cd71f76d76f2fc3415fe5a676" +deps = ["DataAPI", "Future"] +git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414" uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" -version = "1.1.0" +version = "1.2.1" [[PrettyTables]] deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"] @@ -644,11 +644,11 @@ version = "0.11.1" deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" -[[Qt_jll]] +[[Qt5Base_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] -git-tree-sha1 = "588b799506f0b956a7e6e18f767dfa03cf333f26" -uuid = "ede63266-ebff-546c-83e0-1c6fb6d0efc8" -version = "5.15.2+2" +git-tree-sha1 = "16626cfabbf7206d60d84f2bf4725af7b37d4a77" +uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1" +version = "5.15.2+0" [[QuadGK]] deps = ["DataStructures", "LinearAlgebra"] @@ -677,9 +677,9 @@ version = "1.1.1" [[RecipesPipeline]] deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"] -git-tree-sha1 = "c4d54a78e287de7ec73bbc928ce5eb3c60f80b24" +git-tree-sha1 = "7a5026a6741c14147d1cb6daf2528a77ca28eb51" uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" -version = "0.3.1" +version = "0.3.2" [[Reexport]] git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5" @@ -687,27 +687,28 @@ uuid = "189a3867-3050-52da-a836-e630ba90ab69" version = "1.0.0" [[Referenceables]] -git-tree-sha1 = "4a32b1dd124a846580608eb347d4337f873c2499" +deps = ["Adapt"] +git-tree-sha1 = "e681d3bfa49cd46c3c161505caddf20f0e62aaa9" uuid = "42d2dcc6-99eb-4e98-b66c-637b7d73030e" -version = "0.1.0" +version = "0.1.2" [[Requires]] deps = ["UUIDs"] -git-tree-sha1 = "cfbac6c1ed70c002ec6361e7fd334f02820d6419" +git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621" uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.1.2" +version = "1.1.3" [[Rmath]] deps = ["Random", "Rmath_jll"] -git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370" +git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f" uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" -version = "0.6.1" +version = "0.7.0" [[Rmath_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "1b7bf41258f6c5c9c31df8c1ba34c1fc88674957" +git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7" uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" -version = "0.2.2+2" +version = "0.3.0+0" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -739,9 +740,9 @@ uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" [[Showoff]] deps = ["Dates", "Grisu"] -git-tree-sha1 = "ee010d8f103468309b8afac4abb9be2e18ff1182" +git-tree-sha1 = "236dd0ddad6e3764cce8d8b09c0bbba6df2e194f" uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f" -version = "0.3.2" +version = "1.0.2" [[SimpleTraits]] deps = ["InteractiveUtils", "MacroTools"] @@ -764,9 +765,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SpecialFunctions]] deps = ["ChainRulesCore", "OpenSpecFun_jll"] -git-tree-sha1 = "75394dbe2bd346beeed750fb02baa6445487b862" +git-tree-sha1 = "5919936c0e92cff40e57d0ddf0ceb667d42e5902" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "1.2.1" +version = "1.3.0" [[SplittablesBase]] deps = ["Setfield", "Test"] @@ -774,11 +775,17 @@ git-tree-sha1 = "edef25a158db82f4940720ebada14a60ef6c4232" uuid = "171d559e-b47b-412a-8079-5efa626c420e" version = "0.1.13" +[[Static]] +deps = ["IfElse"] +git-tree-sha1 = "ddec5466a1d2d7e58adf9a427ba69763661aacf6" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "0.2.4" + [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "9da72ed50e94dbff92036da395275ed114e04d49" +git-tree-sha1 = "2f01a51c23eed210ff4a1be102c4cc8236b66e5b" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.0.1" +version = "1.1.0" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -786,15 +793,15 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[StatsBase]] deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"] -git-tree-sha1 = "7bab7d4eb46b225b35179632852b595a3162cb61" +git-tree-sha1 = "a83fa3021ac4c5a918582ec4721bc0cf70b495a9" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.33.2" +version = "0.33.4" [[StatsFuns]] deps = ["Rmath", "SpecialFunctions"] -git-tree-sha1 = "3b9f665c70712af3264b61c27a7e1d62055dafd1" +git-tree-sha1 = "ced55fd4bae008a8ea12508314e725df61f0ba45" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "0.9.6" +version = "0.9.7" [[StructArrays]] deps = ["Adapt", "DataAPI", "Tables"] @@ -804,9 +811,9 @@ version = "0.5.0" [[StructTypes]] deps = ["Dates", "UUIDs"] -git-tree-sha1 = "65a43f5218197bc7091b76bc273a5e323a1d7b0d" +git-tree-sha1 = "5eaf731e88587bb72a6c1262c0a014cd1859a08d" uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" -version = "1.2.3" +version = "1.5.2" [[SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] @@ -818,15 +825,15 @@ uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" [[TableTraits]] deps = ["IteratorInterfaceExtensions"] -git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e" +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" -version = "1.0.0" +version = "1.0.1" [[Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"] -git-tree-sha1 = "a716dde43d57fa537a19058d044b495301ba6565" +git-tree-sha1 = "a9ff3dfec713c6677af435d6a6d65f9744feef67" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.3.2" +version = "1.4.1" [[Tar]] deps = ["ArgTools", "SHA"] @@ -850,9 +857,9 @@ version = "1.5.3" [[Transducers]] deps = ["ArgCheck", "BangBang", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"] -git-tree-sha1 = "9550eba57ebc2f7677c4c946aaca56e149ca73ff" +git-tree-sha1 = "721dc0b917a43e6d191f07a6ea095ced164f08ce" uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999" -version = "0.4.59" +version = "0.4.62" [[URIs]] git-tree-sha1 = "7855809b88d7b16e9b029afd17880930626f54a2" diff --git a/CovidAlertVaccinationModel/Project.toml b/CovidAlertVaccinationModel/Project.toml index cdbf10a..41bd14f 100644 --- a/CovidAlertVaccinationModel/Project.toml +++ b/CovidAlertVaccinationModel/Project.toml @@ -14,6 +14,7 @@ ImportAll = "c65182e5-40f4-518f-8165-175b85689199" Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50" NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a" Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0" diff --git a/CovidAlertVaccinationModel/scratch_code/model_cuda.jl b/CovidAlertVaccinationModel/scratch_code/model_cuda.jl new file mode 100644 index 0000000..0d41be4 --- /dev/null +++ b/CovidAlertVaccinationModel/scratch_code/model_cuda.jl @@ -0,0 +1,103 @@ +#testing some GPU implementations of the ABM model + + + +function vaccinate_uniformly!(num_vaccines,u_next,u,graph,population_list,index_vectors) #vaccination_algorithm + vaccination_indices = rand(RNG,1:length(u),num_vaccines) + u_next[vaccination_indices] .= Immunized +end +function agents_step_cuda!(t,u_next,u,population_list,graph,params,index_vectors,vaccination_algorithm!) + @unpack p, recovery_rate, vaccines_per_day,vaccination_start_day = params + u_next .= u#keep state the same if nothing else happens + + if t >= vaccination_start_day + vaccination_algorithm!(Int(floor(vaccines_per_day*length(u))),u_next,u,graph,population_list,index_vectors) + end + @inbounds for i in 1:length(u) + agent_status = u[i] + agent_demo = population_list[i] + if agent_status == Susceptible + @inbounds for j in 1:length(u) + if graph[i,j] != 0 && u[j] == Infected + if rand(RNG) < contact_weight(p,graph[i,j]) + (u_next[i] = Infected) + end + end + end + elseif agent_status == Infected + if rand(RNG) < recovery_rate + u_next[i] = Recovered + end + end + end +end + +function agents_step_gpu!(t,u_next,u,graph,params) + u_next .= u + function kernel(state, _, (t,u_next,u,g,params,randstate)) + @unpack p, recovery_rate = params + i = linear_index(state) + if i <= length(u) + agent_status = u[i] + if agent_status == Susceptible + for j in 1:length(u) + if g[i,j] == 1 && u[j] == Infected && GPUArrays.gpu_rand(Float64, state, randstate) < p + u_next[i] = Infected + end + end + elseif agent_status == Infected + if GPUArrays.gpu_rand(Float64, state, randstate) < recovery_rate + u_next[i] = Recovered + end + end + end + return nothing + end + gpu_call(kernel, u, (t,u_next,u,graph,params, GPUArrays.default_rng(typeof(u)).state)) +end + +function solve_cuda!(u_0,params,steps,agent_model,vaccination_algorithm!; record = false) + solution = vcat([u_0], [fill(Susceptible,length(u_0)) for i in 1:steps]) + + population_list = cu(agent_model.demographics) #list of demographic status for each agent + index_vectors = cu.(agent_model.demographic_index_vectors) #lists of indices of each agent with a given demographic + base_network = cu(agent_model.base_network) #network with households and LTC homes + ws_contacts = cu(agent_model.ws_matrix_list) + rest_contacts = cu(agent_model.rest_matrix_list) + + ws_daily_graph = MixingGraph(index_vectors,ws_contacts.daily,contact_time_distribution_matrix) + rest_daily_graph = MixingGraph(index_vectors,rest_contacts.daily,contact_time_distribution_matrix) + + ws_weekly_graph = MixingGraph(index_vectors,ws_contacts.twice_a_week,contact_time_distribution_matrix) + rest_weekly_graph = MixingGraph(index_vectors,rest_contacts.twice_a_week,contact_time_distribution_matrix) + + + graph = similar(base_network) + graphlist = typeof(graph)[] + for t in 1:steps + day_of_week = mod(t,7) + graph .= 0 + graph .= base_network + if !(day_of_week == 3 || day_of_week == 4) + apply_mixing_graph(graph, ws_daily_graph) + end + apply_mixing_graph(graph, rest_daily_graph) + if rand(RNG)<5/7 + apply_mixing_graph(graph, ws_weekly_graph) + apply_mixing_graph(graph, rest_weekly_graph) + end + + ws_once_graph = MixingGraph(index_vectors,ws_contacts.otherwise,contact_time_distribution_matrix) + rest_once_graph = MixingGraph(index_vectors,rest_contacts.otherwise,contact_time_distribution_matrix) + + apply_mixing_graph(graph, ws_once_graph) + apply_mixing_graph(graph, rest_once_graph) + + if record + push!(graphlist, deepcopy(graph)) + end + agents_step!(t,solution[t+1],solution[t],population_list,graph,params,index_vectors,vaccination_algorithm!) + # advance agent states based on the new network, vaccination process given by vaccination_algorithm!, which is just a function defined as above + end + return solution, graphlist +end \ No newline at end of file diff --git a/CovidAlertVaccinationModel/scratch_code/other_model_implementation.jl b/CovidAlertVaccinationModel/scratch_code/other_model_implementation.jl new file mode 100644 index 0000000..f67028b --- /dev/null +++ b/CovidAlertVaccinationModel/scratch_code/other_model_implementation.jl @@ -0,0 +1,83 @@ + +function agents_step_2!(t,S_indicator,R_indicator, I_indicies,I_indicies_next,population_list,graph,params,index_vectors,vaccination_algorithm!) + @unpack p, recovery_rate, vaccines_per_day,vaccination_start_day = params + # if t >= vaccination_start_day + # vaccination_algorithm!(Int(floor(vaccines_per_day*length(u))),u_next,u,graph,population_list,index_vectors) + # end + + while !isempty(I_indicies) + u = dequeue!(I_indicies) + for v in 1:length(population_list) + if graph[u,v] != 0 && S_indicator[v] == true + if rand(RNG) < contact_weight(p,graph[u,v]) + # display(v) + enqueue!(I_indicies_next,v) + S_indicator[v] = false + end + end + end + if rand(RNG) < recovery_rate + R_indicator[u] = true + else + enqueue!(I_indicies_next,u) + end + end +end +using BenchmarkTools +using DataStructures +function solve_2!(u_0,params,steps,agent_model,vaccination_algorithm!; record = false) + S_indicator = u_0 .== Susceptible + R_indicator = u_0 .== Recovered + I_indicies_next = Queue{Int}() + I_indicies = Queue{Int}() + + for ind in findall(==(Infected), u_0) + enqueue!(I_indicies, ind) + end + + + population_list = agent_model.demographics #list of demographic status for each agent + index_vectors = agent_model.demographic_index_vectors #lists of indices of each agent with a given demographic + base_network = agent_model.base_network #network with households and LTC homes + ws_contacts = agent_model.ws_matrix_list + rest_contacts = agent_model.rest_matrix_list + + ws_daily_graph = MixingGraph(index_vectors,ws_contacts.daily,contact_time_distribution_matrix) + rest_daily_graph = MixingGraph(index_vectors,rest_contacts.daily,contact_time_distribution_matrix) + + ws_weekly_graph = MixingGraph(index_vectors,ws_contacts.twice_a_week,contact_time_distribution_matrix) + rest_weekly_graph = MixingGraph(index_vectors,rest_contacts.twice_a_week,contact_time_distribution_matrix) + + + graph = similar(base_network) + graphlist = typeof(graph)[] + for t in 1:steps + day_of_week = mod(t,7) + graph .= 0 + graph .= base_network + if !(day_of_week == 3 || day_of_week == 4) + apply_mixing_graph(graph, ws_daily_graph) + end + apply_mixing_graph(graph, rest_daily_graph) + if rand(RNG)<5/7 + apply_mixing_graph(graph, ws_weekly_graph) + apply_mixing_graph(graph, rest_weekly_graph) + end + + ws_once_graph = MixingGraph(index_vectors,ws_contacts.otherwise,contact_time_distribution_matrix) + rest_once_graph = MixingGraph(index_vectors,rest_contacts.otherwise,contact_time_distribution_matrix) + + apply_mixing_graph(graph, ws_once_graph) + apply_mixing_graph(graph, rest_once_graph) + + if record + push!(graphlist, deepcopy(graph)) + end + agents_step_2!(t,S_indicator,R_indicator,I_indicies, I_indicies_next,population_list,graph,params,index_vectors,vaccination_algorithm!) + I_indicies, I_indicies_next = I_indicies_next, I_indicies + empty!(I_indicies_next) + + # advance agent states based on the new network, vaccination process given by vaccination_algorithm!, which is just a function defined as above + end + return R_indicator, graphlist +end \ No newline at end of file diff --git a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl index 8562317..76eb7c8 100644 --- a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl +++ b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl @@ -1,5 +1,4 @@ module CovidAlertVaccinationModel -using ZeroWeightedDistributions using LightGraphs using RandomNumbers.Xorshifts using Random @@ -8,7 +7,6 @@ using Plots using DataFrames using Distributions using StatsBase -import StatsBase:mean using Dates using ThreadsX using DelimitedFiles @@ -17,7 +15,6 @@ using NamedTupleTools using NetworkLayout:Stress using NetworkLayout:SFDP -import LightGraphs.add_edge! const PACKAGE_FOLDER = dirname(dirname(pathof(CovidAlertVaccinationModel))) @@ -25,6 +22,7 @@ const RNG = Xoroshiro128Star() include("agents.jl") include("data.jl") +include("contact_vectors.jl") include("mixing_distributions.jl") include("model.jl") include("plotting.jl") @@ -33,31 +31,23 @@ include("graphs.jl") const population = 14.57e6 #population of ontario const color_palette = palette(:seaborn_pastel) #color theme for the plots -const age_bins = [(0.0, 25.0),(25.0,65.0),(65.0,Inf)] -const infection_data = parse_cases_data() -const demographic_distribution = get_canada_demographic_distribution() - +const age_bins = [(0.0, 25.0),(25.0,65.0),(65.0,Inf)] default(dpi = 300) default(framestyle = :box) using BenchmarkTools function main() - # display(size.(agent_model.demographic_index_vectors)) - - # Random.seed!(RNG,1) Random.seed!(RNG,1) agent_model = AgentModel(5000) - steps = 300 - # @btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!); - # Random.seed!(RNG,1) - # agent_model = AgentModel(1_00) - u_0 = get_u_0(length(agent_model.demographics)) - @show length(agent_model.demographics) - @btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!); + steps = 100 + # return agent_model + sol,graphs = solve!(get_parameters(),steps,agent_model,vaccinate_uniformly!); + b1 = @benchmark solve!($get_parameters(),$steps,$agent_model,$vaccinate_uniformly!) + # display(b1) # println("done") # sol,graphs = solve_record!(u_0,get_parameters(),steps,agent_model,vaccinate_uniformly!); # return aggregate_timeseries(sol) - # plot_model_spatial_gif(agent_model.base_network,graphs,sol) + # plot_model_spatial_gif(agent_model.base_network,gddddddddraphs,sol) end end diff --git a/CovidAlertVaccinationModel/src/agents.jl b/CovidAlertVaccinationModel/src/agents.jl index 7e0a1df..e77e0ae 100644 --- a/CovidAlertVaccinationModel/src/agents.jl +++ b/CovidAlertVaccinationModel/src/agents.jl @@ -11,52 +11,57 @@ end Immunized end -struct AgentModel{GraphType,DistMatrixList1,DistMatrixList2,DistMatrixList3} +struct AgentModel{GraphType,L1,L2} demographics::Vector{AgentDemographic} - # is_app_user::Vector{Boolean} - # is_vaccinator::Vector{Boolean} - # app_user_indices::Vector{Int} - # vaccinator_indices::Vector{Int} demographic_index_vectors::Vector{Vector{Int}} base_network::GraphType - workschool_contacts_mean_adjusted_list::DistMatrixList1 - rest_contacts_mean_adjusted_list::DistMatrixList2 - home_contacts_duration_list::DistMatrixList3 + ws_matrix_list::L1 + rest_matrix_list::L2 function AgentModel(num_households) pop_list,base_network,index_vectors = generate_population(num_households) pop_sizes = length.(index_vectors) + map_symmetrize(m_tuple) = map(md -> symmetrize_means(pop_sizes,md), m_tuple) - workschool_contacts_mean_adjusted_list = [symmetrize_means(pop_sizes,initial_workschool_mixing_matrix)] - rest_contacts_mean_adjusted_list = [symmetrize_means(pop_sizes,initial_rest_mixing_matrix)] - home_contacts_duration_list = [symmetrize_means(pop_sizes,contact_time_distribution_matrix)] - - - workplace_data, distancing_data,residential_data = get_distancing_data_sets() - @assert length(workplace_data) == length(distancing_data) == length(residential_data) "uh oh your data is different lengths!" - - for t in 1:length(workplace_data) - push!(workschool_contacts_mean_adjusted_list, adjust_distributions_mean(workschool_contacts_mean_adjusted_list[1], workplace_data[t])) - push!(rest_contacts_mean_adjusted_list, adjust_distributions_mean(rest_contacts_mean_adjusted_list[1], distancing_data[t])) - push!(home_contacts_duration_list, adjust_distributions_mean(home_contacts_duration_list[1], residential_data[t])) - end - + workschool_contacts_mean_adjusted = map_symmetrize(workschool_mixing) + rest_contacts_mean_adjusted = map_symmetrize(rest_mixing) + # vaccinator_indices = sample(RNG,1:length(pop_list),length(pop_list) // 2) return new{ typeof(base_network), - typeof(workschool_contacts_mean_adjusted_list), - typeof(rest_contacts_mean_adjusted_list), - typeof(home_contacts_duration_list) + typeof(workschool_contacts_mean_adjusted), + typeof(rest_contacts_mean_adjusted) }( pop_list, index_vectors, base_network, - workschool_contacts_mean_adjusted_list, - rest_contacts_mean_adjusted_list, - home_contacts_duration_list - ) + workschool_contacts_mean_adjusted, + rest_contacts_mean_adjusted + ) end end +#generate population with households distributed according to household_size_distribution +function generate_population(num_households) + households_composition = sample_household_data(num_households) + households = map( l -> [fill(AgentDemographic(i),n) for (i,n) in enumerate(l)],households_composition) |> + l -> reduce(vcat,l) |> + l -> reduce(vcat,l) + + total_household_pop = sum(sum.(households_composition)) + household_networks = map(LightGraphs.complete_graph, sum.(households_composition)) + + population_list = reduce(vcat,households) + network = reduce(blockdiag,household_networks; init = SimpleGraph()) + index_vectors = [findall(x -> x == AgentDemographic(i), population_list) for i in 1:(AgentDemographic.size-1)] + return (; + population_list, + network, + index_vectors + ) +end + + + # func(a) = println("agent") function Base.show(io::IO, status::AgentStatus) if status == Susceptible diff --git a/CovidAlertVaccinationModel/src/contact_vectors.jl b/CovidAlertVaccinationModel/src/contact_vectors.jl new file mode 100644 index 0000000..64350cb --- /dev/null +++ b/CovidAlertVaccinationModel/src/contact_vectors.jl @@ -0,0 +1,40 @@ +function reindex!(k,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j) + i_index = index_list_i[k] + j_index = index_list_j[k] + csum += j_to_i_contacts[j_index] - i_to_j_contacts[i_index] + + i_to_j_contacts[i_index] = sample_list_i[k] + j_to_i_contacts[j_index] = sample_list_j[k] + csum += i_to_j_contacts[i_index] - j_to_i_contacts[j_index] + + return csum +end + +function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts, j_to_i_contacts) + rand!(RNG,ij_dist,i_to_j_contacts) + rand!(RNG,ji_dist,j_to_i_contacts) + l_i = length(i_to_j_contacts) + l_j = length(j_to_i_contacts) + + csum = sum(i_to_j_contacts) - sum(j_to_i_contacts) + + inner_iter = 10 + index_list_i = zeros(Int,inner_iter) + index_list_j = similar(index_list_i) + sample_list_i = zeros(eltype(i_to_j_contacts),inner_iter) + sample_list_j = similar(sample_list_i) + + while csum != 0 + index_list_i = sample!(RNG,1:l_i,index_list_i) + index_list_j = sample!(RNG,1:l_j,index_list_j) + sample_list_i = rand!(RNG,ij_dist,sample_list_i) + sample_list_j = rand!(RNG,ji_dist,sample_list_j) + for i = 1:inner_iter + if csum != 0 + csum = reindex!(i,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j) + end + end + end + + return nothing +end \ No newline at end of file diff --git a/CovidAlertVaccinationModel/src/graphs.jl b/CovidAlertVaccinationModel/src/graphs.jl index 339c446..56f490f 100644 --- a/CovidAlertVaccinationModel/src/graphs.jl +++ b/CovidAlertVaccinationModel/src/graphs.jl @@ -1,104 +1,26 @@ -#This defines a type for the time-dependent-graph -#The tradeoff is (partially) memory-usage vs speed, so it tries to preallocate as much as possible -#As result, it stores a lot of redundant information -#Julia's LightGraphs.jl package uses edgelists to store their graphs, but here I am experimenting with bitarrays, which seem to be 4-5x faster, even for huge problems. -#Also completely decouple the generation of the contact vectors (the number of contacts each node has in each WorkSchool/Rest layer) from the -#generation of the corresponding graphs -#bitmatrices are much faster at the graph generation but we lose access to neighbors in O(1) -struct MixingGraphs{GraphVector,A} - graph::GraphVector - base_graph::GraphVector - contact_vector_ws::A - contact_vector_rest::A - function MixingGraphs(static_contacts, ws_mixing_matrices, rest_mixing_matrices, index_vectors) - (length(ws_mixing_matrices) == length(rest_mixing_matrices)) || throw(ArgumentError("mixing matrix lists must be of equal length")) - ts_length = length(ws_mixing_matrices) - contact_vector_ws = map(mm -> MixingContacts(index_vectors,mm), ws_mixing_matrices) - contact_vector_rest = map(mm -> MixingContacts(index_vectors,mm), rest_mixing_matrices) - base_graph = convert(BitArray,adjacency_matrix(static_contacts)) - return new{typeof(base_graph),typeof(contact_vector_ws)}(deepcopy(base_graph),base_graph,contact_vector_ws,contact_vector_rest) - end -end - - +using LinearAlgebra +using LightGraphs #A type that defines a matrix of vectors, such that the sum of the ijth vector is equal to the sum of the jith vector struct MixingContacts{V} + total_edges::Int contact_array::V function MixingContacts(index_vectors,mixing_matrix) contacts = map(CartesianIndices(mixing_matrix)) do ind zeros(length(index_vectors[ind[1]])) end + tot = 0 for i in 1:length(mixing_matrix[:,1]), j in 1:i #diagonal generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],contacts[i,j],contacts[j,i] ) + tot += sum(contacts[i,j]) end - new{typeof(contacts)}(contacts) - end -end - -#generate a new mixing graph for the current timestep t -function advance_mixing_graph!(t,mixing_graph,index_vectors) - mixing_graph.graph .= 0 - mixing_graph.graph .= mixing_graph.base_graph - generate_mixing_graph!(mixing_graph.graph, index_vectors,mixing_graph.contact_vector_rest[t],contact_time_distribution_matrix) - generate_mixing_graph!(mixing_graph.graph, index_vectors,mixing_graph.contact_vector_ws[t],contact_time_distribution_matrix) - return nothing -end - - -function update_contacts!(mixing_contacts,mixing_matrix) - for i in 1:length(mixing_matrix[:,1]), j in 1:i #diagonal - generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],mixing_contacts.contact_array[i,j],mixing_contacts.contact_array[j,i] ) - end -end -function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts, j_to_i_contacts) - rand!(RNG,ij_dist,i_to_j_contacts) - rand!(RNG,ji_dist,j_to_i_contacts) - l_i = length(i_to_j_contacts) - l_j = length(j_to_i_contacts) - contacts_sums = sum(i_to_j_contacts) - sum(j_to_i_contacts) - sample_list_length = max(l_i,l_j) #better heuristic for this based on stddev of dist? - index_list_i = sample(RNG,1:l_i,sample_list_length) - index_list_j = sample(RNG,1:l_j,sample_list_length) - sample_list_i = rand(RNG,ij_dist,sample_list_length) - sample_list_j = rand(RNG,ji_dist,sample_list_length) - for k = 1:sample_list_length - if (contacts_sums != 0) - i_index = index_list_i[k] - j_index = index_list_j[k] - contacts_sums += j_to_i_contacts[j_index] - i_to_j_contacts[i_index] - - i_to_j_contacts[i_index] = sample_list_i[k] - j_to_i_contacts[j_index] = sample_list_j[k] - contacts_sums += i_to_j_contacts[i_index] - j_to_i_contacts[j_index] - else - break - end + new{typeof(contacts)}(tot,contacts) end - while contacts_sums != 0 - i_index = sample(RNG,1:l_i) - j_index = sample(RNG,1:l_j) - contacts_sums += j_to_i_contacts[j_index] - i_to_j_contacts[i_index] - i_to_j_contacts[i_index] = rand(RNG,ij_dist) - j_to_i_contacts[j_index] = rand(RNG,ji_dist) - contacts_sums += i_to_j_contacts[i_index] - j_to_i_contacts[j_index] - end - return nothing end -#add a bipartite graph derived from mixing matrices onto g -function generate_mixing_graph!(g,index_vectors,contacts,edge_weight_dists) - - for i in 1:length(index_vectors), j in 1:length(index_vectors) - random_bipartite_graph_fast_CL!(g,index_vectors[i],index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i]) - end - return g -end - -using StatsBase #modify g so that nodes specified in anodes and bnodes are connected by a bipartite graph with expected degrees given by aseq and bseq #implemented from Aksoy, S. G., Kolda, T. G., & Pinar, A. (2017). Measuring and modeling bipartite graphs with community structure -#simple algorithm, might produce parallel edges for small graphs -function random_bipartite_graph_fast_CL!(g::SimpleGraph,anodes,bnodes,aseq,bseq) +#simple algorithm, lightgraphs does not allow parallel edges +function random_bipartite_graph_fast_CL!(g::SimpleGraph,anodes,bnodes,aseq,bseq, weights_dict) lena = length(aseq) lenb = length(bseq) m = Int(sum(aseq)) @@ -108,40 +30,65 @@ function random_bipartite_graph_fast_CL!(g::SimpleGraph,anodes,bnodes,aseq,bseq) for k in 1:m add_edge!(g,astubs[k],bstubs[k]) end - return g + return g end -#same algorithm but with a bitarray, this cannot produce parallel edges -function random_bipartite_graph_fast_CL!(g::T,anodes,bnodes,aseq,bseq) where T<:AbstractArray - lena = length(aseq) - lenb = length(bseq) - m = Int(sum(aseq)) - @assert sum(aseq) == sum(bseq) "degree sequences must have equal sum" - astubs = sample(RNG,anodes,StatsBase.weights(aseq./m), m) - bstubs = sample(RNG,bnodes,StatsBase.weights(bseq./m), m) - for i in 1:m - g[astubs[i],bstubs[i]] = 1 - g[bstubs[i],astubs[i]] = 1 + +struct MixingGraph{G, V, M} + g::G + weights::V + covid_alert_time::M + function MixingGraph(population_list,index_vectors,mixing_matrix) + + contacts = MixingContacts(index_vectors,mixing_matrix) + + g = Graph(length(population_list)) + + weights = Dict{Tuple{Int,Int},UInt8}() + for i in 1:length(index_vectors), j in 1:i #diagonal + random_bipartite_graph_fast_CL!(g,index_vectors[i],index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i],weights) + end + + covid_alert_time = zeros(nv(g)) + return new{typeof(g),typeof(weights),typeof(covid_alert_time)}( + g, + weights, + covid_alert_time + ) + end + function MixingGraph(g::SimpleGraph) + weights = Dict{Tuple{Int,Int},UInt8}() + covid_alert_time = zeros(nv(g)) + return new{typeof(g),typeof(weights),typeof(covid_alert_time)}( + g, + weights, + covid_alert_time + ) end - return g end -#generate population with households distributed according to household_size_distribution -function generate_population(num_households) - households_composition = sample_household_data(num_households) - households = map( l -> [fill(AgentDemographic(i),n) for (i,n) in enumerate(l)],households_composition) |> - l -> reduce(vcat,l) |> - l -> reduce(vcat,l) +function MixingGraph(g::SimpleGraph,population_list,contact_time_distribution_matrix) + mg = MixingGraph(g) + sample_mixing_graph!(mg, population_list, contact_time_distribution_matrix) + return mg +end +function MixingGraph(population_list,index_vectors,mixing_matrix,contact_time_distribution_matrix) + mg = MixingGraph(population_list,index_vectors,mixing_matrix) + sample_mixing_graph!(mg, population_list, contact_time_distribution_matrix) + return mg +end - total_household_pop = sum(sum.(households_composition)) - household_networks = map(LightGraphs.complete_graph, sum.(households_composition)) +function sample_mixing_graph!(mixing_graph,population_demographics, contact_time_distributions) + for (k,e) in enumerate(edges(mixing_graph.g)) + i = src(e) + j = dst(e) + demo_i = Int(population_demographics[i]) + demo_j = Int(population_demographics[j]) - population_list = reduce(vcat,households) - network = reduce(blockdiag,household_networks; init = SimpleGraph()) - index_vectors = [findall(x -> x == AgentDemographic(i), population_list) for i in 1:(AgentDemographic.size-1)] - return (; - population_list, - network, - index_vectors - ) + contact_time = rand(RNG, contact_time_distributions[demo_i,demo_j]) + mixing_graph.covid_alert_time[i] += contact_time + mixing_graph.covid_alert_time[j] += contact_time + mixing_graph.weights[(i,j)] = contact_time + mixing_graph.weights[(j,i)] = contact_time + end end \ No newline at end of file diff --git a/CovidAlertVaccinationModel/src/mixing_distributions.jl b/CovidAlertVaccinationModel/src/mixing_distributions.jl index 60ad4a0..9002c01 100644 --- a/CovidAlertVaccinationModel/src/mixing_distributions.jl +++ b/CovidAlertVaccinationModel/src/mixing_distributions.jl @@ -1,18 +1,9 @@ - -function adjust_distributions_mean(distribution_matrix::AbstractArray{T},mean_shift_percentage) where T<:ZWDist - return map(distribution_matrix) do dist - new_mean = mean(dist)*(1 + mean_shift_percentage) - return from_mean(typeof(dist),dist.α, new_mean) - end -end - function adjust_distributions_mean(distribution_matrix::AbstractArray{T},mean_shift_percentage) where T<:Distribution return map(distribution_matrix) do dist new_mean = mean(dist)*(1 + mean_shift_percentage) return from_mean(typeof(dist), new_mean) end end - function symmetrize_means(population_sizes,mixing_matrix::Matrix{<:Distribution}) mixing_means = mean.(mixing_matrix) symmetrized_means = zeros((length(population_sizes),length(population_sizes))) @@ -30,10 +21,6 @@ end function from_mean(::Type{Poisson{T}},μ) where T return Poisson(μ) end -function from_mean(::Type{ZWDist{DistType,T}},α,μ) where {DistType <: Distribution{Univariate,W} where W, T} - return ZWDist(α,from_mean(DistType,μ/(1-α))) -end -StatsBase.mean(d::ZWDist{Dist,T}) where {Dist,T} = (1 - d.α)*StatsBase.mean(d.base_dist) function from_mean_workschool(_,μ) return Poisson(μ) end @@ -59,5 +46,3 @@ const workschool_mixing, rest_mixing = load_mixing_matrices() const contact_time_distribution_matrix = [Geometric() for i in 1:(AgentDemographic.size-1), j in 1:(AgentDemographic.size-1)] - - diff --git a/CovidAlertVaccinationModel/src/model.jl b/CovidAlertVaccinationModel/src/model.jl index d2ef5a9..3ff6785 100644 --- a/CovidAlertVaccinationModel/src/model.jl +++ b/CovidAlertVaccinationModel/src/model.jl @@ -5,8 +5,7 @@ function get_u_0(nodes) return u_0 end -function contact_weight(p, agent_v,agent_w) - contact_time = rand(RNG,contact_time_distribution_matrix[Int(agent_v),Int(agent_w)]) +function contact_weight(p, contact_time) return 1 - (1-p)^contact_time end @@ -18,20 +17,18 @@ function hotspotting!(num_vaccines,u_next,u,graph,population_list,index_vectors) # vaccination_indices = rand(RNG,1:length(u),num_vaccines) # u_next[vaccination_indices] .= Immunized end -function agents_step!(t,u_next,u,population_list,graph,params,index_vectors,vaccination_algorithm!) +function agents_step!(t,u_next,u,population_list,network_list,params,index_vectors,vaccination_algorithm!) @unpack p, recovery_rate, vaccines_per_day,vaccination_start_day = params - u_next .= u#keep state the same if nothing else happens - - if t >= vaccination_start_day - vaccination_algorithm!(Int(floor(vaccines_per_day*length(u))),u_next,u,graph,population_list,index_vectors) - end + u_next .= u @inbounds for i in 1:length(u) agent_status = u[i] agent_demo = population_list[i] if agent_status == Susceptible - @inbounds for j in 1:length(u) #LightGraphs.neighbors(graph,i) - if graph[i,j] == 1 && u[j] == Infected && rand(RNG) < contact_weight(p, agent_demo,population_list[j]) #only need contact weights for infected subnetwork - (u_next[i] = Infected) + for mixing_graph in network_list + for j in neighbors(mixing_graph.g,i) + if u[j] == Infected && rand(RNG) < contact_weight(p,mixing_graph.weights[(i,j)]) + u_next[i] = Infected + end end end elseif agent_status == Infected @@ -43,60 +40,61 @@ function agents_step!(t,u_next,u,population_list,graph,params,index_vectors,vacc end -function solve!(u_0,params,steps,agent_model,vaccination_algorithm!) - solution = vcat([u_0], [fill(Susceptible,length(u_0)) for i in 1:steps]) +function solve!(params,steps,agent_model,vaccination_algorithm!; record = false) population_list = agent_model.demographics #list of demographic status for each agent + popsize = length(population_list) index_vectors = agent_model.demographic_index_vectors #lists of indices of each agent with a given demographic - base_network = agent_model.base_network #static network, households and LTC homes - ws_contacts = agent_model.workschool_contacts_mean_adjusted_list - rest_contacts = agent_model.rest_contacts_mean_adjusted_list + ws_contacts = agent_model.ws_matrix_list + rest_contacts = agent_model.rest_matrix_list + p = params.p + u_0 = get_u_0(length(agent_model.demographics)) - @show length(ws_contacts) - mixing_graphs = MixingGraphs(base_network, ws_contacts, rest_contacts, index_vectors) - for t in 1:steps - advance_mixing_graph!(t,mixing_graphs,index_vectors) - agents_step!(t,solution[t+1],solution[t],population_list,mixing_graphs.graph,params,index_vectors,vaccination_algorithm!) - #advance agent states based on the new network, vaccination process given by vaccination_algorithm!, which is just a function defined as above - end - - return solution #mixing_graphs.graph_list -end + home_static_edges = MixingGraph(agent_model.base_network,population_list,contact_time_distribution_matrix) #network with households and LTC homes + ws_static_edges = MixingGraph(population_list,index_vectors,ws_contacts.daily,contact_time_distribution_matrix) + ws_weekly_edges = MixingGraph(population_list,index_vectors,ws_contacts.twice_a_week,contact_time_distribution_matrix) + ws_daily_edges_vector = [MixingGraph(population_list,index_vectors,ws_contacts.otherwise,contact_time_distribution_matrix) for i in 1:steps] + + rest_static_edges = MixingGraph(population_list,index_vectors,rest_contacts.daily,contact_time_distribution_matrix) + rest_weekly_edges = MixingGraph(population_list,index_vectors,rest_contacts.twice_a_week,contact_time_distribution_matrix) + rest_daily_edges_vector = [MixingGraph(population_list,index_vectors,rest_contacts.otherwise,contact_time_distribution_matrix) for i in 1:steps] + -function solve_record!(u_0,params,steps,agent_model,vaccination_algorithm!) solution = vcat([u_0], [fill(Susceptible,length(u_0)) for i in 1:steps]) - population_list = agent_model.demographics #list of demographic status for each agent - index_vectors = agent_model.demographic_index_vectors #lists of indices of each agent with a given demographic - base_network = agent_model.base_network #static network, households and LTC homes - ws_contacts = agent_model.workschool_contacts_mean_adjusted_list - rest_contacts = agent_model.rest_contacts_mean_adjusted_list - mixing_graphs = MixingGraphs(base_network, ws_contacts, rest_contacts, index_vectors) - graphs = SimpleGraph[] + network_lists = [ + [home_static_edges,rest_static_edges] for i in 1:steps + ] + + for (t,l) in enumerate(network_lists) + day_of_week = mod(t,7) + if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess + push!(l, ws_static_edges) + end + if rand(RNG)<5/7 + push!(l, ws_weekly_edges) + push!(l, rest_weekly_edges) + end + push!(l,ws_daily_edges_vector[t]) + push!(l,rest_daily_edges_vector[t]) + end + for t in 1:steps - advance_mixing_graph!(t,mixing_graphs,index_vectors) - g_matrix_graph = SimpleGraph(deepcopy(mixing_graphs.graph)) - push!(graphs,g_matrix_graph) - agents_step!(t,solution[t+1],solution[t],population_list,mixing_graphs.graph,params,index_vectors,vaccination_algorithm!) + network_list = network_lists[t] #advance agent states based on the new network, vaccination process given by vaccination_algorithm!, which is just a function defined as above + agents_step!(t,solution[t+1],solution[t],population_list,network_list,params,index_vectors,vaccination_algorithm!) end - - return solution,graphs + return solution, network_lists end function get_parameters() params = ( p = 0.05, - recovery_rate = 0.05, + recovery_rate = 0.1, vaccines_per_day = 0.0, #percentage of population vaccinated per day vaccination_start_day = 60 ) return params end -function fit_parameters() - - -end - diff --git a/CovidAlertVaccinationModel/test/runtests.jl b/CovidAlertVaccinationModel/test/runtests.jl index ae7b6c5..b8ec529 100644 --- a/CovidAlertVaccinationModel/test/runtests.jl +++ b/CovidAlertVaccinationModel/test/runtests.jl @@ -15,8 +15,8 @@ dem_cat = AgentDemographic.size -1 #covidalert @testset "mixing matrices, size: $sz" for (m,sz) in zip(agent_models,model_sizes) - ws_dist = m.workschool_contacts_mean_adjusted_list - r_dist = m.rest_contacts_mean_adjusted_list + ws_dist = m.ws_matrix_list + r_dist = m.rest_matrix_list index_vec =m.demographic_index_vectors @testset "workschool" for i in dem_cat, j in dem_cat for t in 1:length(ws_dist) @@ -35,7 +35,7 @@ function vac_rate_test(model,vac_strategy, vac_rate; rng = Xoroshiro128Plus()) u_0 = get_u_0(length(model.demographics)) params = merge(get_parameters(),(vaccines_per_day = vac_rate,)) steps = 300 - sol1 = solve!(u_0,params,steps,model,vac_strategy); + sol1,_ = solve!(u_0,params,steps,model,vac_strategy); total_infections = count(x->x == AgentStatus(3),sol1[end]) display(total_infections) return total_infections @@ -43,14 +43,13 @@ end function infection_rate_test(model, inf_parameter; rng = Xoroshiro128Plus()) - u_0 = get_u_0(length(model.demographics)) params = merge(get_parameters(),(p = inf_parameter,)) steps = 300 - display(params) - sol1 = solve!(u_0,params,steps,model,vaccinate_uniformly!); + # display(params) + sol1,_ = solve!(params,steps,model,vaccinate_uniformly!); total_infections = count(x->x == AgentStatus(3),sol1[end]) - display(total_infections) + # display(total_infections) return total_infections end @@ -60,12 +59,12 @@ function test_comparison(f,xpts,comparison) return all(comparison(ypts[i],ypts[i+1]) for i in 1:length(ypts)-1) end -@testset "vaccination efficacy $sz" for (m,sz) in zip(deepcopy(agent_models),model_sizes) - @show vac_rate_test(m,vaccination_strategies[1],vaccination_rates[1]) - @testset "strategy" for strat in vaccination_strategies - @test test_comparison(x->vac_rate_test(m,strat,x),vaccination_rates,>) - end -end +# @testset "vaccination efficacy $sz" for (m,sz) in zip(deepcopy(agent_models),model_sizes) +# @show vac_rate_test(m,vaccination_strategies[1],vaccination_rates[1]) +# @testset "strategy" for strat in vaccination_strategies +# @test test_comparison(x->vac_rate_test(m,strat,x),vaccination_rates,>) +# end +# end @testset "infection efficacy $sz" for (m,sz) in zip(deepcopy(agent_models),model_sizes) @test test_comparison(x->infection_rate_test(m,x),infection_rates,<) diff --git a/IntervalsModel/src/IntervalsModel.jl b/IntervalsModel/src/IntervalsModel.jl index 936a9d8..1edbc94 100644 --- a/IntervalsModel/src/IntervalsModel.jl +++ b/IntervalsModel/src/IntervalsModel.jl @@ -1,7 +1,5 @@ module IntervalsModel export main, do_hh, do_ws, do_rest - - using Intervals using CSV using DataFrames @@ -10,7 +8,6 @@ using Random using KernelDensity using StatsBase using Distributions -using ZeroWeightedDistributions const PACKAGE_FOLDER = dirname(dirname(pathof(IntervalsModel))) using DataStructures:OrderedDict using Serialization diff --git a/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl b/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl index 4d2c1ab..3eda876 100644 --- a/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl +++ b/ZeroWeightedDistributions/src/ZeroWeightedDistributions.jl @@ -22,6 +22,10 @@ struct ZWDist{BaseDistType <: Sampleable,T} <: Sampleable{Univariate,T} end end +# function from_mean(::Type{ZWDist{DistType,T}},α,μ) where {DistType <: Distribution{Univariate,W} where W, T} +# return ZWDist(α,from_mean(DistType,μ/(1-α))) +# end +StatsBase.mean(d::ZWDist{Dist,T}) where {Dist,T} = (1 - d.α)*StatsBase.mean(d.base_dist) function Distributions.pdf(d::ZWDist, x) -- GitLab