diff --git a/.gitignore b/.gitignore index dd0d139..91d3439 100644 --- a/.gitignore +++ b/.gitignore @@ -54,4 +54,7 @@ Icon .AppleDesktop Network Trash Folder Temporary Items -.apdisk \ No newline at end of file +.apdisk + +.vscode +.vscode/* \ No newline at end of file diff --git a/docs/Manifest.toml b/docs/Manifest.toml deleted file mode 100644 index aeee0c6..0000000 --- a/docs/Manifest.toml +++ /dev/null @@ -1,512 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.10.0" -manifest_format = "2.0" -project_hash = "ccf1569f817ef6784727f5a958ab1673869da75b" - -[[deps.ANSIColoredPrinters]] -git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" -uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9" -version = "0.0.1" - -[[deps.AbstractTrees]] -git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177" -uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" -version = "0.4.5" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.1" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[deps.Calculus]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad" -uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" -version = "0.5.1" - -[[deps.CodecZlib]] -deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73" -uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.4" - -[[deps.Combinatorics]] -git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860" -uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" -version = "1.0.2" - -[[deps.Compat]] -deps = ["TOML", "UUIDs"] -git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.14.0" -weakdeps = ["Dates", "LinearAlgebra"] - - [deps.Compat.extensions] - CompatLinearAlgebraExt = "LinearAlgebra" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+1" - -[[deps.DataAPI]] -git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" -uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.16.0" - -[[deps.DataStructures]] -deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "1fb174f0d48fe7d142e1109a10636bc1d14f5ac2" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.17" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[deps.Distributions]] -deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "7c302d7a5fec5214eb8a5a4c466dcf7a51fcf169" -uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.107" - - [deps.Distributions.extensions] - DistributionsChainRulesCoreExt = "ChainRulesCore" - DistributionsDensityInterfaceExt = "DensityInterface" - DistributionsTestExt = "Test" - - [deps.Distributions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" - Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" - -[[deps.Documenter]] -deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"] -git-tree-sha1 = "4a40af50e8b24333b9ec6892546d9ca5724228eb" -uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.3.0" - -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.DualNumbers]] -deps = ["Calculus", "NaNMath", "SpecialFunctions"] -git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566" -uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" -version = "0.6.8" - -[[deps.Expat_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c" -uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.5.0+0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" - -[[deps.FillArrays]] -deps = ["LinearAlgebra", "Random"] -git-tree-sha1 = "5b93957f6dcd33fc343044af3d48c215be2562f1" -uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "1.9.3" -weakdeps = ["PDMats", "SparseArrays", "Statistics"] - - [deps.FillArrays.extensions] - FillArraysPDMatsExt = "PDMats" - FillArraysSparseArraysExt = "SparseArrays" - FillArraysStatisticsExt = "Statistics" - -[[deps.Fleck]] -deps = ["Combinatorics", "DataStructures", "Distributions", "Documenter", "Logging", "Random", "SafeTestsets"] -path = ".." -uuid = "5bb9b785-358c-4fee-af0f-b94a146244a8" -version = "0.1.0" - -[[deps.Git]] -deps = ["Git_jll"] -git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e" -uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" -version = "1.3.1" - -[[deps.Git_jll]] -deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "12945451c5d0e2d0dca0724c3a8d6448b46bbdf9" -uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" -version = "2.44.0+1" - -[[deps.HypergeometricFunctions]] -deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] -git-tree-sha1 = "f218fe3736ddf977e0e772bc9a586b2383da2685" -uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.23" - -[[deps.IOCapture]] -deps = ["Logging", "Random"] -git-tree-sha1 = "8b72179abc660bfab5e28472e019392b97d0985c" -uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.4" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[deps.IrrationalConstants]] -git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.2" - -[[deps.JLLWrappers]] -deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.5.0" - -[[deps.JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" -uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.4" - -[[deps.LazilyInitializedFields]] -git-tree-sha1 = "8f7f3cabab0fd1800699663533b6d5cb3fc0e612" -uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf" -version = "1.2.2" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.4" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.4.0+0" - -[[deps.LibGit2]] -deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] -uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.6.4+0" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[deps.Libiconv_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" -uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.17.0+0" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[deps.LogExpFunctions]] -deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "18144f3e9cbe9b15b070288eef858f71b291ce37" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.27" - - [deps.LogExpFunctions.extensions] - LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" - LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" - LogExpFunctionsInverseFunctionsExt = "InverseFunctions" - - [deps.LogExpFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[deps.MarkdownAST]] -deps = ["AbstractTrees", "Markdown"] -git-tree-sha1 = "465a70f0fc7d443a00dcdc3267a497397b8a3899" -uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391" -version = "0.1.2" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+1" - -[[deps.Missings]] -deps = ["DataAPI"] -git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" -uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.1.0" - -[[deps.Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2023.1.10" - -[[deps.NaNMath]] -deps = ["OpenLibm_jll"] -git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" -uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.0.2" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+2" - -[[deps.OpenLibm_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+2" - -[[deps.OpenSSL_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "60e3045590bd104a16fefb12836c00c0ef8c7f8c" -uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.13+0" - -[[deps.OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" -uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" - -[[deps.OrderedCollections]] -git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.3" - -[[deps.PCRE2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" -version = "10.42.0+1" - -[[deps.PDMats]] -deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65" -uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.31" - -[[deps.Parsers]] -deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" -uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.8.1" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.10.0" - -[[deps.PrecompileTools]] -deps = ["Preferences"] -git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f" -uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.2.0" - -[[deps.Preferences]] -deps = ["TOML"] -git-tree-sha1 = "9e8fed0505b0c15b4c1295fd59ea47b411c019cf" -uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.4.2" - -[[deps.Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" - -[[deps.QuadGK]] -deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "9b23c31e76e333e6fb4c1595ae6afa74966a729e" -uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.9.4" - -[[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.Random]] -deps = ["SHA"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[deps.Reexport]] -git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" -uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "1.2.2" - -[[deps.RegistryInstances]] -deps = ["LazilyInitializedFields", "Pkg", "TOML", "Tar"] -git-tree-sha1 = "ffd19052caf598b8653b99404058fce14828be51" -uuid = "2792f1a3-b283-48e8-9a74-f99dce5104f3" -version = "0.1.0" - -[[deps.Rmath]] -deps = ["Random", "Rmath_jll"] -git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b" -uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" -version = "0.7.1" - -[[deps.Rmath_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da" -uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" -version = "0.4.0+0" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -version = "0.7.0" - -[[deps.SafeTestsets]] -git-tree-sha1 = "81ec49d645af090901120a1542e67ecbbe044db3" -uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -version = "0.1.0" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[deps.Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[deps.SortingAlgorithms]] -deps = ["DataStructures"] -git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" -uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.2.1" - -[[deps.SparseArrays]] -deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -version = "1.10.0" - -[[deps.SpecialFunctions]] -deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d" -uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.3.1" - - [deps.SpecialFunctions.extensions] - SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" - - [deps.SpecialFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - -[[deps.Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -version = "1.10.0" - -[[deps.StatsAPI]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" -uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.7.0" - -[[deps.StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "1d77abd07f617c4868c33d4f5b9e1dbb2643c9cf" -uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.34.2" - -[[deps.StatsFuns]] -deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] -git-tree-sha1 = "cef0472124fab0695b58ca35a77c6fb942fdab8a" -uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "1.3.1" - - [deps.StatsFuns.extensions] - StatsFunsChainRulesCoreExt = "ChainRulesCore" - StatsFunsInverseFunctionsExt = "InverseFunctions" - - [deps.StatsFuns.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.SuiteSparse]] -deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] -uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" - -[[deps.SuiteSparse_jll]] -deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] -uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "7.2.1+1" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - -[[deps.Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[deps.TranscodingStreams]] -git-tree-sha1 = "54194d92959d8ebaa8e26227dbe3cdefcdcd594f" -uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.10.3" -weakdeps = ["Random", "Test"] - - [deps.TranscodingStreams.extensions] - TestExt = ["Test", "Random"] - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[deps.Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+1" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+1" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.52.0+1" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+2" diff --git a/docs/Project.toml b/docs/Project.toml index fa059a8..a5d7782 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,8 @@ [deps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Fleck = "5bb9b785-358c-4fee-af0f-b94a146244a8" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/make.jl b/docs/make.jl index 213db55..d18dfca 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,15 @@ using Fleck using Documenter +using Literate + +example_base = joinpath(dirname(@__FILE__), "src") +println("example base is $example_base") +Literate.markdown( + joinpath(example_base, "simple_board.jl"), + example_base, + name="mainloop", + execute=true + ) makedocs(; modules=[Fleck], @@ -13,10 +23,17 @@ makedocs(; ), pages=[ "Home" => "index.md", - "Structure" => "objects.md", - "Develop" => "develop.md", - "Delay Equations" => "delay.md", - "Vector Addition Systems" => "vas.md", + "Guide" => [ + "guide.md", + "mainloop.md", + "distributions.md" + ], + "Manual" => [ + "Structure" => "objects.md", + "Delay Equations" => "delay.md", + "Develop" => "develop.md", + "Vector Addition Systems" => "vas.md", + ], "Reference" => "reference.md" ], ) diff --git a/docs/notes/DistPlots.ipynb b/docs/notes/DistPlots.ipynb new file mode 100644 index 0000000..3e29efb --- /dev/null +++ b/docs/notes/DistPlots.ipynb @@ -0,0 +1,523 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "471a5dda", + "metadata": {}, + "outputs": [], + "source": [ + "using Plots\n", + "using Distributions\n", + "using LaTeXStrings\n", + "using QuadGK\n", + "using ForwardDiff" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "aefdaf97", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Plot{Plots.GRBackend() n=1}" + ] + } + ], + "source": [ + "x = 0:0.01:15\n", + "y = pdf.(Gamma(5.0, 1.0), x)\n", + "p = plot(x, y, label=\"Gamma(5,1)\")\n", + "xlabel!(p, \"Time\")\n", + "ylabel!(p, \"Probability Density Function\")\n", + "title!(p, \"PDF of a Gamma Distribution\")\n", + "savefig(p, \"gammapdf.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7cbc9372", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "shiftpdf (generic function with 1 method)" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "survival = ccdf\n", + "hazard(dist, x) = pdf(dist, x) / survival(dist, x)\n", + "conditional_survival(dist, intermediate, x) = survival(dist, x) / survival(dist, intermediate)\n", + "shiftpdf(dist, intermediate, x) = hazard(dist, x) * conditional_survival(dist, intermediate, x)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "c70e5181", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\"/home/adolgert/dev/Fleck/docs/notes/gammahazard.png\"" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = 0:0.01:15\n", + "y = hazard.(Gamma(5.0, 1.0), x)\n", + "p = plot(x, y, label=\"Gamma(5,1)\")\n", + "xlabel!(p, \"Time\")\n", + "ylabel!(p, \"Hazard Rate [per unit time]\")\n", + "title!(p, \"Hazard of a Gamma Distribution\")\n", + "savefig(p, \"gammahazard.png\")" + ] + }, + { + "cell_type": "markdown", + "id": "91d3514a", + "metadata": {}, + "source": [ + "OK, now shift this distribution and see how it changes." + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "c8e47a9b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\"/home/adolgert/dev/Fleck/docs/notes/shiftedgamma.png\"" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t=5.0\n", + "\n", + "x1 = 5:0.01:15\n", + "y1 = shiftpdf.(Gamma(5.0, 1.0), t, x1)\n", + "plot(x1, y1, label=\"Shifted Gamma(5,1)\")\n", + "x = 0:0.01:15\n", + "y = pdf.(Gamma(5.0, 1.0), x)\n", + "plot!(x, y, label=\"Gamma(5,1)\")\n", + "vline!([t], label=\"New start time\")\n", + "xlabel!(\"Time\")\n", + "ylabel!(\"Probability Density Function\")\n", + "title!(\"Change to PDF After Not Firing\")\n", + "savefig(\"shiftedgamma.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "3340aa16", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\"/home/adolgert/dev/Fleck/docs/notes/individual_distributions.png\"" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "com_dists = [\n", + " Exponential(10),\n", + " Weibull(1.8, 8),\n", + " Gamma(5.0, 1.5)\n", + "]\n", + "com_diff_dists = [\n", + " Exponential(10.0),\n", + " Weibull(1.8, 8),\n", + " Gamma{ForwardDiff.Dual}(5.0, 1.5)\n", + "]\n", + "com_labels = [\"Exponential(10)\" \"Weibull(1.8, 8)\" \"Gamma(5,1.5)\"]\n", + "com_x = 0:0.01:15\n", + "com_hazard = stack([hazard.(d, com_x) for d in com_dists])\n", + "a = plot(com_x, com_hazard,\n", + " label=com_labels,\n", + " title=\"Hazards of Individual Distributions\")\n", + "com_pdf = stack([pdf.(d, com_x) for d in com_dists])\n", + "b = plot(com_x, com_pdf,\n", + " label=[\"Exponential(10)\" \"Weibull(1.8, 8)\" \"Gamma(5,1.5)\"],\n", + " title=\"PDFs of Individual Distributions\"\n", + " )\n", + "plot(a, b, layout=(1, 2), size=(1000, 400))\n", + "savefig(\"individual_distributions.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "80280377", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "tp = [0.38787017893511294, 0.35009751360605557, 0.26211976998246983]\n" + ] + }, + { + "data": { + "text/plain": [ + "\"/home/adolgert/dev/Fleck/docs/notes/conditional_on_which.png\"" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function marginal_probability_event(dists, which, t)\n", + " core_matrix = pdf(dists[which], t)\n", + " for i in eachindex(dists)\n", + " if i != which\n", + " core_matrix *= survival(dists[i], t)\n", + " end\n", + " end\n", + " core_matrix\n", + "end\n", + "\n", + "total_probability(which, lim) = quadgk(t -> marginal_probability_event(com_dists, which, t), 0, lim, rtol=1e-3)[1]\n", + "tp = total_probability.([1, 2, 3], Inf)\n", + "@show tp\n", + "bb = bar(reshape(com_labels, 3), tp, legend=:none, ylabel=\"Probability of Each Clock\")\n", + "com_pdfs = zeros(Float64, length(com_x), length(com_dists))\n", + "for pdfidx in eachindex(com_dists)\n", + " for (ix, x) in enumerate(com_x)\n", + " com_pdfs[ix, pdfidx] = marginal_probability_event(com_dists, pdfidx, x)\n", + " end\n", + "end\n", + "bp = plot(com_x, com_pdfs, labels=com_labels, title=\"PDFs of Competing Clocks\")\n", + "plot(bb, bp, size=(1000, 400))\n", + "savefig(\"conditional_on_which.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "d78e9f35", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\"/home/adolgert/dev/Fleck/docs/notes/conditional_on_when.png\"" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function total_pdf(t)\n", + " factor1 = zero(Float64)\n", + " for dist in com_dists\n", + " factor1 += hazard(dist, t)\n", + " end\n", + " factor2 = one(Float64)\n", + " for dist in com_dists\n", + " factor2 *= survival(dist, t)\n", + " end\n", + " return factor1 * factor2\n", + "end\n", + "a = plot(com_x, total_pdf.(com_x), legend=:none, title=\"PDF for First to Fire\")\n", + "\n", + "function com_cond_prob(i, t)\n", + " hazard(com_dists[i], t) / sum(hazard.(com_dists, t))\n", + "end\n", + "com_rel_prob = stack([com_cond_prob.(i, com_x) for i in eachindex(com_dists)])\n", + "b = plot(com_x, com_rel_prob, labels=com_labels, ylabel=\"Conditional Probability\",\n", + " title=\"Which Event Fires Given Firing Time\")\n", + "plot(a, b, size=(1000, 400))\n", + "savefig(\"conditional_on_when.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "602ad494", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\"/home/adolgert/dev/Fleck/docs/notes/competinghazards.png\"" + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = 0:0.01:15\n", + "yw = hazard.(Weibull(1.8, 8), x)\n", + "plot(x, yw, label=\"Weibull(2, 1)\")\n", + "yg = hazard.(Gamma(5.0, 1.5), x)\n", + "plot!(x, yg, label=\"Gamma(5,1.5)\")\n", + "ye = hazard.(Exponential(4), x)\n", + "plot!(x, ye, label=\"Exponential(4)\")\n", + "xlabel!(\"Time\")\n", + "ylabel!(\"Hazard [per unit time]\")\n", + "title!(\"Three Hazard Rates\")\n", + "savefig(\"competinghazards.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": 109, + "id": "59c99c70", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", + "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/dev/Fleck/docs/notes/Project.toml`\n", + " \u001b[90m[f6369f11] \u001b[39m\u001b[92m+ ForwardDiff v0.10.36\u001b[39m\n", + "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/dev/Fleck/docs/notes/Manifest.toml`\n", + " \u001b[90m[bbf7d656] \u001b[39m\u001b[92m+ CommonSubexpressions v0.3.0\u001b[39m\n", + " \u001b[90m[163ba53b] \u001b[39m\u001b[92m+ DiffResults v1.1.0\u001b[39m\n", + " \u001b[90m[b552c78f] \u001b[39m\u001b[92m+ DiffRules v1.15.1\u001b[39m\n", + " \u001b[90m[f6369f11] \u001b[39m\u001b[92m+ ForwardDiff v0.10.36\u001b[39m\n", + " \u001b[90m[1e83bf80] \u001b[39m\u001b[92m+ StaticArraysCore v1.4.2\u001b[39m\n", + "\u001b[32m\u001b[1mPrecompiling\u001b[22m\u001b[39m project...\n", + "\u001b[32m ✓ \u001b[39mForwardDiff\n", + " 1 dependency successfully precompiled in 2 seconds. 178 already precompiled.\n" + ] + } + ], + "source": [ + "using Pkg\n", + "Pkg.add(\"ForwardDiff\")" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "842cbbe0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1×3 Matrix{Int64}:\n", + " 1 2 3" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a=[1 2 3]" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "6d04fe6b", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", + "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/dev/Fleck/docs/notes/Project.toml`\n", + " \u001b[90m[1fd47b50] \u001b[39m\u001b[92m+ QuadGK v2.9.4\u001b[39m\n", + "\u001b[32m\u001b[1m No Changes\u001b[22m\u001b[39m to `~/dev/Fleck/docs/notes/Manifest.toml`\n" + ] + } + ], + "source": [ + "using Pkg\n", + "Pkg.add(\"QuadGK\")" + ] + }, + { + "cell_type": "code", + "execution_count": 83, + "id": "13d05d0e", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAASR0lEQVR4nO3da4yU9b3A8ZnZmYW9VAE1sGaPq7soGm09KpcaFVcBY9/QqBS0Gosm1pimDbyxSqO22jQYozno8XJc24o0RcVqI1SspTGeXg7WC8ISXS5VRBcSoii4UtxlZs6LPeUYKzqrO8+z7O/zefXM+s/z/2Uzu1+feWaHbLlczgBAVLm0BwCANAkhAKEJIQChCSEAoQkhAKEJIQChDTiEixYtWrRoUTVGAYDkZQf0d4SPP/741Vdf3dzcvGbNmurNBACJGcAV4bvvvnvDDTf86Ec/qt40AJCwfOVL582bd+211xYKhepNAwAJqzSETz31VHd39+WXX7506dLPXnnnnXc+8sgjTU1N/Q+z2ezChQuPOuqoLzXmP234+5aXtu4clFMNHaVSKZPJ5HLD7Y1Lxx3RMPGkCWlPcbAql8ulUqmmpibtQRhCyuVysVjM5wdwAUM+n89ms5+zppIT7dq1a/78+StXrvzc02UymZ07dx522GGzZ8/uf1goFJqamgbr5/n6G25a8d5hmfFnDMrZqKIdm497Y/lrL/1P2nMcrIrFYrFYrK2tTXsQhpBSqdTb2ztixIi0BzmYVJKtikL48ssvv/XWWzNmzMhkMj09Pe+//35bW9v69evr6ur+dXFDQ8Pxxx+/P4SDK5vNZo6ZnJl4UTVOzmB644XMG8uH32VuYsrlci6X8w3kEzwrqqGiEJ511lnd3d39x4899thdd9313HPPjRw5spqDAUASKgphPp8fPXp0/3FDQ0NNTc3+hwBwUBvwTdeZM2dOnTq1GqMAQPIGHMLGxsbGxsZqjAIAyXPTFYDQhBCA0IQQgNCEEIDQhBCA0HxmHcBws2fPnquuuqqvry/tQQZNQ0PDL37xi0o+L+0LEEKA4ea99957+umn77vvvrQHGTSXXHJJR0dHlT5wXAgBhqG6urpvfetbaU8xaL797W9X7+TuEQIQmhACEJoQAhCaEAIQmhACEJp3jQKEMOWsc9avW5vYdrmamvvu+o9LL700sR2/MCEECOHvr7+xZ/6qzJjmZLbLr7ilu7u7wsVr1qxZvXr1m2++OWvWrIkTJ1Z1sH8lhABh1B2aqR+d0F6FkZWvXbBgwahRo/70pz8dd9xxyYfQPUIAEvK73/2uq6tr/8OVK1e++uqr/QdLly5taWlJZSohBCAhGzduvOGGG/qPe3p6LrnkkkKhkO5IGSEEIDFXXHHFH/7wh23btmUymaVLl06ePPnYY49NeyghBCApo0aNuuCCCxYvXpzJZDo6Oq6++uq0J8pkhBCAJF1zzTUdHR1r1qzZunXrzJkz0x4nkxFCAJI0efLkMWPGzJ0798orrxwKNwgzQghAwq655prOzs4rr7zy419pa2t7+eWXr7vuura2tmeffTbJefwdIUAI2Wym8Mxt2bpDktmuvOG/s+2zPvU/7d2797zzzhs/fvz+r9x2220/+9nP9j9sbGys+nwfI4QAIfzn7bdu2bIlse1yZ8+64IILPvHFd95557e//e1Pf/rThx9++ONfT7h8nyCEACHMmTMn7REyH3744ebNm++///6zzz477Vn+nxACkJCWlpaFCxemPcUnebMMAKEJIQChCSEAoQkhAKF5swzAcFNTU7Njx47k/2G/6imXy9lstkonF0KA4WbcuHEvvvhiX19f2oMMmvr6+pqamiqdXAgBhqGvfe1raY9w0HCPEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNDyFa7btGnTAw88sGnTpvr6+nPPPfc73/lOTU1NVScDgARUGsLt27cffvjhZ5555vvvv3/zzTe/8cYbt9xyS1UnA4AEVBrCqVOnTp06tf+4WCzec889QgjAMDDge4Td3d1PPPHEueeeW41pACBhlV4RZjKZzZs3T548+b333ps8efLixYsPtGzTpk2///3vX3zxxf6H2Wx20aJFRx999JcctN++ffsG5TwkoFQq9fT0pD3FwapYLPb29haLxbQHOVg9uuyxOx94KO0pBl+5XMpmh9ubHE+aMP6+O++o0snr6+tzuc/5jg0ghOPHj9+5c2dPT8/8+fNnzZq1atWqT13W3Nz89a9//aqrrup/mM1mJ0yYUCgUKt/oM+TzAxiYdOVyucbGxrSnOFj1h7Curi7tQQ5Wa9Z1ri0cmzntgrQH4fPs3vH2Uzem+7tiwF1pbGycN2/eySefXCwWP/WNo3V1dS0tLdOnTx+M8QC+qLHjMydMS3sIPs+7W9OeoOJ7hG+99Va5XO4/fvLJJydMmODPJwAYBiq9Ily4cOGTTz559NFHb9++PZPJ/PrXv67mVACQkEpDePfdd19//fVvv/32mDFj2traXA4CMDwM4B5hc3Nzc3Nz9UYBgOQNt7fhAsCACCEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQCh5Stct3379l/+8pfPP/98qVRqb2//3ve+N3LkyKpOBgAJqDSEzz333NatW+fOnVtbW3vjjTd2dnY++OCD1RwMAJJQaQgvvvjiiy++uP+4sbHxm9/8ZtVGAoDkfJF7hBs3bmxpaRn0UQAgeZVeEe735ptvLliw4KGHHjrQgq6urhUrVqxatar/YaFQ6OjoaG1t/eIzfsy+ffsG5TwkoFQq9fT0JLDRvf/Vcf+DSxLYKEnlcqacKeey2bQHGWSnT5l0z3/cnsBGfX19CezCoChnytX7XVFfX5/Lfc4l38BCuH379hkzZixYsOAb3/jGgda0traef/75P/jBD/ofZrPZE088saamZkAbHUg+P+Byk5ZcLtfY2JjARute7drccn7mFC/XD3nvbNm76uZknhWFQiGBXRgU2Uw2mWfFgQygKzt27Jg+ffrcuXPnz5//Gctqa2vHjh172mmnfenZoGJj/i3TcmraQ/B5CnVpTwCfotJ7hO+888706dNnz569YMGCqg4EAEmqNIS/+tWvOjs7Fy1aNOaf9uzZU9XJACABlYZw3rx55XJ558fU19dXdTIASICPWAMgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEIDQhBCA0IQQgNCEEILR8hes++uijRx999KWXXtq2bdtdd901duzYqo4FAMmo9Ipw9+7djz76aGNj47Jlyz788MOqzgQAiak0hEccccTy5ctvvPHGqk4DAAlzjxCA0Cq9R1i511577fHHH1+2bFn/w0Kh8Mgjj4wfP35QTr5v375BOQ8JKJVKH3zwQQIb9fX1JbALg6JUTuhZ0dvrWXHQKJfL1XtW1NfX19TUfPaawQ/hscceO3v27Ouuu27/V4455phsNjsoJ8/nB39gqiSXy33lK19JYKNCoZDxS+8gkcsm9KyorS0ksAuDIpvNJvOsOJDB70o+nx81alRra+ugnxkABp17hACENoAQjhkzZsSIEZlMpq2tLZvN+iMKAIaBAbw0unPnzurNAQCp8NIoAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChDSCES5YsmTJlyimnnLJo0aLqDQQAScpXuO7Pf/7zvHnzfvOb3xxyyCEXXnjhuHHj5syZU9XJACABlV4R3nvvvd/97nfb29tPPfXUa6+99t57763qWACQjEpDuH79+okTJ/YfT5o0qbOzs2ojAUBysuVyuZJ1Rx555JIlS6ZNm5bJZDZt2jRhwoSPPvqoUCj868rzzjvvj3/8Yy73f4nNZrNjx4791JVfwK7eck/TydnCyEE52xBRLBbL5XI+X+nL1AeHYt+IbZ2HFfYlsNUHHxV3N/17trYugb0SM0yfFftqt3UeXuhLYKuevb27mk7NjqhPYK/ElIqlUqmULwyvZ0WpWOjuPKLQW6XTX3rppbfccstnr6n0G3rooYf29PT0H3/wwQcNDQ0HatsTTzyxevXq/T/A+Xy+qampwl1iGp6/8vhySqVSsVgcrP+DZHgol8t9fX21tbVpD3IwqSRAlf7ybW1t3bhxY//xxo0bW1tbD7SyoaGh/8IRAIa+Su8RXnbZZT//+c93797d29t79913X3bZZVUdCwCSUWkI58yZ097eftRRRzU1NY0bN+773/9+VccCgGRU+maZfv/4xz+KxWJjY2P1BgKAJA0shAAwzPisURharrjiir/85S9pTwGBCCEMLV1dXbt27Up7CgjE366lbN26da+88ko+n586dWpzc3Pa4zBUrFq1qru7u729vaWlJe1ZGCo2btz4t7/9LZfLtbe3H3nkkWmPM3y4IkzT2rVr58+f39XVtXr16kmTJq1evTrtiRgSbrrppiVLlrzyyiuTJk3661//mvY4DAm33377tGnT1q5d+/zzz//kJz9Je5xhxZtlhooHHnhg1apVDz/8cNqDkLLTTz/9pJNO6ujoyGQy99xzz7Jly5599tm0hyJlW7Zs+epXv7pu3bpjjjkm7VmGIS+Npqm3t/fHP/7xM888s3fv3r6+vkMPPTTtiRgSZsyY0X8wffr0H/7wh+kOw1DwwgsvnHjiiSpYJV4aTdMdd9zR2dn5zDPPrF+//rbbbtu3L4mPqGboKxaL/Qflcnn/59cTWS6XK5VKaU8xbPkZS9OmTZvOOOOMMWPGZDKZ5cuXpz0OQ8XTTz/df7By5copU6akOwxDwZQpU7q6ujZs2JD2IMOTl0bTdNFFF11++eXbt29//fXX9+zZk/Y4DBWvv/76hRdeOHr06BUrVjz11FNpj0P6mpubb7311nPOOWfmzJnFYrG3t3fx4sVpDzV8eLNMyjZs2LBmzZq2trbjjz/+7bffPuGEE9KeiJR1dXU1NTW9+uqr27ZtO+OMM8aNG5f2RAwVW7dufeGFF2pra88888zRo0enPc7wIYQAhOYeIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAoQkhAKEJIQChCSEAof0vqdPWO60qfIkAAAAASUVORK5CYII=", + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/html": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "execution_count": 83, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bar([\"a\", \"b\", \"c\"], [4,2,3])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0541b211", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/src/assets/competinghazards.png b/docs/src/assets/competinghazards.png new file mode 100644 index 0000000..81e76bf Binary files /dev/null and b/docs/src/assets/competinghazards.png differ diff --git a/docs/src/assets/conditional_on_when.png b/docs/src/assets/conditional_on_when.png new file mode 100644 index 0000000..af1648b Binary files /dev/null and b/docs/src/assets/conditional_on_when.png differ diff --git a/docs/src/assets/conditional_on_which.png b/docs/src/assets/conditional_on_which.png new file mode 100644 index 0000000..5734856 Binary files /dev/null and b/docs/src/assets/conditional_on_which.png differ diff --git a/docs/src/assets/gammahazard.png b/docs/src/assets/gammahazard.png new file mode 100644 index 0000000..7b59d21 Binary files /dev/null and b/docs/src/assets/gammahazard.png differ diff --git a/docs/src/assets/gammapdf.png b/docs/src/assets/gammapdf.png new file mode 100644 index 0000000..d9a54a2 Binary files /dev/null and b/docs/src/assets/gammapdf.png differ diff --git a/docs/src/assets/individual_distributions.png b/docs/src/assets/individual_distributions.png new file mode 100644 index 0000000..7bc6bdf Binary files /dev/null and b/docs/src/assets/individual_distributions.png differ diff --git a/docs/src/assets/shiftedgamma.png b/docs/src/assets/shiftedgamma.png new file mode 100644 index 0000000..9402c2a Binary files /dev/null and b/docs/src/assets/shiftedgamma.png differ diff --git a/docs/src/delay.md b/docs/src/delay.md index ee1bcc8..54e9b7e 100644 --- a/docs/src/delay.md +++ b/docs/src/delay.md @@ -60,18 +60,3 @@ So we need something more general to encompass these kinds of models. Counting p 2. [Leier, Andre, and Tatiana T. Marquez-Lago. "Delay chemical master equation: direct and closed-form solutions." Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471.2179 (2015): 20150049.](https://doi.org/10.1098/rspa.2015.0049) 3. [Brett, Tobias, and Tobias Galla. "Generating functionals and Gaussian approximations for interruptible delay reactions." Physical Review E 92.4 (2015): 042105.](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.92.042105) 4. [Nisbet, R. M., and W. S. C. Gurney. "The systematic formulation of population models for insects with dynamically varying instar duration." Theoretical Population Biology 23.1 (1983): 114-135.](https://www.sciencedirect.com/science/article/pii/0040580983900084) - - -## Attic - -Random things that may be useful to plunder for ideas: - - 1. https://github.com/slwu89/stochastic_ng - 2. https://github.com/slwu89/mosymodel/tree/master/premonition - 3. https://github.com/dd-harp/delay - 4. https://github.com/dd-harp/euler-abm - 5. https://github.com/dd-harp/mini-mash - 6. https://www.overleaf.com/project/5e3d95a0f8136c00010af62f - 7. https://www.overleaf.com/project/5bb3e39b92734f6c26b6fbd4 - 8. https://www.overleaf.com/project/5f249f60d9e45f0001f710cf - 9. https://www.overleaf.com/project/5de6ab381e175e0001c8c6bd \ No newline at end of file diff --git a/docs/src/distributions.md b/docs/src/distributions.md new file mode 100644 index 0000000..fbd357c --- /dev/null +++ b/docs/src/distributions.md @@ -0,0 +1,96 @@ +# How Fleck Decides What Happens Next + +Fleck simulates generalized semi-Markov processes (GSMP). Every event in a generalized semi-Markov process is chosen as the result of a competion among clocks to see which fires next. + +In a *process-oriented* simulation (like [SimJulia](https://simjuliajl.readthedocs.io/en/stable/welcome.html)), control flow is based on tasks. Each task performs some action on the state, rolls the dice, and sets a wake-up time. It might wake up as expected and possible execute code, or it might be interrupted by another task's actions. In contrast, an *event-oriented* simulation using Fleck will create a set of possible next events, assign a probability distribution for *when* each can happen, and the timing of which happens first determines *which* next event happens. Let's look at how a probability distribution describes the time for an event to happen and then how they compete in Fleck. + +## Distributions in Time + +Let's say you have a cold. You know you aren't going to recover immediately, but, as days go by, you're more and more sure you'll recover soon. You can read this graph below as the probability to recover given that you have not yet recovered. + +![](assets/gammahazard.png) + +It starts at zero, meaning there's no way you'll recover as soon as you're sick. It gets more likely over time that you're at the tail end of being sick. This is called a *hazard rate,* and the hazard rate shown is that of a Gamma distribution, commonly used to describe the rate of recovery for a population of individuals who are sick. + +If, instead, you want to see the number of people who recover on any given day, that is called a probability distribution function. + +![](assets/gammapdf.png) + +Where the hazard rate is an instantaneous quantity at a point in time, the probability distribution function (pdf) integrates over all possible future times. If we call the hazard rate $\lambda(t)$ and call the pdf $f$, we get this relationship. + +```math +f(t) = \lambda(t) e^{-\int_0^t \lambda(s)ds} +``` + +The graph of the pdf tells us that the most likely time for this event is a little before time 5, in whatever units. You will see graphs of [pdfs on Wikipedia](https://en.wikipedia.org/wiki/Gamma_distribution#/media/File:Gamma_distribution_pdf.svg) because this is how people usually think about the probability an event happens at some time. + +A simulation, however, has multiple events possible at any one time. One event may happen, and then other events need to restart. Let's ask, if you still have a cold on day 5, what is the probability distribution function for when you will recover? + +![](assets/shiftedgamma.png) + +The probability distribution function changes now that you know you didn't recover earlier than day 5. On the other hand, the hazard rate for recovery from the cold will be unchanged. Using the same hazard rate, we can recalculate the pdf from the new time $t_0=5$. + +```math +f(t;t>t_0) = \lambda(t) e^{-\int_{t_0}^t \lambda(s)ds} +``` + +In some sense, the underlying nature of when a particular event will happen is described more by the hazard rate, whereas the distribution function tells us about ensembles of events. + +The hazard rate is related to the well known cumulative distribution function (CDF) by an integral. The CDF tells us what is the overall probability the event occured some time in the interval ``[t_0,t_1)``. + +```math +F(t_0;t_1) = 1 - e^{-\int_{t_0}^{t_1} \lambda(s) ds} +``` + +Equally important for simulation is the survival function (sometimes called the complementary cumulative distribution function), which is the probability the event will not occur until after ``t_1``. + +```math +S(t_1) = 1 - F(t_1;t_1) = e^{-\int_{t_0}^{t_1} \lambda(s) ds} +``` + +## Competition + +### Individual Distributions + +Let's think of a moment when there are three possible next events. There is a Gamma distribution for when you recover from a cold, a Weibull distribution for when you decide to take medicine for the cold, and an Exponential distribution for when your Mom calls you. Each one is described by a distribution in time, and we can think of them as three hazard rates. + +![](assets/individual_distributions.png) + +The separate hazard rates are what we put into the simulation. Given their competition, the hazard rates will remain unchanged, but the pdfs will change. + +### Total Probability + +Each of the three clock distributions above corresponds to a unique event ``E_i``, which has a probability that it will be the first to fire. We calculate this probability by marginalizing over the other events, which ends up being an integral over the distribution, multiplied by the survivals of the other events. + +```math +P[E_i] = \int_0^\infty f_i(t) \prod_{j\ne i} S_j(t) dt +``` + +That gives the graph on the left. + +![](assets/conditional_on_which.png) + +The graph on the right shows the conditional distribution in time for each event, given that it was the one that fired, so it is $P[t_i | E_i]$. You can see that these distributions don't match the distributions for the individual events. They are modified by competition. + + +### Total Time + +From the graph above, if we pick a time, $t_1=10$, we can read from the graph three hazard rates, $(\lambda_1(t_1),\lambda_2(t_1), \lambda_3(t_1))$. Each hazard rate is the rate, per unit time, of that event. We know that, if the simulation makes it to $t=10$ without any event happening, the conditional probability for any one of those events is the ratio of hazard rates. + +```math +P[E_1|t=t_1] = \frac{\lambda_1(t_1)}{\lambda_1(t_1)+\lambda_2(t_1)+\lambda_3(t_1)} +``` + +If all distributions are Exponential, then they all have a constant hazard rate, and that conditional probability stays the same over time. If any distributions are non-Exponential, then that conditional probability changes. + +![](assets/conditional_on_when.png) + +On the left of this graph is the pdf for the first event of the three to fire. We can see this as a marginal $P[t]$ and then the right-hand graph as the conditional $P[E_i|t]$. + +## Specification of a Simulation + +We simulated using competing clocks when we measure using survival analysis. One is the inverse of the other. + +A classic example of survival analysis is a hospital study where participants enter a trial for a drug, and doctors observe the date at which each patient recovers. How would you draw a probability distribution function for such a trial? You would draw a histogram where each week shows a count of how many recovered. What about the hazard rate for such a trial? Here, you would look only at the number of participants remaining in the trial and ask what fraction recover each week. That's a rate of recovery given that participants have not yet recovered. + +There are many disciplines that use survival analysis to measure rates of processes. Epidemiology measures rates of contagion, death, and recovery. Demography measures time to have children or move to another country. Reliability studies the time for a part to break and be fixed. Drug discovery measures the time to approve a chemical for the next trial stage. For traffic flow, it's the time to the next intersection. For package delivery, it's the time to the next package scan. All of these do use, or can use, survival analysis to observe a system to simulate. diff --git a/docs/src/guide.md b/docs/src/guide.md new file mode 100644 index 0000000..3fb98b7 --- /dev/null +++ b/docs/src/guide.md @@ -0,0 +1,17 @@ +# Bag of Clocks + +We can introduce this style of simulation with a description that has no equations. + +We use the term discrete event simulation when code tracks the effects of individual events in time. We classify discrete event simulations by the pace of advances in time. If time advances in steps of the same size, it is a *discrete-time* simulation. If time advances to the next event, whenever that event might happen, it is *continuous-time* simulation. This library supports continuous-time discrete-event simulations. + +We write discrete event simulations by defining the state of a system and events that can act on that state. The state of a system can be the count of each chemical species for a chemical reaction system. It can be different numbers of chickens at different weights in different chicken coops for a farming simulation. It can be the age of each widget in an engine for a reliability simulation. These states are called the *physical state* of the system and vary widely in what they represent. + +Even though nearly every simulation represents a unique system, it is often possible to find enough in common among systems that we can make libraries to help create, simulate, and observe changes of state. For instance, a large class of physical states are counts of things, whether they are molecules, vehicles, or chickens. These states can be represented in code as vectors of integers, and there is a long history of models represented by vectors of integers. Most of this work focuses on using chemical simulation codes to represent not only chemicals but also other systems like vehicles. An older, more general, form of this simulation is the vector addition system, which a later section of this manual describes. However, states don't need to be vectors of integers. + +We usually want to simulate systems that have real-valued observations such as ages, temperatures, and distances. These continuous parts of the state of the system mean that the system rarely returns to exactly the same state. There are now an infinite number of possible states of the system, which changes how we would think of the statistics of the system, but it doesn't much change how we handle simulation, because simulation is still driven by events. + +An *event* is an observable change of state at a time. It is a function where the argument is the current state of the system and the return value is the new state of the system. We say an event occurs, happens, or fires when we apply that function to the state. How that state changes is up to the model and how it defines the state, but what about *when* that state changes? + +Given the current state of the system, there is a set of possible next events which could happen. We call the possible next events *enabled.* Think of the moment in time just after the simulation has begun or an event has just happened. At this moment, there may be multiple enabled events. Which one is next depends on how often that event happens. If our simulation includes both the radioactive decay of an unstable element and decay of iron, the unstable element will usually decay first. We describe the rates of these events using probability distributions in time. Each event has an associated continuous univariate distribution where the variable is time. + +We can think of all event distributions in a model as a *bag of clocks.* The next one to ring is the next event to occur. When that event occurs, the state changes. When the state changes, it will enable some events and disable others. Enabled events are added to the bag of clocks. Disabled ones are removed. This library holds the bag of clocks. diff --git a/docs/src/index.md b/docs/src/index.md index 84bd371..30145be 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,8 +18,7 @@ If you want to write a simulation engine that tracks the state of a system, and * Generalized stochastic Petri nets. * Generalized semi-Markov Processes. -In statistical terms, this library is a sampler for generalized semi-Markov processes. - +Many of these mathematical objects provide intricate structure to tell us what events are enabled in any particular state (e.g; Petri nets). Fleck organizes the stochastic firing times (clocks) of enabled events so we can sample from clocks to find out what happens next and update them after the state changes. In statistical terms, this library is a sampler for generalized semi-Markov processes. ## Usage @@ -41,3 +40,5 @@ If I make a quick simulation for myself, I sample distributions the moment an ev * I'm looking at rare events, so I want to use splitting techniques and [importance sampling](https://en.wikipedia.org/wiki/Importance_sampling). * Performance matters (which it often doesn't), so I would like to try different samplers on my problem. + + * I want to focus on developing and testing my *model* not my *simulation algorithm*; Fleck is designed and tested with care to ensure correctness. \ No newline at end of file diff --git a/docs/src/mainloop.md b/docs/src/mainloop.md new file mode 100644 index 0000000..7bbed84 --- /dev/null +++ b/docs/src/mainloop.md @@ -0,0 +1,191 @@ +```@meta +EditURL = "simple_board.jl" +``` + +# Sample Main Loop + +## Introduction +Let's walk through a small simulation so that we can see how Fleck +could appear in the main loop. +This models individuals wandering around on a checkerboard, where no +two individuals can occupy the same space. First, let's import libraries. + +````julia +import Base +using Distributions +using Random +using SparseArrays +using Test +using Fleck +```` + +There is a physical state for the model, separate from the state of +the bag of clocks and separate from the time. Most of the board is empty, +so let's use a spare matrix to represent the locations of individuals. + +````julia +mutable struct PhysicalState + board::SparseMatrixCSC{Int64, Int64} +end + +# They can move in any of four directions. +@enum Direction Up Left Down Right +const DirectionDelta = Dict( + Up => CartesianIndex(-1, 0), + Left => CartesianIndex(0, -1), + Down => CartesianIndex(1, 0), + Right => CartesianIndex(0, 1), + ); +```` + +The simulation, itself, carries the state of the clocks in the sampler, as +well as the physical state. We'll call it a finite state machine (FSM) +because it has the traits of a +[Moore machine](https://en.wikipedia.org/wiki/Moore_machine). + +````julia +mutable struct SimulationFSM{Sampler} + physical::PhysicalState + sampler::Sampler + when::Float64 + rng::Xoshiro +end +```` + +## Main Loop +The main loop will ask what event happens next to the state. When that +event changes the state, the loop will update the set of possible next +events by disabling outdated ones and enabling new ones. The calls to +Fleck are: + +* `next(sampler, current time, random number generator (RNG))` +* `enable!(sampler, event ID, distribution, current time, start time of clock, RNG)` +* `disable!(sampler, event ID, current time)` + +````julia +function run(event_count) + # An event is (who moves, where they start, which direction) + EventKey = Tuple{Int,CartesianIndex{2},Direction} + Sampler = CombinedNextReaction{EventKey,Float64} + physical = PhysicalState(zeros(Int, 10, 10)) + @test showable(MIME("text/plain"), physical) + sim = SimulationFSM{Sampler}( + physical, + Sampler(), + 0.0, + Xoshiro(2947223) + ) + initialize!(sim.physical, 9, sim.rng) + current_events = allowed_moves(sim.physical) + for event_id in current_events + enable!(sim.sampler, event_id, Weibull(1.0), 0.0, 0.0, sim.rng) + end + + for i in 1:event_count + (when, what) = next(sim.sampler, sim.when, sim.rng) + if isfinite(when) && !isnothing(what) + sim.when = when + move!(sim.physical, what) + next_events = allowed_moves(sim.physical) + for remove_event in setdiff(current_events, next_events) + disable!(sim.sampler, remove_event, when) + end + for add_event in setdiff(next_events, current_events) + enable!(sim.sampler, add_event, Weibull(1.0), when, when, sim.rng) + end + current_events = next_events + @show (when, what) + end + end +end +```` + +```` +run (generic function with 1 method) +```` + +For this checkerboard with wandering individuals, the allowed moves are +moves by any individual to any free square on the board. We will denote +each allowed move by a tuple, (who moves, where they are, which direction). + +````julia +function allowed_moves(physical::PhysicalState) + allowed = Set{Tuple{Int,CartesianIndex{2},Direction}}() + row, col, value = findnz(physical.board) + for ind_idx in eachindex(value) + location = CartesianIndex((row[ind_idx], col[ind_idx])) + for (direction, offset) in DirectionDelta + if checkbounds(Bool, physical.board, location + offset) + if physical.board[location + offset] == 0 + push!(allowed, (value[ind_idx], location, direction)) + end + end + end + end + return allowed +end +```` + +```` +allowed_moves (generic function with 1 method) +```` + +## Changes to the state of the board +We set up the board with an initializer that places individuals at random. +We move one individual at a time, when their next event happens. + +````julia +function initialize!(physical::PhysicalState, individuals::Int, rng) + physical.board .= 0 + dropzeros!(physical.board) + locations = zeros(CartesianIndex{2}, individuals) + for ind_idx in 1:individuals + loc = rand(rng, CartesianIndices(physical.board)) + while physical.board[loc] != 0 + loc = rand(rng, CartesianIndices(physical.board)) + end + locations[ind_idx] = loc + physical.board[loc] = ind_idx + end +end + + +function move!(physical::PhysicalState, event_id) + (individual, previous_location, direction) = event_id + next_location = previous_location + DirectionDelta[direction] + # This sets the previous board value to zero. + SparseArrays.dropstored!(physical.board, previous_location.I...) + physical.board[next_location] = individual +end +```` + +```` +move! (generic function with 1 method) +```` + +## How it runs +A run of this simulation produces a sequence of moves, no two happening +at the same time because this is a continuous-time simulation. + +````julia +run(10) +```` + +```` +(when, what) = (0.07319364933011555, (7, CartesianIndex(7, 6), Main.var"##225".Left)) +(when, what) = (0.10866577949385807, (9, CartesianIndex(1, 3), Main.var"##225".Right)) +(when, what) = (0.15079187330778177, (9, CartesianIndex(1, 4), Main.var"##225".Left)) +(when, what) = (0.1682581650543986, (9, CartesianIndex(1, 3), Main.var"##225".Right)) +(when, what) = (0.16874877793434204, (3, CartesianIndex(1, 8), Main.var"##225".Left)) +(when, what) = (0.18888751327810988, (3, CartesianIndex(1, 7), Main.var"##225".Right)) +(when, what) = (0.20347320994043128, (4, CartesianIndex(9, 1), Main.var"##225".Down)) +(when, what) = (0.2118081725791219, (1, CartesianIndex(7, 10), Main.var"##225".Down)) +(when, what) = (0.2152994448757914, (3, CartesianIndex(1, 8), Main.var"##225".Right)) +(when, what) = (0.2634608206635203, (8, CartesianIndex(2, 10), Main.var"##225".Down)) + +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/objects.md b/docs/src/objects.md index c9380c8..00f2fc7 100644 --- a/docs/src/objects.md +++ b/docs/src/objects.md @@ -19,7 +19,7 @@ For a simulation, we think of the distribution as being over time, and it's over In Julia, the [univariate distributions](https://juliastats.org/Distributions.jl/stable/univariate/) don't carry their enabling time as a parameter, so we store it separately. -## Competing Clock +## Competing Clocks Let's say you make a simulation of rabbits eating kibble. A rabbit just picked up a piece of kibble, and the simulation decides there is a Gamma-distributed time at which it will be ready to eat the next bit of kibble. Meanwhile, another rabbit eats the last bit of kibble. That makes it impossible for the first rabbit to eat any kibble, so we say that its event is interrupted. That interruption was a result of the rabbits competing, but what the simulation sees is distributions competing to fire first, and we call them competing clocks. diff --git a/docs/src/reference.md b/docs/src/reference.md index 08e925d..ce82b02 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -4,9 +4,70 @@ CurrentModule = Fleck # Reference -```@index + * [Interface](@ref) + * [Samplers](@ref) + * [Algorithms](@ref) + + +## Interface + + +### Use a Sampler + +```@docs +enable! +disable! +next ``` -```@autodocs -Modules = [Fleck] +### Query a Sampler + +```@docs +getindex +keys +keytype +length +``` + +## Samplers + +### Sampler Supertype + +```@docs +SSA +``` + +### Sampler Types + +```@docs +FirstReaction +FirstToFire +DirectCall +CombinedNextReaction +consume_survival +sampling_space +``` + +### Sampling Helpers + +```@docs +ChatReaction +CommonRandomRecorder +DebugWatcher +MultiSampler +SingleSampler +TrackWatcher +``` + +## Algorithms + +```@docs +CumSumPrefixSearch +BinaryTreePrefixSearch +KeyedKeepPrefixSearch +KeyedRemovalPrefixSearch +choose +setindex! +rand +set_multiple! ``` diff --git a/docs/src/rules.md b/docs/src/rules.md new file mode 100644 index 0000000..95f018d --- /dev/null +++ b/docs/src/rules.md @@ -0,0 +1,7 @@ +# Rules and Guidelines + +Even though this library is just handling a bag of event clocks, there are some rules for how to call the interface so that those clocks remain consistent. There are rules that keep the simulation running, and there are rules that guarantee repeatable simulations that have the correct statistical likelihood. + + * Some simulations treat events as maybe changing state, maybe not changing state. + * Some simulations think of an event as probabilistically changing state. + * For a simulation in continuous time, only one event can happen at any time. diff --git a/docs/src/simple_board.jl b/docs/src/simple_board.jl new file mode 100644 index 0000000..e083bd7 --- /dev/null +++ b/docs/src/simple_board.jl @@ -0,0 +1,148 @@ +# # Sample Main Loop + +# ## Introduction +# Let's walk through a small simulation so that we can see how Fleck +# could appear in the main loop. +# This models individuals wandering around on a checkerboard, where no +# two individuals can occupy the same space. First, let's import libraries. + +import Base +using Distributions +using Random +using SparseArrays +using Test +using Fleck + +# There is a physical state for the model, separate from the state of +# the bag of clocks and separate from the time. Most of the board is empty, +# so let's use a spare matrix to represent the locations of individuals. + +mutable struct PhysicalState + board::SparseMatrixCSC{Int64, Int64} +end + +## They can move in any of four directions. +@enum Direction Up Left Down Right +const DirectionDelta = Dict( + Up => CartesianIndex(-1, 0), + Left => CartesianIndex(0, -1), + Down => CartesianIndex(1, 0), + Right => CartesianIndex(0, 1), + ); + +# The simulation, itself, carries the state of the clocks in the sampler, as +# well as the physical state. We'll call it a finite state machine (FSM) +# because it has the traits of a +# [Moore machine](https://en.wikipedia.org/wiki/Moore_machine). + +mutable struct SimulationFSM{Sampler} + physical::PhysicalState + sampler::Sampler + when::Float64 + rng::Xoshiro +end + +# ## Main Loop +# The main loop will ask what event happens next to the state. When that +# event changes the state, the loop will update the set of possible next +# events by disabling outdated ones and enabling new ones. The calls to +# Fleck are: +# +# * `next(sampler, current time, random number generator (RNG))` +# * `enable!(sampler, event ID, distribution, current time, start time of clock, RNG)` +# * `disable!(sampler, event ID, current time)` +# +# There are a lot of samplers in Fleck to choose from. This example uses `CombinedNextReaction` +# algorithm, which has good performance for a variety of distributions. Samplers in Fleck +# require two type parameters, a key type for clocks and the type used to represent time. +# In this case, the clock key type fully represents an event, giving the ID of the individual, +# where they start, and which direction they may move. + +const ClockKey = Tuple{Int,CartesianIndex{2},Direction} + +function run(event_count) + Sampler = CombinedNextReaction{ClockKey,Float64} + physical = PhysicalState(zeros(Int, 10, 10)) + @test showable(MIME("text/plain"), physical) + sim = SimulationFSM{Sampler}( + physical, + Sampler(), + 0.0, + Xoshiro(2947223) + ) + initialize!(sim.physical, 9, sim.rng) + current_events = allowed_moves(sim.physical) + for event_id in current_events + enable!(sim.sampler, event_id, Weibull(1.0), 0.0, 0.0, sim.rng) + end + + for i in 1:event_count + (when, what) = next(sim.sampler, sim.when, sim.rng) + if isfinite(when) && !isnothing(what) + sim.when = when + move!(sim.physical, what) + next_events = allowed_moves(sim.physical) + for remove_event in setdiff(current_events, next_events) + disable!(sim.sampler, remove_event, when) + end + for add_event in setdiff(next_events, current_events) + enable!(sim.sampler, add_event, Weibull(1.0), when, when, sim.rng) + end + current_events = next_events + @show (when, what) + end + end +end; + +# For this checkerboard with wandering individuals, the allowed moves are +# moves by any individual to any free square on the board. The set of allowed +# moves is precisely the set of enabled clocks, so it stores `ClockKey`s. + +function allowed_moves(physical::PhysicalState) + allowed = Set{ClockKey}() + row, col, value = findnz(physical.board) + for ind_idx in eachindex(value) + location = CartesianIndex((row[ind_idx], col[ind_idx])) + for (direction, offset) in DirectionDelta + if checkbounds(Bool, physical.board, location + offset) + if physical.board[location + offset] == 0 + push!(allowed, (value[ind_idx], location, direction)) + end + end + end + end + return allowed +end; + +# ## Changes to the state of the board +# We set up the board with an initializer that places individuals at random. +# We move one individual at a time, when their next event happens. + +function initialize!(physical::PhysicalState, individuals::Int, rng) + physical.board .= 0 + dropzeros!(physical.board) + locations = zeros(CartesianIndex{2}, individuals) + for ind_idx in 1:individuals + loc = rand(rng, CartesianIndices(physical.board)) + while physical.board[loc] != 0 + loc = rand(rng, CartesianIndices(physical.board)) + end + locations[ind_idx] = loc + physical.board[loc] = ind_idx + end +end; + + +function move!(physical::PhysicalState, event_id) + (individual, previous_location, direction) = event_id + next_location = previous_location + DirectionDelta[direction] + ## This sets the previous board value to zero. + SparseArrays.dropstored!(physical.board, previous_location.I...) + physical.board[next_location] = individual +end; + +# ## How it runs +# A run of this simulation produces a sequence of moves, no two happening +# at the same time because this is a continuous-time simulation. + +run(10) diff --git a/docs/src/vas.md b/docs/src/vas.md index 880a2f0..10df7c2 100644 --- a/docs/src/vas.md +++ b/docs/src/vas.md @@ -1,29 +1,28 @@ # Vector Addition System -A [vector addition system](https://en.wikipedia.org/wiki/Vector_addition_system) (VAS) is a simple mathematical language to model systems. It is simpler than the Petri net language, but has the advantage of allowing for more free-form simulation design, as intensity functions can be arbitrary functions of state. +One way to understand how to use Fleck is to look at a very simple simulation. A [vector addition system](https://en.wikipedia.org/wiki/Vector_addition_system) (VAS) is a lot like chemical simulations, but it's free of some of the assumptions about chemical rates. The physical state of the system is a vector of integers. Such a simple physical state can make it easier to understand possible complications in how we define events. -```mermaid -classDiagram +There are variations on how to define a vector addition system. Let's begin with one and introduce variations as they become interesting. -class VectorAdditionSystem - VectorAdditionSystem: +Matrix take - VectorAdditionSystem: +Matrix give - VectorAdditionSystem: +Vector rates +A vector addition system is + * a physical state $p$ which is a set of $d$ integers, labeled $(p_1, p_2,\cdots, p_d)$. -class VectorAdditionModel - VectorAdditionModel: +VectorAdditionSystem vas - VectorAdditionModel: +Vector state - VectorAdditionModel: +Float64 when + * a system time, $t$ + + * a finite set of events $E$ where each event has an enabling rule, a distribution, and a firing rule. + + - a rule for when the event is enabled, where the rule is an invariant on the physical state, such as $r_i(p)\ge 0$. When the invariant is met, the event is enabled. When the invariant is not met, the event is disabled. For a vector addition system, this rule must be expressible as a matrix multiplication, so $M\cdot \vec{p} .\ge 0$, where the $.\ge$ means we are checking elementwise that every element is non-negative. + + - a distribution in time, where the distribution is determined at enabling time and is a function of the physical state and system time at enabling. The distribution in time can have shocks, such as a delta function, but it may not have a shock at time zero, must have a measure-zero probability of firing at the same moment it is enabled. + + - a rule for how the event changes physical state when it happens. The domain and co-domain of the rule are physical states of the system. Let's limit this model to require all events to modify the state, so that we exclude events that don't modify the state. As with enabling, this must be expressible as a matrix operation on the physical states to produce a new physical state. + + From the rules, we see that an event's state is either disabled or enabled at an enabling time $t_e$. + +How would we implement this system? It's set up so that enabling rules and firing rules can be done with linear algebra. We can write this in Julia code -class VectorAdditionFSM - VectorAdditionFSM: +VectorAdditionModel vam - VectorAdditionFSM: +Any sampler -end -VectorAdditionModel *-- VectorAdditionSystem -VectorAdditionFSM *-- VectorAdditionModel -``` Simulation using the VAS requires concrete implementation of the interface: diff --git a/docs/src/walk_board.jl b/docs/src/walk_board.jl new file mode 100644 index 0000000..9460277 --- /dev/null +++ b/docs/src/walk_board.jl @@ -0,0 +1,158 @@ +# # Sample Main Loop + +# Let's walk through a small simulation so that we can see how and when Fleck +# could appear in the main loop. We're going to split this up here in a +# way that is more simplistic than the Fleck interface requires. + +import Base +using Distributions +using Random +using Test +using Fleck + +# There is a physical state for the model, separate from the state of +# the bag of clocks and separate from the time. + +mutable struct PhysicalState + board::Matrix{Int} +end + +function Base.show(io::IO, physical::PhysicalState) + for row_idx in 1:size(physical.board, 1) + println(io, join([string(x) for x in physical.board[row_idx, :]], "")) + end +end + + +@enum Direction Up Left Down Right +const DirectionDelta = Dict( + Up => CartesianIndex(-1, 0), + Left => CartesianIndex(0, -1), + Down => CartesianIndex(1, 0), + Right => CartesianIndex(0, 1), + ) +const OppositeDirection = Dict( + Up => Down, + Down => Up, + Left => Right, + Right => Left +) + +function over_neighbors(f::Function, physical::PhysicalState, location) + for (direction, offset) in DirectionDelta + if checkbounds(Bool, physical.board, location + offset) + f(direction, location + offset) + end + end +end + + +## Use separate initializer in case you want to re-initialize the same memory. +function initialize!(physical::PhysicalState, sampler, individuals::Int, rng) + physical.board .= 0 + locations = zeros(CartesianIndex{2}, individuals) + for ind_idx in 1:individuals + loc = rand(rng, CartesianIndices(physical.board)) + while physical.board[loc] != 0 + loc = rand(rng, CartesianIndices(physical.board)) + end + locations[ind_idx] = loc + physical.board[loc] = ind_idx + end + # During initialization, tell the sampler each individual can move to + # any blank spot next to them. + for mover_idx in 1:individuals + location = locations[mover_idx] + over_neighbors(physical, location) do direction, neighbor_location + if physical.board[neighbor_location] == 0 + rate = Weibull(1.0) + enable!(sampler, (mover_idx, direction), rate, 0.0, 0.0, rng) + end + end + end +end + + +function move!(physical::PhysicalState, sampler, event_id, when, rng) + (individual, direction) = event_id + previous_location = findfirst(x -> x == individual, physical.board) + next_location = previous_location + DirectionDelta[direction] + + # The individual who moves will disable all of its previous jumps. + over_neighbors(physical, previous_location) do direction, neighbor_location + if physical.board[neighbor_location] == 0 + disable!(sampler, (individual, direction), when) + end + end + # If that individual had neighbors, those neighbors are now free to move + # to the newly-vacated spot. + over_neighbors(physical, previous_location) do direction, neighbor_location + previous_neighbor = physical.board[neighbor_location] + if previous_neighbor != 0 + rate = Weibull(1.0) + now_move = OppositeDirection[direction] + enable!(sampler, (previous_neighbor, now_move), rate, when, when, rng) + end + end + + physical.board[previous_location] = 0 + physical.board[next_location] = individual + + # The individual moved to a new spot that may block neighbor movements. + over_neighbors(physical, next_location) do direction, neighbor_location + new_neighbor = physical.board[neighbor_location] + if new_neighbor != 0 + no_longer_move = OppositeDirection[direction] + disable!(sampler, (new_neighbor, no_longer_move), when) + end + end + # The individual now has a new set of jumps it can do, depending on + # free spots nearby. + over_neighbors(physical, next_location) do direction, neighbor_location + if physical.board[neighbor_location] == 0 + rate = Weibull(1.0) + enable!(sampler, (individual, direction), rate, when, when, rng) + end + end +end + + +# The simulation, itself, carries the state of the clocks, as well as the +# physical state. We'll call it a finite state machine (FSM) because it has +# the traits of a [Moore machine](https://en.wikipedia.org/wiki/Moore_machine). + +mutable struct SimulationFSM{Sampler} + physical::PhysicalState + sampler::Sampler + when::Float64 + rng::Xoshiro +end + + +function run() + EventKey = Tuple{Int,Direction} + TimeType = Float64 + Sampler = CombinedNextReaction{EventKey,TimeType} + physical = PhysicalState(zeros(Int, 10, 10)) + @test showable(MIME("text/plain"), physical) + sim = SimulationFSM{Sampler}( + physical, + Sampler(), + 0.0, + Xoshiro(2947223) + ) + initialize!(sim.physical, sim.sampler, 9, sim.rng) + for i in 1:100 + (when, what) = next(sim.sampler, sim.when, sim.rng) + if isfinite(when) + sim.when = when + move!(sim.physical, sim.sampler, what, sim.when, sim.rng) + @show (when, what) + else + println("no possible moves") + break + end + end +end + +run() diff --git a/src/prefixsearch/binarytreeprefixsearch.jl b/src/prefixsearch/binarytreeprefixsearch.jl index 653780c..5ac0b22 100644 --- a/src/prefixsearch/binarytreeprefixsearch.jl +++ b/src/prefixsearch/binarytreeprefixsearch.jl @@ -4,6 +4,9 @@ using Distributions: Uniform using Logging """ + BinaryTreePrefixSearch{T}(N=32) + +This stores hazard rates to make them faster for the Direct method to sample. This is a binary tree where the leaves are values and the nodes are sums of those values. It is meant to make it easier to find the leaf such that the sum of it and all previous @@ -35,11 +38,7 @@ time_type(ps::BinaryTreePrefixSearch{T}) where {T} = T time_type(::Type{BinaryTreePrefixSearch{T}}) where {T} = T -""" - ceil_log2(v::Integer) - -Integer log2, rounding up. -""" +# Integer log2, rounding up. function ceil_log2(v::Integer) r = 0 power_of_two = ((v & (v - 1)) == 0) ? 0 : 1 diff --git a/src/prefixsearch/cumsumprefixsearch.jl b/src/prefixsearch/cumsumprefixsearch.jl index 28f3428..02bc1ec 100644 --- a/src/prefixsearch/cumsumprefixsearch.jl +++ b/src/prefixsearch/cumsumprefixsearch.jl @@ -3,6 +3,15 @@ using Random using Distributions: Uniform +""" + CumSumPrefixSearch{T}() + +This stores hazard rates in order to make it easier for the Direct +method to sample them. This version is the dumbest possible, but it can +be faster when there are few hazards enabled. It uses a simple array +and, each time the Direct method samples, this evaluates the cumulative +sum of the array. +""" struct CumSumPrefixSearch{T<:Real} array::Vector{T} cumulant::Vector{T} diff --git a/src/sample/combinednr.jl b/src/sample/combinednr.jl index 4c2e910..2c3accf 100644 --- a/src/sample/combinednr.jl +++ b/src/sample/combinednr.jl @@ -99,17 +99,21 @@ end This combines Next Reaction Method and Modified Next Reaction Method. The Next Reaction Method is from Gibson and Bruck in their 2000 paper called -``Efficient Exact Stochastic Simulation of Chemical Systems with Many Species +"Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels." The Modified Next Reaction Method is from David F. Anderson's -2007 paper, ``A modified Next Reaction Method for simulating chemical systems +2007 paper, "A modified Next Reaction Method for simulating chemical systems with time dependent propensities and delays." Both methods reuse draws of random numbers. The former works by accumulating survival of a distribution in a linear space and the latter works by accumulating survival of a distribution in a log space. -Each distribution is more precise being sampled in either a linear space -or a log space. This sampler chooses which space to use depending on the -type of the `UnivariateDistribution`. Defaults are set for those distributions +Each enabled clock specifies a univariate distribution from the `Distributions` +package. Every distribution is more precise being sampled in the manner +of the Next Reaction method (linear space) or the manner of the Modified +Next Reaction method (log space). This sampler +chooses which space to use depending on the +type of the `UnivariateDistribution` and based on performance timings that +are done during package testing. Defaults are set for those distributions included in the `Distributions.jl` package. If you want to add a distribution, then define: diff --git a/src/sample/direct.jl b/src/sample/direct.jl index 187e70a..ded305d 100644 --- a/src/sample/direct.jl +++ b/src/sample/direct.jl @@ -6,26 +6,40 @@ export DirectCall, enable!, disable!, next """ - DirectCall{K,T} + DirectCall{KeyType,TimeType}() DirectCall is responsible for sampling among Exponential distributions. It samples using the Direct method. In this case, there is no optimization to that Direct method, so we call it DirectCall because it recalculates everything every time you call it. -The type `K` is the type of an identifier for each transition. This identifier +The type `KeyType` is the type of an identifier for each transition. This identifier is usually a nominal integer but can be a any key that identifies it, such as -a string or tuple of integers. Instances of type `K` are used as keys in a -dictionary. +a string or tuple of integers. Instances of type `KeyType` are used as keys in a +dictionary. The type `TimeType` is usually `Float64`. + +The algorithm for the Direct Method relies heavily on what data structure +it uses to maintain a list of hazard rates, such that it can know the sum +of those hazards and index into them using a random value. This struct +has a default constructor that chooses a data structure for you, but there +are several options. # Example ```julia -DirectCall{K,T}() where {K,T} = - DirectCall{K,CumSumPrefixSearch{T}}(CumSumPrefixSearch(T)) +sampler_string_keys = DirectCall{String,Float64}() +sampler_inttuple_keys = DirectCall{(Int,Int), Float64}() +``` +If we know that our simulation will only use a small number of different +clock keys, then it would make sense to use a data structure that disables +clocks by zeroing them out, instead of removing them from the list. This +will greatly reduce memory churn. We can do that by changing the +underlying data structure. -DirectCall{K,T}() where {K,T} = - DirectCall{K,BinaryTreePrefixSearch{T}}(BinaryTreePrefixSearch(T)) +```julia +prefix_tree = BinaryTreePrefixSearch{T}() +keyed_prefix_tree = KeyedKeepPrefixSearch{K,typeof(prefix_tree)}(prefix_tree) +sampler_noremove = DirectCall{K,T,typeof(keyed_prefix_tree)}(keyed_prefix_tree) ``` """ diff --git a/src/sample/firstreaction.jl b/src/sample/firstreaction.jl index 22ca99a..0cddfd0 100644 --- a/src/sample/firstreaction.jl +++ b/src/sample/firstreaction.jl @@ -6,9 +6,19 @@ using Logging export FirstReaction, ChatReaction """ + FirstReaction{KeyType,TimeType}() + +This classic, first reaction method is the simplest and most assuredly correct +sampler. It can also be very fast when there are only a few clocks to sample. + Classic First Reaction method for Exponential and non-Exponential -distributions. Every time you sample, go to each distribution and ask when it -would fire. Then take the soonest and throw out the rest until the next sample. +distributions. Every time you sample, this goes to each distribution and asks +when it would fire. Then it takes the soonest and throws out the rest of the +sampled times until the next sample. + +One interesting property of this sampler is that you can call `next()` +multiple times in order to get a distribution of next firing clocks and their +times to fire. """ struct FirstReaction{K,T} <: SSA{K,T} # This other class already stores the current set of distributions, so use it. @@ -70,7 +80,8 @@ end """ This sampler can help if it's the first time you're trying a model. It checks -all of the things and uses Julia's logger to communicate them. +all of the things and uses Julia's logger to communicate them. It samples +using the first reaction algorithm. """ mutable struct ChatReaction{K,T} # This other class already stores the current set of distributions, so use it. diff --git a/src/sample/firsttofire.jl b/src/sample/firsttofire.jl index 64c5c85..1d4cdff 100644 --- a/src/sample/firsttofire.jl +++ b/src/sample/firsttofire.jl @@ -5,10 +5,11 @@ export FirstToFire """ FirstToFire{KeyType,TimeType}() -Construct a FirstToFire. -As soon as a distribution is enabled, this draws a value from the distribution. -The soonest to fire wins. When a clock is disabled, its future firing time is -removed from the list. There is no memory of previous firing times. +This sampler is often the fastest for non-exponential distributions. +When a clock is first enabled, this sampler asks the clock when it would +fire and saves that time in a sorted heap of future times. Then it works +through the heap, one by one. When a clock is disabled, its future firing time +is removed from the list. There is no memory of previous firing times. """ struct FirstToFire{K,T} <: SSA{K,T} firing_queue::MutableBinaryMinHeap{OrderedSample{K,T}} diff --git a/src/sample/interface.jl b/src/sample/interface.jl index 343e82a..5c138e2 100644 --- a/src/sample/interface.jl +++ b/src/sample/interface.jl @@ -6,6 +6,19 @@ export enable!, disable!, next, getindex, keys, length, keytype +""" + enable!(sampler, clock, distribution, enablingtime, currenttime, RNG) + +Tell the sampler to start a clock. + + * `sampler::SSA{KeyType,TimeType}` - The sampler to tell. + * `clock::KeyType` - The ID of the clock. Can be a string, integer, tuple, etc. + * `distribution::Distributions.UnivariateDistribution` + * `enablingtime::TimeType` - The zero time for the clock's distribution, in absolute time. Usually equal to `when`. + * `when::TimeType` - The current time of the simulation. + * `rng::AbstractRNG` - A random number generator. + +""" function enable!( sampler::SSA{K,T}, clock::K, @@ -16,27 +29,47 @@ function enable!( ) where {K,T} end +""" + disable!(sampler, clock, when) + +Tell the sampler to forget a clock. We include the current simulation time +because some Next Reaction methods use this to optimize sampling. +""" function disable!(sampler::SSA{K,T}, clock::K, when::T) where {K,T} end +""" + next(sampler, when, rng) + +Ask the sampler for what happens next, in the form of +`(when, which)::Tuple{TimeType,KeyType}`. `rng` is a random number generator. +""" function next(sampler::SSA{K,T}, when::T, rng::AbstractRNG) where {K,T} end """ - Return stored state for a particular clock. If the clock does not exist, + getindex(sampler, clock::KeyType) + +Return stored state for a particular clock. If the clock does not exist, a `KeyError` will be thrown. """ function Base.getindex(sampler::SSA{K,T}, clock::K) where {K,T} end """ - Return all stored clocks as a vector. + keys(sampler) + +Return all stored clocks as a vector. """ function Base.keys(sampler::SSA) end """ - Return the number of stored clocks. + length(sampler)::Int64 + +Return the number of stored clocks. """ function Base.length(sampler::SSA) end """ - Return the type of clock keys. + keytype(sampler) + +Return the type of clock keys. """ Base.keytype(::SSA{K,T}) where {K,T} = K diff --git a/src/sample/nrtransition.jl b/src/sample/nrtransition.jl index f46d32a..17eddf0 100644 --- a/src/sample/nrtransition.jl +++ b/src/sample/nrtransition.jl @@ -1,9 +1,7 @@ import Base: ==, <, > -""" -A record of a transition and the time. -It's sortable by time. Immutable. -""" +# A record of a transition and the time. +# It's sortable by time. Immutable. struct OrderedSample{K,T} key::K time::T diff --git a/src/sample/sampler.jl b/src/sample/sampler.jl index f3211df..9073526 100644 --- a/src/sample/sampler.jl +++ b/src/sample/sampler.jl @@ -3,11 +3,15 @@ export SingleSampler, MultiSampler export SSA, enable!, disable!, sample! """ - SSA{Key,Time} + SSA{KeyType,TimeType} This abstract type represents a stochastic simulation algorithm (SSA). It is parametrized by the clock ID, or key, and the type used for the time, which -is typically a Float64. +is typically a Float64. The type of the key can be anything you would use +as a dictionary key. This excludes mutable values but includes a wide range +of identifiers useful for simulation. For instance, it could be a `String`, +but it could be a `Tuple{Int64,Int64,Int64}`, so that it indexes into a +complicated simulation state. """ abstract type SSA{Key,Time} end