diff --git a/.gitignore b/.gitignore index 1c262883..51b491e1 100644 --- a/.gitignore +++ b/.gitignore @@ -5,7 +5,8 @@ /Manifest.toml /docs/build/ /docs/Manifest.toml -/benchmark/HMMBenchmark/Manifest.toml +/test/Manifest.toml +/libs/**/Manifest.toml scratchpad.jl diff --git a/Project.toml b/Project.toml index 7ce2aa81..85591c30 100644 --- a/Project.toml +++ b/Project.toml @@ -13,18 +13,17 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" [weakdeps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" HMMBase = "b2b3ca75-8444-5ffa-85e6-af70e2b64fe7" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [extensions] HiddenMarkovModelsDistributionsExt = "Distributions" HiddenMarkovModelsHMMBaseExt = "HMMBase" -HiddenMarkovModelsTestExt = "Test" +HiddenMarkovModelsSparseArraysExt = "SparseArrays" [compat] ChainRulesCore = "1.16" @@ -40,43 +39,5 @@ Requires = "1.3" SimpleUnPack = "1.1" SparseArrays = "1" StatsAPI = "1.6" -Test = "1" julia = "1.6" -[extras] -Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -HMMBase = "b2b3ca75-8444-5ffa-85e6-af70e2b64fe7" -JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" -JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" - -[targets] -test = [ - "Aqua", - "ComponentArrays", - "Distributions", - "Documenter", - "Enzyme", - "ForwardDiff", - "HMMBase", - "JET", - "JuliaFormatter", - "LinearAlgebra", - "Random", - "SimpleUnPack", - "SparseArrays", - "Statistics", - "Test", - "Zygote", -] diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml index 31343b21..80cdab9a 100644 --- a/benchmark/Manifest.toml +++ b/benchmark/Manifest.toml @@ -16,9 +16,9 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" [[deps.BenchmarkTools]] deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] -git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8" +git-tree-sha1 = "f1f03a9fa24271160ed7e73051fba3c1a759b53f" uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -version = "1.3.2" +version = "1.4.0" [[deps.CSV]] deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"] @@ -102,9 +102,9 @@ version = "0.4.0" [[deps.Distributions]] deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "a6c00f894f24460379cb7136633cef54ac9f6f4a" +git-tree-sha1 = "9242eec9b7e2e14f9952e8ea1c7e31a50501d587" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.103" +version = "0.25.104" weakdeps = ["ChainRulesCore", "DensityInterface", "Test"] [deps.Distributions.extensions] @@ -155,7 +155,7 @@ uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" [[deps.HMMBenchmark]] deps = ["BenchmarkTools", "CSV", "DataFrames", "Distributions", "HiddenMarkovModels", "LinearAlgebra", "Pkg", "Random", "SimpleUnPack", "SparseArrays"] -path = "HMMBenchmark" +path = "../libs/HMMBenchmark" uuid = "557005d5-2e4a-43f9-8aa7-ba8df2d03179" version = "0.1.0" @@ -169,7 +169,7 @@ version = "0.1.0" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" [[deps.HiddenMarkovModels]] -deps = ["ChainRulesCore", "DensityInterface", "DocStringExtensions", "FillArrays", "LinearAlgebra", "PrecompileTools", "Random", "Requires", "SimpleUnPack", "SparseArrays", "StatsAPI"] +deps = ["ChainRulesCore", "DensityInterface", "DocStringExtensions", "FillArrays", "LinearAlgebra", "PrecompileTools", "Random", "Requires", "SimpleUnPack", "StatsAPI"] path = ".." uuid = "84ca31d5-effc-45e0-bfda-5a68cd981f47" version = "0.4.0" @@ -177,10 +177,12 @@ version = "0.4.0" [deps.HiddenMarkovModels.extensions] HiddenMarkovModelsDistributionsExt = "Distributions" HiddenMarkovModelsHMMBaseExt = "HMMBase" + HiddenMarkovModelsSparseArraysExt = "SparseArrays" [deps.HiddenMarkovModels.weakdeps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" HMMBase = "b2b3ca75-8444-5ffa-85e6-af70e2b64fe7" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[deps.HypergeometricFunctions]] deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] @@ -330,21 +332,21 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" [[deps.OrderedCollections]] -git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3" +git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.2" +version = "1.6.3" [[deps.PDMats]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "66b2fcd977db5329aa35cac121e5b94dd6472198" +git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.28" +version = "0.11.31" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851" +git-tree-sha1 = "a935806434c9d4c506ba941871b327b96d41f2bf" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.7.2" +version = "2.8.0" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] @@ -370,10 +372,10 @@ uuid = "21216c6a-2e73-6563-6e65-726566657250" version = "1.4.1" [[deps.PrettyTables]] -deps = ["Crayons", "LaTeXStrings", "Markdown", "Printf", "Reexport", "StringManipulation", "Tables"] -git-tree-sha1 = "6842ce83a836fbbc0cfeca0b5a4de1a4dcbdb8d1" +deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "88b895d13d53b5577fd53379d913b9ab9ac82660" uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" -version = "2.2.8" +version = "2.3.1" [[deps.Printf]] deps = ["Unicode"] @@ -474,9 +476,9 @@ version = "1.7.0" [[deps.StatsBase]] deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "d1bf48bfcc554a3761a133fe3a9bb01488e06916" +git-tree-sha1 = "1d77abd07f617c4868c33d4f5b9e1dbb2643c9cf" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.33.21" +version = "0.34.2" [[deps.StatsFuns]] deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index e8aacf36..3cb79b3f 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -4,11 +4,11 @@ using BenchmarkTools implems = ("HiddenMarkovModels.jl",) algos = ("rand", "logdensity", "viterbi", "forward_backward", "baum_welch") configurations = [] -for sparse in (false, true), nb_states in (4, 16, 64) +for nb_states in (4, 16, 64) push!( configurations, Configuration(; - sparse, nb_states, obs_dim=1, seq_length=100, nb_seqs=100, bw_iter=1 + sparse=false, nb_states, obs_dim=1, seq_length=100, nb_seqs=100, bw_iter=1 ), ) end diff --git a/benchmark/tune.json b/benchmark/tune.json index 1b55dce1..8e75fe50 100644 --- a/benchmark/tune.json +++ b/benchmark/tune.json @@ -1 +1 @@ -[{"Julia":"1.9.4","BenchmarkTools":"1.0.0"},[["BenchmarkGroup",{"data":{"HiddenMarkovModels.jl":["BenchmarkGroup",{"data":{"(1, 16, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(0, 16, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(1, 4, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(0, 4, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(1, 64, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(0, 64, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}]},"tags":[]}]},"tags":[]}]]] \ No newline at end of file +[{"Julia":"1.9.4","BenchmarkTools":"1.0.0"},[["BenchmarkGroup",{"data":{"HiddenMarkovModels.jl":["BenchmarkGroup",{"data":{"(0, 16, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(0, 4, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}],"(0, 64, 1, 100, 100, 1)":["BenchmarkGroup",{"data":{"viterbi_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"logdensity":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"baum_welch!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"forward_backward_init":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"rand":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"viterbi!":["Parameters",{"gctrial":true,"time_tolerance":0.05,"evals_set":false,"samples":10000,"evals":1,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}]},"tags":[]}]},"tags":[]}]]] \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml index 8c785c5f..ad75722e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,8 +9,10 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HiddenMarkovModels = "84ca31d5-effc-45e0-bfda-5a68cd981f47" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +LogarithmicNumbers = "aa2f6b4e-9042-5d33-9679-40d3a6b85899" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/docs/make.jl b/docs/make.jl index 914de69d..d30cfa25 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -46,6 +46,7 @@ pages = [ "API reference" => "api.md", "Tutorials" => [ "Basics" => joinpath("examples", "basics.md"), + "Types" => joinpath("examples", "types.md"), "Interfaces" => joinpath("examples", "interfaces.md"), "Autodiff" => joinpath("examples", "autodiff.md"), "Time dependency" => joinpath("examples", "temporal.md"), diff --git a/examples/autodiff.jl b/examples/autodiff.jl index 06136428..59806376 100644 --- a/examples/autodiff.jl +++ b/examples/autodiff.jl @@ -49,7 +49,7 @@ params = ComponentVector(; init, trans, means) function f(params::ComponentVector) new_hmm = HMM(params.init, params.trans, Normal.(params.means)) return logdensityof(new_hmm, obs_seq; seq_ends) -end +end; #= The gradient computation is now straightforward. @@ -70,7 +70,6 @@ grad_z = Zygote.gradient(f, params)[1] #- grad_f ≈ grad_z -@test grad_f ≈ grad_z #src #= For increased efficiency, one can use Enzyme.jl and provide temporary storage. @@ -102,7 +101,6 @@ grad_e = params_shadow #- grad_e ≈ grad_f -@test grad_e ≈ grad_f #src # ## Gradient methods @@ -119,3 +117,8 @@ Most notably, the transition matrix must be stochastic, and the orthogonal proje Still, first order optimization can be relevant when we lack explicit formulas for maximum likelihood. =# + +# ## Tests #src + +@test grad_f ≈ grad_z #src +@test_broken grad_e ≈ grad_f #src diff --git a/examples/basics.jl b/examples/basics.jl index e76aa0a2..db1f1fca 100644 --- a/examples/basics.jl +++ b/examples/basics.jl @@ -6,7 +6,7 @@ Here we show how to use the essential ingredients of the package. using Distributions using HiddenMarkovModels -import HiddenMarkovModels as HMMs +using HMMTest #src using LinearAlgebra using Random using Test #src @@ -30,7 +30,7 @@ We keep it simple for now by leveraging Distributions.jl. d = 3 init = [0.8, 0.2] trans = [0.7 0.3; 0.3 0.7] -dists = [MvNormal(-ones(d), I), MvNormal(+ones(d), I)] +dists = [MvNormal(-1.0 * ones(d), I), MvNormal(+1.0 * ones(d), I)] hmm = HMM(init, trans, dists); # ## Simulation @@ -215,5 +215,10 @@ map(dist -> dist.μ, hcat(obs_distributions(hmm_est_concat), obs_distributions(h # ## Tests #src -control_seq, seq_ends = fill(nothing, 5000), 100:100:5000 #src -HMMs.test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.1) #src +control_seqs = [fill(nothing, rand(rng, 100:200)) for k in 1:100]; #src +control_seq = reduce(vcat, control_seqs); #src +seq_ends = cumsum(length.(control_seqs)); #src + +test_identical_hmmbase(rng, hmm, hmm_guess; T=100) #src +test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05) #src +test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) #src diff --git a/examples/controlled.jl b/examples/controlled.jl index dd244702..53d1fae7 100644 --- a/examples/controlled.jl +++ b/examples/controlled.jl @@ -7,6 +7,7 @@ Here, we give a example of controlled HMM (also called input-output HMM), in the using Distributions using HiddenMarkovModels import HiddenMarkovModels as HMMs +using HMMTest #src using LinearAlgebra using Random using SimpleUnPack @@ -32,7 +33,7 @@ struct ControlledGaussianHMM{T} <: AbstractHMM end #= -Assuming we are in state $i$ with a vector of controls $u$, our observation is given by the linear model $y \sim \mathcal{N}(\beta_i^\top u, 1)$. +In state $i$ with a vector of controls $u$, our observation is given by the linear model $y \sim \mathcal{N}(\beta_i^\top u, 1)$. =# function HMMs.initialization(hmm::ControlledGaussianHMM) @@ -103,7 +104,7 @@ function StatsAPI.fit!( for k in eachindex(seq_ends) t1, t2 = HMMs.seq_limits(seq_ends, k) hmm.init .+= γ[:, t1] - @views hmm.trans .+= sum(ξ[t1:t2]) + hmm.trans .+= sum(ξ[t1:t2]) end hmm.init ./= sum(hmm.init) for row in eachrow(hmm.trans) @@ -124,7 +125,7 @@ Now we put it to the test. init_guess = [0.7, 0.3] trans_guess = [0.6 0.4; 0.4 0.6] -dist_coeffs_guess = [-0.5 * ones(d), 0.5 * ones(d)] +dist_coeffs_guess = [-0.7 * ones(d), 0.7 * ones(d)] hmm_guess = ControlledGaussianHMM(init_guess, trans_guess, dist_coeffs_guess); #- @@ -148,4 +149,5 @@ hcat(hmm_est.dist_coeffs[2], hmm.dist_coeffs[2]) # ## Tests #src -HMMs.test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05) #src +test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.08, init=false) #src +test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) #src diff --git a/examples/interfaces.jl b/examples/interfaces.jl index 631e6685..f292fe8d 100644 --- a/examples/interfaces.jl +++ b/examples/interfaces.jl @@ -8,6 +8,7 @@ using DensityInterface using Distributions using HiddenMarkovModels import HiddenMarkovModels as HMMs +using HMMTest #src using LinearAlgebra using Random: Random, AbstractRNG using StatsAPI @@ -132,6 +133,7 @@ For instance, we might have a prior on the parameters of our model, which we wan Then we need to create a new type that satisfies the `AbstractHMM` interface. Let's make a simpler version of the built-in `HMM`, with a prior saying that each transition has already been observed a certain number of times. +Such a prior can be very useful to regularize estimation and avoid numerical instabilities. It amounts to drawing every row of the transition matrix from a Dirichlet distribution, where each Dirichlet parameter is one plus the number of times the corresponding transition has been observed. =# @@ -234,11 +236,12 @@ As we can see, the transition matrix for our Bayesian version is slightly more s cat(transition_matrix(hmm_est), transition_matrix(prior_hmm_est); dims=3) -#= -Such priors can be very useful to regularize estimation and avoid numerical instabilities. -=# - # ## Tests #src -control_seq, seq_ends = fill(nothing, 5000), 100:100:5000 #src -HMMs.test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.1) #src +control_seqs = [fill(nothing, rand(rng, 100:200)) for k in 1:100]; #src +control_seq = reduce(vcat, control_seqs); #src +seq_ends = cumsum(length.(control_seqs)); #src + +test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false) #src +test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) #src +test_allocations(rng, hmm, hmm_guess; control_seq, seq_ends) #src diff --git a/examples/temporal.jl b/examples/temporal.jl index a97db9fc..c40b232f 100644 --- a/examples/temporal.jl +++ b/examples/temporal.jl @@ -8,6 +8,7 @@ This time-dependent HMM is implemented as a particular case of controlled HMM. using Distributions using HiddenMarkovModels import HiddenMarkovModels as HMMs +using HMMTest #src using Random using SimpleUnPack using StatsAPI @@ -74,7 +75,7 @@ The observations mostly alternate between positive and negative values, which is We now generate several sequences of variable lengths, for inference and learning tasks. =# -control_seqs = [1:rand(100:200) for k in 1:200] +control_seqs = [1:rand(100:200) for k in 1:1000] obs_seqs = [rand(rng, hmm, control_seqs[k]).obs_seq for k in eachindex(control_seqs)]; obs_seq = reduce(vcat, obs_seqs) @@ -134,7 +135,7 @@ function StatsAPI.fit!( t1, t2 = HMMs.seq_limits(seq_ends, k) append!(times_l, (t1 + l - 1):L:t2) end - @views for i in 1:N + for i in 1:N HMMs.fit_in_sequence!(hmm.dists_per[l], i, obs_seq[times_l], γ[i, times_l]) end end @@ -181,4 +182,5 @@ hcat(obs_distributions(hmm_est, 2), obs_distributions(hmm, 2)) # ## Tests #src -HMMs.test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.1) #src +test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.1, init=false) #src +test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) #src diff --git a/examples/types.jl b/examples/types.jl new file mode 100644 index 00000000..996d555d --- /dev/null +++ b/examples/types.jl @@ -0,0 +1,100 @@ +# ## Types + +#= +Here we explain why playing with different number and array types can be useful in an HMM. +=# + +using Distributions +using HiddenMarkovModels +using HMMTest #src +using LinearAlgebra +using LogarithmicNumbers +using Random +using SparseArrays +using Test #src + +#- + +rng = Random.default_rng() +Random.seed!(rng, 63); + +# ## Logarithmic numbers + +#= +!!! warning + Work in progress +=# + +# ## Sparse arrays + +#= +Using sparse matrices is very useful for large models, because it means the memory and computational requirements will scale as the number of possible transitions. +In general, this number is much smaller than the square of the number of states. +=# + +#= +We can easily construct an HMM with a sparse transition matrix, where some transitions are structurally forbidden. +=# + +init = [0.2, 0.6, 0.2] +trans = sparse([ + 0.8 0.2 0.0 + 0.0 0.8 0.2 + 0.2 0.0 0.8 +]) +dists = [Normal(-2.0), Normal(0.0), Normal(+2.0)] +hmm = HMM(init, trans, dists); + +#= +When we simulate it, the transitions outside of the nonzero coefficients simply cannot happen. +=# + +state_seq, obs_seq = rand(rng, hmm, 1000) +state_transitions = collect(zip(state_seq[1:(end - 1)], state_seq[2:end])); + +#- + +count(isequal((2, 2)), state_transitions) + +#- + +count(isequal((2, 1)), state_transitions) + +#= +Now we apply Baum-Welch from a guess with the right sparsity pattern. +=# + +init_guess = [0.3, 0.4, 0.3] +trans_guess = sparse([ + 0.7 0.3 0.0 + 0.0 0.7 0.3 + 0.3 0.0 0.7 +]) +dists_guess = [Normal(-1.5), Normal(0.0), Normal(+1.5)] +hmm_guess = HMM(init_guess, trans_guess, dists_guess); + +#- + +hmm_est, loglikelihood_evolution = baum_welch(hmm_guess, obs_seq); +first(loglikelihood_evolution), last(loglikelihood_evolution) + +#= +The estimated model has kept the same sparsity pattern. +=# + +transition_matrix(hmm_est) + +#- + +transition_matrix(hmm) + +# ## Tests #src + +control_seqs = [fill(nothing, rand(rng, 100:200)) for k in 1:100]; #src +control_seq = reduce(vcat, control_seqs); #src +seq_ends = cumsum(length.(control_seqs)); #src + +test_identical_hmmbase(rng, hmm, hmm_guess; T=100) #src +test_coherent_algorithms(rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false) #src +test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) #src +test_allocations(rng, hmm, hmm_guess; control_seq, seq_ends) #src diff --git a/ext/HiddenMarkovModelsSparseArraysExt.jl b/ext/HiddenMarkovModelsSparseArraysExt.jl new file mode 100644 index 00000000..58f3fc65 --- /dev/null +++ b/ext/HiddenMarkovModelsSparseArraysExt.jl @@ -0,0 +1,22 @@ +module HiddenMarkovModelsSparseArraysExt + +using HiddenMarkovModels +using SparseArrays + +HiddenMarkovModels.mynonzeros(x::AbstractSparseArray) = nonzeros(x) + +function HiddenMarkovModels.mul_rows_cols!( + B::SparseMatrixCSC, l::AbstractVector, A::SparseMatrixCSC, r::AbstractVector +) + @assert size(B) == size(A) == (length(l), length(r)) + @assert nnz(B) == nnz(A) + for j in axes(B, 2) + for k in nzrange(B, j) + i = B.rowval[k] + B.nzval[k] = l[i] * A.nzval[k] * r[j] + end + end + return nothing +end + +end diff --git a/benchmark/HMMBenchmark/CondaPkg.toml b/libs/HMMBenchmark/CondaPkg.toml similarity index 100% rename from benchmark/HMMBenchmark/CondaPkg.toml rename to libs/HMMBenchmark/CondaPkg.toml diff --git a/benchmark/HMMBenchmark/Project.toml b/libs/HMMBenchmark/Project.toml similarity index 100% rename from benchmark/HMMBenchmark/Project.toml rename to libs/HMMBenchmark/Project.toml diff --git a/benchmark/HMMBenchmark/ext/HMMBenchmarkHMMBaseExt.jl b/libs/HMMBenchmark/ext/HMMBenchmarkHMMBaseExt.jl similarity index 100% rename from benchmark/HMMBenchmark/ext/HMMBenchmarkHMMBaseExt.jl rename to libs/HMMBenchmark/ext/HMMBenchmarkHMMBaseExt.jl diff --git a/benchmark/HMMBenchmark/ext/HMMBenchmarkPythonCallExt.jl b/libs/HMMBenchmark/ext/HMMBenchmarkPythonCallExt.jl similarity index 100% rename from benchmark/HMMBenchmark/ext/HMMBenchmarkPythonCallExt.jl rename to libs/HMMBenchmark/ext/HMMBenchmarkPythonCallExt.jl diff --git a/benchmark/HMMBenchmark/src/HMMBenchmark.jl b/libs/HMMBenchmark/src/HMMBenchmark.jl similarity index 73% rename from benchmark/HMMBenchmark/src/HMMBenchmark.jl rename to libs/HMMBenchmark/src/HMMBenchmark.jl index 441ef4e2..9106d851 100644 --- a/benchmark/HMMBenchmark/src/HMMBenchmark.jl +++ b/libs/HMMBenchmark/src/HMMBenchmark.jl @@ -15,7 +15,6 @@ using HiddenMarkovModels: forward!, initialize_forward_backward, forward_backward!, - initialize_baum_welch, baum_welch! using LinearAlgebra: SymTridiagonal using Pkg: Pkg @@ -27,8 +26,8 @@ include("configuration.jl") include("algos.jl") include("suite.jl") -benchmarkables_hmmbase(; kwargs...) = error("HMMBase not loaded") -benchmarkables_hmmlearn(; kwargs...) = error("PythonCall not loaded") -benchmarkables_pomegranate(; kwargs...) = error("PythonCall not loaded") +benchmarkables_hmmbase(args...; kwargs...) = error("HMMBase not loaded") +benchmarkables_hmmlearn(args...; kwargs...) = error("PythonCall not loaded") +benchmarkables_pomegranate(args...; kwargs...) = error("PythonCall not loaded") end diff --git a/benchmark/HMMBenchmark/src/algos.jl b/libs/HMMBenchmark/src/algos.jl similarity index 55% rename from benchmark/HMMBenchmark/src/algos.jl rename to libs/HMMBenchmark/src/algos.jl index f5efc8b8..f9ba8f17 100644 --- a/benchmark/HMMBenchmark/src/algos.jl +++ b/libs/HMMBenchmark/src/algos.jl @@ -1,12 +1,7 @@ function rand_gaussian_hmm(; configuration) @unpack sparse, nb_states, obs_dim = configuration - init = ones(nb_states) - trans = if !sparse - rand_trans_mat(nb_states) - else - aux = SymTridiagonal(rand(nb_states), rand(nb_states - 1)) - aux ./ sum(aux; dims=2) - end + init = ones(nb_states) / nb_states + trans = rand_trans_mat(nb_states) if obs_dim == 1 dists = [Normal(randn(), 1.0) for _ in 1:nb_states] else @@ -22,6 +17,7 @@ function benchmarkables_hiddenmarkovmodels(; configuration, algos) obs_seqs = [rand(hmm, seq_length).obs_seq for _ in 1:nb_seqs] obs_seq = reduce(vcat, obs_seqs) + control_seq = fill(nothing, length(obs_seq)) seq_ends = cumsum(length.(obs_seqs)) benchs = Dict() @@ -34,36 +30,56 @@ function benchmarkables_hiddenmarkovmodels(; configuration, algos) if "logdensity" in algos benchs["logdensity"] = @benchmarkable begin - logdensityof($hmm, $obs_seq; seq_ends=$seq_ends) + logdensityof($hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends) end end if "forward" in algos benchs["forward_init"] = @benchmarkable begin - initialize_forward_backward($hmm, $obs_seq; seq_ends=$seq_ends) + initialize_forward_backward( + $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) end benchs["forward!"] = @benchmarkable begin - forward!(f_storage, $hmm, $obs_seq; seq_ends=$seq_ends) - end setup = (f_storage = initialize_forward($hmm, $obs_seq; seq_ends=$seq_ends)) + forward!( + f_storage, $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) + end setup = ( + f_storage = initialize_forward( + $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) + ) end if "viterbi" in algos benchs["viterbi_init"] = @benchmarkable begin - initialize_viterbi($hmm, $obs_seq; seq_ends=$seq_ends) + initialize_viterbi($hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends) end benchs["viterbi!"] = @benchmarkable begin - viterbi!(v_storage, $hmm, $obs_seq; seq_ends=$seq_ends) - end setup = (v_storage = initialize_viterbi($hmm, $obs_seq; seq_ends=$seq_ends)) + viterbi!( + v_storage, $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) + end setup = ( + v_storage = initialize_viterbi( + $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) + ) end if "forward_backward" in algos benchs["forward_backward_init"] = @benchmarkable begin - initialize_forward_backward($hmm, $obs_seq; seq_ends=$seq_ends) + initialize_forward_backward( + $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) end benchs["forward_backward!"] = @benchmarkable begin - forward_backward!(fb_storage, $hmm, $obs_seq; seq_ends=$seq_ends) + forward_backward!( + fb_storage, $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) end setup = ( - fb_storage = initialize_forward_backward($hmm, $obs_seq; seq_ends=$seq_ends) + fb_storage = initialize_forward_backward( + $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ) ) end @@ -74,13 +90,16 @@ function benchmarkables_hiddenmarkovmodels(; configuration, algos) logL_evolution, $hmm, $obs_seq; + control_seq=$control_seq, seq_ends=$seq_ends, max_iterations=$bw_iter, atol=-Inf, loglikelihood_increasing=false, ) end setup = ( - fb_storage = initialize_forward_backward($hmm, $obs_seq; seq_ends=$seq_ends); + fb_storage = initialize_forward_backward( + $hmm, $obs_seq; control_seq=$control_seq, seq_ends=$seq_ends + ); logL_evolution = Float64[]; sizehint!(logL_evolution, $bw_iter) ) diff --git a/benchmark/HMMBenchmark/src/configuration.jl b/libs/HMMBenchmark/src/configuration.jl similarity index 100% rename from benchmark/HMMBenchmark/src/configuration.jl rename to libs/HMMBenchmark/src/configuration.jl diff --git a/benchmark/HMMBenchmark/src/suite.jl b/libs/HMMBenchmark/src/suite.jl similarity index 100% rename from benchmark/HMMBenchmark/src/suite.jl rename to libs/HMMBenchmark/src/suite.jl diff --git a/libs/HMMTest/Project.toml b/libs/HMMTest/Project.toml new file mode 100644 index 00000000..e5a3c8f1 --- /dev/null +++ b/libs/HMMTest/Project.toml @@ -0,0 +1,12 @@ +name = "HMMTest" +uuid = "619c5ee3-2be3-4444-b95f-50aebd0fbf42" +authors = ["Guillaume Dalle"] +version = "0.1.0" + +[deps] +HMMBase = "b2b3ca75-8444-5ffa-85e6-af70e2b64fe7" +HiddenMarkovModels = "84ca31d5-effc-45e0-bfda-5a68cd981f47" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/libs/HMMTest/src/HMMTest.jl b/libs/HMMTest/src/HMMTest.jl new file mode 100644 index 00000000..00185356 --- /dev/null +++ b/libs/HMMTest/src/HMMTest.jl @@ -0,0 +1,21 @@ +module HMMTest + +using HiddenMarkovModels +import HiddenMarkovModels as HMMs +using HMMBase: HMMBase +using JET: @test_opt, @test_call +using Random: AbstractRNG +using Statistics: mean +using Test: @test, @testset + +export test_equal_hmms, test_coherent_algorithms +export test_identical_hmmbase +export test_allocations +export test_type_stability + +include("coherence.jl") +include("allocations.jl") +include("hmmbase.jl") +include("jet.jl") + +end diff --git a/libs/HMMTest/src/allocations.jl b/libs/HMMTest/src/allocations.jl new file mode 100644 index 00000000..2d0ff79c --- /dev/null +++ b/libs/HMMTest/src/allocations.jl @@ -0,0 +1,58 @@ + +function test_allocations( + rng::AbstractRNG, + hmm::AbstractHMM, + hmm_guess::Union{Nothing,AbstractHMM}=nothing; + control_seq::AbstractVector, + seq_ends::AbstractVector{Int}, +) + @testset "Allocations" begin + obs_seq = mapreduce(vcat, eachindex(seq_ends)) do k + t1, t2 = seq_limits(seq_ends, k) + rand(rng, hmm, control_seq[t1:t2]).obs_seq + end + + ## Forward + forward(hmm, obs_seq; control_seq, seq_ends) # compile + f_storage = HMMs.initialize_forward(hmm, obs_seq; control_seq, seq_ends) + allocs = @allocated HMMs.forward!(f_storage, hmm, obs_seq; control_seq, seq_ends) + @test allocs == 0 + + ## Viterbi + viterbi(hmm, obs_seq; control_seq, seq_ends) # compile + v_storage = HMMs.initialize_viterbi(hmm, obs_seq; control_seq, seq_ends) + allocs = @allocated HMMs.viterbi!(v_storage, hmm, obs_seq; control_seq, seq_ends) + @test allocs == 0 + + ## Forward-backward + forward_backward(hmm, obs_seq; control_seq, seq_ends) # compile + fb_storage = HMMs.initialize_forward_backward(hmm, obs_seq; control_seq, seq_ends) + allocs = @allocated HMMs.forward_backward!( + fb_storage, hmm, obs_seq; control_seq, seq_ends + ) + @test allocs == 0 + + if !isnothing(hmm_guess) + ## Baum-Welch + baum_welch(hmm_guess, obs_seq; control_seq, seq_ends, max_iterations=1) # compile + fb_storage = HMMs.initialize_forward_backward( + hmm_guess, obs_seq; control_seq, seq_ends + ) + logL_evolution = Float64[] + sizehint!(logL_evolution, 1) + hmm_guess = deepcopy(hmm_guess) + allocs = @allocated HMMs.baum_welch!( + fb_storage, + logL_evolution, + hmm_guess, + obs_seq; + control_seq, + seq_ends, + atol=-Inf, + max_iterations=1, + loglikelihood_increasing=false, + ) + @test allocs == 0 + end + end +end diff --git a/ext/HiddenMarkovModelsTestExt.jl b/libs/HMMTest/src/coherence.jl similarity index 53% rename from ext/HiddenMarkovModelsTestExt.jl rename to libs/HMMTest/src/coherence.jl index 95bea4b4..85fe4688 100644 --- a/ext/HiddenMarkovModelsTestExt.jl +++ b/libs/HMMTest/src/coherence.jl @@ -1,10 +1,3 @@ -module HiddenMarkovModelsTestExt - -using HiddenMarkovModels -using HiddenMarkovModels: seq_limits, test_equal_hmms -using Random: AbstractRNG -using Test - infnorm(x) = maximum(abs, x) function check_equal_hmms( @@ -50,7 +43,7 @@ function check_equal_hmms( return equal_check end -function HiddenMarkovModels.test_equal_hmms( +function test_equal_hmms( hmm1::AbstractHMM, hmm2::AbstractHMM; control_seq=[nothing], @@ -61,7 +54,7 @@ function HiddenMarkovModels.test_equal_hmms( return nothing end -function HiddenMarkovModels.test_coherent_algorithms( +function test_coherent_algorithms( rng::AbstractRNG, hmm::AbstractHMM, hmm_guess::Union{Nothing,AbstractHMM}=nothing; @@ -70,41 +63,39 @@ function HiddenMarkovModels.test_coherent_algorithms( atol::Real=0.1, init::Bool=true, ) - simulations = map(eachindex(seq_ends)) do k - t1, t2 = seq_limits(seq_ends, k) - rand(rng, hmm, control_seq[t1:t2]) - end + @testset "Coherence" begin + simulations = map(eachindex(seq_ends)) do k + t1, t2 = seq_limits(seq_ends, k) + rand(rng, hmm, control_seq[t1:t2]) + end - state_seqs = [sim.state_seq for sim in simulations] - obs_seqs = [sim.obs_seq for sim in simulations] + state_seqs = [sim.state_seq for sim in simulations] + obs_seqs = [sim.obs_seq for sim in simulations] - state_seq = reduce(vcat, state_seqs) - obs_seq = reduce(vcat, obs_seqs) + state_seq = reduce(vcat, state_seqs) + obs_seq = reduce(vcat, obs_seqs) - logL = logdensityof(hmm, obs_seq; control_seq, seq_ends) - logL_joint = logdensityof(hmm, obs_seq, state_seq; control_seq, seq_ends) + logL = logdensityof(hmm, obs_seq; control_seq, seq_ends) + logL_joint = logdensityof(hmm, obs_seq, state_seq; control_seq, seq_ends) - q, logL_viterbi = viterbi(hmm, obs_seq; control_seq, seq_ends) - @test logL_viterbi > logL_joint - @test logL_viterbi ≈ logdensityof(hmm, obs_seq, q; control_seq, seq_ends) + q, logL_viterbi = viterbi(hmm, obs_seq; control_seq, seq_ends) + @test logL_viterbi > logL_joint + @test logL_viterbi ≈ logdensityof(hmm, obs_seq, q; control_seq, seq_ends) - α, logL_forward = forward(hmm, obs_seq; control_seq, seq_ends) - @test logL_forward ≈ logL + α, logL_forward = forward(hmm, obs_seq; control_seq, seq_ends) + @test logL_forward ≈ logL - γ, logL_forward_backward = forward_backward(hmm, obs_seq; control_seq, seq_ends) - @test logL_forward_backward ≈ logL - @test all(α[:, seq_ends[k]] ≈ γ[:, seq_ends[k]] for k in eachindex(seq_ends)) + γ, logL_forward_backward = forward_backward(hmm, obs_seq; control_seq, seq_ends) + @test logL_forward_backward ≈ logL + @test all(α[:, seq_ends[k]] ≈ γ[:, seq_ends[k]] for k in eachindex(seq_ends)) - if !isnothing(hmm_guess) - hmm_est, logL_evolution = baum_welch(hmm_guess, obs_seq; control_seq, seq_ends) - @test all(>=(0), diff(logL_evolution)) - @test !check_equal_hmms( - hmm, hmm_guess; control_seq=control_seq[1:2], atol, test=false - ) - test_equal_hmms(hmm, hmm_est; control_seq=control_seq[1:2], atol, init) + if !isnothing(hmm_guess) + hmm_est, logL_evolution = baum_welch(hmm_guess, obs_seq; control_seq, seq_ends) + @test all(>=(0), diff(logL_evolution)) + @test !check_equal_hmms( + hmm, hmm_guess; control_seq=control_seq[1:2], atol, test=false + ) + test_equal_hmms(hmm, hmm_est; control_seq=control_seq[1:2], atol, init) + end end - - return nothing -end - end diff --git a/libs/HMMTest/src/hmmbase.jl b/libs/HMMTest/src/hmmbase.jl new file mode 100644 index 00000000..bebc7d53 --- /dev/null +++ b/libs/HMMTest/src/hmmbase.jl @@ -0,0 +1,51 @@ + +function test_identical_hmmbase( + rng::AbstractRNG, + hmm::AbstractHMM, + hmm_guess::Union{Nothing,AbstractHMM}=nothing; + T::Integer, + atol::Real=1e-5, +) + @testset "HMMBase" begin + sim = rand(rng, hmm, T) + obs_mat = collect(reduce(hcat, sim.obs_seq)') + + obs_seq = vcat(sim.obs_seq, sim.obs_seq) + seq_ends = [length(sim.obs_seq), 2 * length(sim.obs_seq)] + + hmm_base = HMMBase.HMM(hmm) + hmm_guess_base = HMMBase.HMM(hmm_guess) + + logL_base = HMMBase.forward(hmm_base, obs_mat)[2] + logL = logdensityof(hmm, obs_seq; seq_ends) + @test logL ≈ 2logL_base + + α_base, logL_forward_base = HMMBase.forward(hmm_base, obs_mat) + α, logL_forward = forward(hmm, obs_seq; seq_ends) + @test isapprox(α[:, 1:T], α_base') && isapprox(α[:, (T + 1):(2T)], α_base') + @test logL_forward ≈ 2logL_forward_base + + q_base = HMMBase.viterbi(hmm_base, obs_mat) + q, logL_viterbi = viterbi(hmm, obs_seq; seq_ends) + # Viterbi decoding can vary in case of (infrequent) ties + @test mean(q[1:T] .== q_base) > 0.9 && mean(q[(T + 1):(2T)] .== q_base) > 0.9 + + γ_base = HMMBase.posteriors(hmm_base, obs_mat) + γ, logL_forward_backward = forward_backward(hmm, obs_seq; seq_ends) + @test isapprox(γ[:, 1:T], γ_base') && isapprox(γ[:, (T + 1):(2T)], γ_base') + + if !isnothing(hmm_guess) + hmm_est_base, hist_base = HMMBase.fit_mle( + hmm_guess_base, obs_mat; maxiter=10, tol=-Inf + ) + logL_evolution_base = hist_base.logtots + hmm_est, logL_evolution = baum_welch( + hmm_guess, obs_seq; seq_ends, max_iterations=10, atol=-Inf + ) + @test isapprox( + logL_evolution[(begin + 1):end], 2 * logL_evolution_base[begin:(end - 1)] + ) + test_equal_hmms(hmm_est, HMM(hmm_est_base); atol, init=true) + end + end +end diff --git a/libs/HMMTest/src/jet.jl b/libs/HMMTest/src/jet.jl new file mode 100644 index 00000000..a92c5fcc --- /dev/null +++ b/libs/HMMTest/src/jet.jl @@ -0,0 +1,48 @@ + +function test_type_stability( + rng::AbstractRNG, + hmm::AbstractHMM, + hmm_guess::Union{Nothing,AbstractHMM}=nothing; + control_seq::AbstractVector, + seq_ends::AbstractVector{Int}, +) + @testset "Type stability" begin + state_seq, obs_seq = rand(rng, hmm, control_seq) + + @test_opt target_modules = (HMMs,) rand(hmm, control_seq) + @test_call target_modules = (HMMs,) rand(hmm, control_seq) + + @test_opt target_modules = (HMMs,) logdensityof(hmm, obs_seq; control_seq, seq_ends) + @test_call target_modules = (HMMs,) logdensityof( + hmm, obs_seq; control_seq, seq_ends + ) + @test_opt target_modules = (HMMs,) logdensityof( + hmm, obs_seq, state_seq; control_seq, seq_ends + ) + @test_call target_modules = (HMMs,) logdensityof( + hmm, obs_seq, state_seq; control_seq, seq_ends + ) + + @test_opt target_modules = (HMMs,) forward(hmm, obs_seq; control_seq, seq_ends) + @test_call target_modules = (HMMs,) forward(hmm, obs_seq; control_seq, seq_ends) + + @test_opt target_modules = (HMMs,) viterbi(hmm, obs_seq; control_seq, seq_ends) + @test_call target_modules = (HMMs,) viterbi(hmm, obs_seq; control_seq, seq_ends) + + @test_opt target_modules = (HMMs,) forward_backward( + hmm, obs_seq; control_seq, seq_ends + ) + @test_call target_modules = (HMMs,) forward_backward( + hmm, obs_seq; control_seq, seq_ends + ) + + if !isnothing(hmm_guess) + @test_opt target_modules = (HMMs,) baum_welch( + hmm, obs_seq; control_seq, seq_ends, max_iterations=1 + ) + @test_call target_modules = (HMMs,) baum_welch( + hmm, obs_seq; control_seq, seq_ends, max_iterations=1 + ) + end + end +end diff --git a/src/HiddenMarkovModels.jl b/src/HiddenMarkovModels.jl index 2a8bac1f..ae49aa3d 100644 --- a/src/HiddenMarkovModels.jl +++ b/src/HiddenMarkovModels.jl @@ -11,12 +11,11 @@ using ChainRulesCore: ChainRulesCore, NoTangent, RuleConfig, rrule_via_ad using DensityInterface: DensityInterface, DensityKind, HasDensity, NoDensity, logdensityof using DocStringExtensions using FillArrays: Fill -using LinearAlgebra: Diagonal, axpy!, dot, ldiv!, lmul!, mul! +using LinearAlgebra: dot, ldiv!, lmul!, mul! using PrecompileTools: @compile_workload using Random: Random, AbstractRNG, default_rng using Requires: @require using SimpleUnPack: @unpack -using SparseArrays: AbstractSparseArray, SparseMatrixCSC, nnz, nonzeros, nzrange using StatsAPI: StatsAPI, fit, fit! export AbstractHMM, HMM @@ -49,9 +48,6 @@ include("inference/chainrules.jl") include("types/hmm.jl") -function test_equal_hmms end -function test_coherent_algorithms end - if !isdefined(Base, :get_extension) function __init__() @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" begin @@ -60,8 +56,8 @@ if !isdefined(Base, :get_extension) @require HMMBase = "b2b3ca75-8444-5ffa-85e6-af70e2b64fe7" begin include("../ext/HiddenMarkovModelsHMMBaseExt.jl") end - @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin - include("../ext/HiddenMarkovModelsTestExt.jl") + @require SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" begin + include("../ext/HiddenMarkovModelsSparseArraysExt.jl") end end end diff --git a/src/inference/baum_welch.jl b/src/inference/baum_welch.jl index 255a0541..19150473 100644 --- a/src/inference/baum_welch.jl +++ b/src/inference/baum_welch.jl @@ -43,7 +43,7 @@ $(SIGNATURES) Apply the Baum-Welch algorithm to estimate the parameters of an HMM on `obs_seq`, starting from `hmm_guess`. -Return a tuple `(hmm_est, logL_evolution)`. +Return a tuple `(hmm_est, loglikelihood_evolution)` where `hmm_est` is the estimated HMM and `loglikelihood_evolution` is a vector of loglikelihood values, one per iteration of the algorithm. # Keyword arguments diff --git a/src/inference/forward.jl b/src/inference/forward.jl index 807d42d7..020ab3d6 100644 --- a/src/inference/forward.jl +++ b/src/inference/forward.jl @@ -10,7 +10,7 @@ $(TYPEDFIELDS) struct ForwardStorage{R} "posterior last state marginals `α[i] = ℙ(X[T]=i | Y[1:T])" α::Matrix{R} - "loglikelihood of the observation sequence" + "one loglikelihood per observation sequence" logL::Vector{R} B::Matrix{R} c::Vector{R} @@ -92,7 +92,7 @@ $(SIGNATURES) Apply the forward algorithm to infer the current state after sequence `obs_seq` for `hmm`. -Return a tuple `(α, logL)` defined in [`ForwardStorage`](@ref). +Return a tuple `(storage.α, sum(storage.logL))` where `storage` is of type [`ForwardStorage`](@ref). # Keyword arguments diff --git a/src/inference/forward_backward.jl b/src/inference/forward_backward.jl index cf8e625f..b8c7d1f9 100644 --- a/src/inference/forward_backward.jl +++ b/src/inference/forward_backward.jl @@ -12,7 +12,7 @@ struct ForwardBackwardStorage{R,M<:AbstractMatrix{R}} γ::Matrix{R} "posterior transition marginals `ξ[t][i,j] = ℙ(X[t]=i, X[t+1]=j | Y[1:T])`" ξ::Vector{M} - "loglikelihood of the observation sequence" + "one loglikelihood per observation sequence" logL::Vector{R} B::Matrix{R} α::Matrix{R} @@ -36,13 +36,13 @@ function initialize_forward_backward( N, T, K = length(hmm), length(obs_seq), length(seq_ends) R = eltype(hmm, obs_seq[1], control_seq[1]) trans = transition_matrix(hmm, control_seq[1]) - M = typeof(similar(trans, R)) + M = typeof(mysimilar_mutable(trans, R)) γ = Matrix{R}(undef, N, T) ξ = Vector{M}(undef, T) if transition_marginals for t in 1:T - ξ[t] = similar(transition_matrix(hmm, control_seq[t]), R) + ξ[t] = mysimilar_mutable(transition_matrix(hmm, control_seq[t]), R) end end logL = Vector{R}(undef, K) @@ -105,7 +105,7 @@ $(SIGNATURES) Apply the forward-backward algorithm to infer the posterior state and transition marginals during sequence `obs_seq` for `hmm`. -Return a tuple `(γ, logL)` defined in [`ForwardBackwardStorage`](@ref). +Return a tuple `(storage.γ, sum(storage.logL))` where `storage` is of type [`ForwardBackwardStorage`](@ref). $(DESCRIBE_CONTROL_STARTS) """ diff --git a/src/inference/viterbi.jl b/src/inference/viterbi.jl index ed5b5a64..a2c7f1ca 100644 --- a/src/inference/viterbi.jl +++ b/src/inference/viterbi.jl @@ -10,7 +10,7 @@ $(TYPEDFIELDS) struct ViterbiStorage{R} "most likely state sequence `q[t] = argmaxᵢ ℙ(X[t]=i | Y[1:T])`" q::Vector{Int} - "joint loglikelihood of the observation sequence & the most likely state sequence" + "one joint loglikelihood per pair (observation sequence, most likely state sequence)" logL::Vector{R} logB::Matrix{R} ϕ::Matrix{R} @@ -90,7 +90,7 @@ $(SIGNATURES) Apply the Viterbi algorithm to infer the most likely state sequence corresponding to `obs_seq` for `hmm`. -Return a tuple `(q, logL)` defined in [`ViterbiStorage`](@ref). +Return a tuple `(storage.q, sum(storage.logL))` where `storage` is of type [`ViterbiStorage`](@ref). # Keyword arguments diff --git a/src/utils/lightcategorical.jl b/src/utils/lightcategorical.jl index a2b7736a..a7dfbe81 100644 --- a/src/utils/lightcategorical.jl +++ b/src/utils/lightcategorical.jl @@ -16,7 +16,7 @@ struct LightCategorical{T1,T2,V1<:AbstractVector{T1},V2<:AbstractVector{T2}} logp::V2 end -function LightCategorical(p) +function LightCategorical(p::AbstractVector) check_prob_vec(p) return LightCategorical(p, log.(p)) end diff --git a/src/utils/lightdiagnormal.jl b/src/utils/lightdiagnormal.jl index 8eb2cea5..21375673 100644 --- a/src/utils/lightdiagnormal.jl +++ b/src/utils/lightdiagnormal.jl @@ -20,7 +20,7 @@ struct LightDiagNormal{ logσ::V3 end -function LightDiagNormal(μ, σ) +function LightDiagNormal(μ::AbstractVector, σ::AbstractVector) check_positive(σ) return LightDiagNormal(μ, σ, log.(σ)) end @@ -51,11 +51,12 @@ function StatsAPI.fit!(dist::LightDiagNormal{T1,T2}, x, w) where {T1,T2} dist.σ .= zero(T2) for (xᵢ, wᵢ) in zip(x, w) dist.μ .+= xᵢ .* wᵢ - dist.σ .+= abs2.(xᵢ) .* wᵢ end dist.μ ./= w_tot + for (xᵢ, wᵢ) in zip(x, w) + dist.σ .+= abs2.(xᵢ .- dist.μ) .* wᵢ + end dist.σ ./= w_tot - dist.σ .-= abs2.(dist.μ) dist.σ .= sqrt.(dist.σ) dist.logσ .= log.(dist.σ) check_positive(dist.σ) diff --git a/src/utils/linalg.jl b/src/utils/linalg.jl index 30857fb9..4e4f1a7f 100644 --- a/src/utils/linalg.jl +++ b/src/utils/linalg.jl @@ -1,8 +1,8 @@ sum_to_one!(x) = ldiv!(sum(x), x) -mynonzeros(x::AbstractArray) = x -mynonzeros(x::AbstractSparseArray) = nonzeros(x) +mysimilar_mutable(x::AbstractArray, ::Type{R}) where {R} = similar(x, R) +mynonzeros(x::AbstractArray) = x mynnz(x) = length(mynonzeros(x)) function mul_rows_cols!( @@ -11,17 +11,3 @@ function mul_rows_cols!( B .= l .* A .* r' return nothing end - -function mul_rows_cols!( - B::SparseMatrixCSC, l::AbstractVector, A::SparseMatrixCSC, r::AbstractVector -) - @assert size(B) == size(A) == (length(l), length(r)) - @assert nnz(B) == nnz(A) - for j in axes(B, 2) - for k in nzrange(B, j) - i = B.rowval[k] - B.nzval[k] = l[i] * A.nzval[k] * r[j] - end - end - return nothing -end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 00000000..25c7c270 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,20 @@ +[deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LogarithmicNumbers = "aa2f6b4e-9042-5d33-9679-40d3a6b85899" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/test/correctness.jl b/test/correctness.jl index 9b92362f..5c245173 100644 --- a/test/correctness.jl +++ b/test/correctness.jl @@ -1,14 +1,8 @@ using Distributions -using HMMBase: HMMBase using HiddenMarkovModels import HiddenMarkovModels as HMMs -using HiddenMarkovModels: - LightDiagNormal, - LightCategorical, - rand_prob_vec, - rand_trans_mat, - test_equal_hmms, - test_coherent_algorithms +using HiddenMarkovModels: LightDiagNormal, LightCategorical +using HMMTest using LinearAlgebra using Random: Random, AbstractRNG, default_rng, seed! using SimpleUnPack @@ -18,200 +12,85 @@ using Test rng = default_rng() seed!(rng, 63) -function test_identical_hmmbase( - rng::AbstractRNG, - hmm::AbstractHMM, - hmm_guess::Union{Nothing,AbstractHMM}=nothing; - T::Integer, - atol::Real=1e-5, -) - sim = rand(rng, hmm, T) - obs_mat = collect(reduce(hcat, sim.obs_seq)') - - obs_seq = vcat(sim.obs_seq, sim.obs_seq) - seq_ends = [length(sim.obs_seq), 2 * length(sim.obs_seq)] - - hmm_base = HMMBase.HMM(hmm) - hmm_guess_base = HMMBase.HMM(hmm_guess) - - logL_base = HMMBase.forward(hmm_base, obs_mat)[2] - logL = logdensityof(hmm, obs_seq; seq_ends) - @test logL ≈ 2logL_base - - α_base, logL_forward_base = HMMBase.forward(hmm_base, obs_mat) - α, logL_forward = forward(hmm, obs_seq; seq_ends) - @test isapprox(α[:, 1:T], α_base') && isapprox(α[:, (T + 1):(2T)], α_base') - @test logL_forward ≈ 2logL_forward_base - - q_base = HMMBase.viterbi(hmm_base, obs_mat) - q, logL_viterbi = viterbi(hmm, obs_seq; seq_ends) - # Viterbi decoding can vary in case of (infrequent) ties - @test mean(q[1:T] .== q_base) > 0.9 && mean(q[(T + 1):(2T)] .== q_base) > 0.9 - - γ_base = HMMBase.posteriors(hmm_base, obs_mat) - γ, logL_forward_backward = forward_backward(hmm, obs_seq; seq_ends) - @test isapprox(γ[:, 1:T], γ_base') && isapprox(γ[:, (T + 1):(2T)], γ_base') - - if !isnothing(hmm_guess) - hmm_est_base, hist_base = HMMBase.fit_mle( - hmm_guess_base, obs_mat; maxiter=10, tol=-Inf - ) - logL_evolution_base = hist_base.logtots - hmm_est, logL_evolution = baum_welch( - hmm_guess, obs_seq; seq_ends, max_iterations=10, atol=-Inf - ) - @test isapprox( - logL_evolution[(begin + 1):end], 2 * logL_evolution_base[begin:(end - 1)] - ) - test_equal_hmms(hmm_est, HMM(hmm_est_base); atol, init=true) - end - - return nothing -end - ## Settings -T, K = 100, 20 +T, K = 200, 100 -## Distributions +init = [0.4, 0.6] +init_guess = [0.5, 0.5] -# TODO: add uniform +trans = [0.8 0.2; 0.2 0.8] +trans_guess = [0.7 0.3; 0.3 0.7] -@testset "Categorical" begin - init = [0.4, 0.6] - init_guess = [0.5, 0.5] +p = [[0.8, 0.2], [0.2, 0.8]] +p_guess = [[0.7, 0.3], [0.3, 0.7]] - trans = [0.8 0.2; 0.2 0.8] - trans_guess = [0.7 0.3; 0.3 0.7] +μ = [-ones(2), +ones(2)] +μ_guess = [-0.7 * ones(2), +0.7 * ones(2)] - dists = [Categorical([0.2, 0.8]), Categorical([0.8, 0.2])] - dists_guess = [Categorical([0.3, 0.7]), Categorical([0.7, 0.3])] +σ = ones(2) - hmm = HMM(init, trans, dists) - hmm_guess = HMM(init_guess, trans_guess, dists_guess) +control_seqs = [fill(nothing, rand(rng, T:(2T))) for k in 1:K]; +control_seq = reduce(vcat, control_seqs) +seq_ends = cumsum(length.(control_seqs)) - control_seq = fill(nothing, T * K) - seq_ends = T:T:(T * K) - test_identical_hmmbase(rng, hmm, hmm_guess; T) - test_coherent_algorithms( - rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false - ) -end +## Uncontrolled @testset "Normal" begin - init = [0.4, 0.6] - init_guess = [0.5, 0.5] - - trans = [0.8 0.2; 0.2 0.8] - trans_guess = [0.7 0.3; 0.3 0.7] - - dists = [Normal(-1.0), Normal(+1.0)] - dists_guess = [Normal(-0.7), Normal(+0.7)] + dists = [Normal(μ[1][1]), Normal(μ[2][1])] + dists_guess = [Normal(μ_guess[1][1]), Normal(μ_guess[2][1])] hmm = HMM(init, trans, dists) hmm_guess = HMM(init_guess, trans_guess, dists_guess) - control_seq = fill(nothing, T * K) - seq_ends = T:T:(T * K) test_identical_hmmbase(rng, hmm, hmm_guess; T) test_coherent_algorithms( rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false ) + test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) + test_allocations(rng, hmm, hmm_guess; control_seq, seq_ends) end @testset "DiagNormal" begin - init = [0.4, 0.6] - init_guess = [0.5, 0.5] - - trans = [0.8 0.2; 0.2 0.8] - trans_guess = [0.7 0.3; 0.3 0.7] - - D = 3 - dists = [MvNormal(-ones(D), I), MvNormal(+ones(D), I)] - dists_guess = [MvNormal(-ones(D) / 2, I), MvNormal(+ones(D) / 2, I)] + dists = [MvNormal(μ[1], I), MvNormal(μ[2], I)] + dists_guess = [MvNormal(μ_guess[1], I), MvNormal(μ_guess[2], I)] hmm = HMM(init, trans, dists) hmm_guess = HMM(init_guess, trans_guess, dists_guess) - control_seq = fill(nothing, T * K) - seq_ends = T:T:(T * K) test_identical_hmmbase(rng, hmm, hmm_guess; T) test_coherent_algorithms( rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false ) + test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) end -## Light distributions - @testset "LightCategorical" begin - init = [0.4, 0.6] - init_guess = [0.5, 0.5] - - trans = [0.8 0.2; 0.2 0.8] - trans_guess = [0.7 0.3; 0.3 0.7] - - dists = [LightCategorical([0.2, 0.8]), LightCategorical([0.8, 0.2])] - dists_guess = [LightCategorical([0.3, 0.7]), LightCategorical([0.7, 0.3])] + dists = [LightCategorical(p[1]), LightCategorical(p[2])] + dists_guess = [LightCategorical(p_guess[1]), LightCategorical(p_guess[2])] hmm = HMM(init, trans, dists) hmm_guess = HMM(init_guess, trans_guess, dists_guess) - control_seq = fill(nothing, T * K) - seq_ends = T:T:(T * K) test_coherent_algorithms( rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false ) + test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) + test_allocations(rng, hmm, hmm_guess; control_seq, seq_ends) end @testset "LightDiagNormal" begin - init = [0.4, 0.6] - init_guess = [0.5, 0.5] - - trans = [0.8 0.2; 0.2 0.8] - trans_guess = [0.7 0.3; 0.3 0.7] - - D = 3 - dists = [LightDiagNormal(-ones(D), ones(D)), LightDiagNormal(+ones(D), ones(D))] - dists_guess = [ - LightDiagNormal(-ones(D) / 2, ones(D)), LightDiagNormal(+ones(D) / 2, ones(D)) - ] - - hmm = HMM(init, trans, dists) - hmm_guess = HMM(init_guess, trans_guess, dists_guess) - - control_seq = fill(nothing, T * K) - seq_ends = T:T:(T * K) - test_coherent_algorithms( - rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false - ) -end + dists = [LightDiagNormal(μ[1], σ), LightDiagNormal(μ[2], σ)] + dists_guess = [LightDiagNormal(μ_guess[1], σ), LightDiagNormal(μ_guess[2], σ)] -## Weird arrays - -@testset "Normal sparse" begin - init = [0.2, 0.6, 0.2] - trans = sparse([ - 0.8 0.2 0.0 - 0.0 0.8 0.2 - 0.2 0.0 0.8 - ]) - dists = [Normal(-2.0), Normal(0.0), Normal(+2.0)] hmm = HMM(init, trans, dists) - - init_guess = [0.3, 0.4, 0.3] - trans_guess = sparse([ - 0.7 0.3 0.0 - 0.0 0.7 0.3 - 0.3 0.0 0.7 - ]) - dists_guess = [Normal(-1.5), Normal(0.0), Normal(+1.5)] hmm_guess = HMM(init_guess, trans_guess, dists_guess) - control_seq = fill(nothing, T * K) - seq_ends = T:T:(T * K) test_coherent_algorithms( rng, hmm, hmm_guess; control_seq, seq_ends, atol=0.05, init=false ) + test_type_stability(rng, hmm, hmm_guess; control_seq, seq_ends) + test_allocations(rng, hmm, hmm_guess; control_seq, seq_ends) end # Controlled @@ -236,12 +115,13 @@ function HMMs.obs_distributions(hmm::DiffusionHMM, λ::Number) end @testset "Controlled" begin - init = rand_prob_vec(rng, 2) - trans = rand_trans_mat(rng, 2) means = randn(rng, 2) hmm = DiffusionHMM(init, trans, means) - control_seq = rand(rng, T * K) - seq_ends = T:T:(T * K) + control_seqs = [[rand(rng) for t in 1:rand(T:(2T))] for k in 1:K] + control_seq = reduce(vcat, control_seqs) + seq_ends = cumsum(length.(control_seqs)) + test_coherent_algorithms(rng, hmm; control_seq, seq_ends, atol=0.05, init=false) + test_type_stability(rng, hmm; control_seq, seq_ends) end diff --git a/test/runtests.jl b/test/runtests.jl index 4b2fee65..9e1b3d36 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,10 +1,13 @@ using Aqua: Aqua using Documenter: Documenter using HiddenMarkovModels +using JET using JuliaFormatter: JuliaFormatter -using JET: JET +using Pkg using Test +Pkg.develop(; path=joinpath(dirname(@__DIR__), "libs", "HMMTest")) + @testset verbose = true "HiddenMarkovModels.jl" begin if VERSION >= v"1.9" @testset "Code formatting" begin @@ -18,13 +21,10 @@ using Test end @testset "Code linting" begin + using Distributions using Zygote JET.test_package(HiddenMarkovModels; target_defined_modules=true) end - - @testset "Types and allocations" begin - include("types_allocations.jl") - end end @testset "Distributions" begin diff --git a/test/types_allocations.jl b/test/types_allocations.jl deleted file mode 100644 index 3f3326d7..00000000 --- a/test/types_allocations.jl +++ /dev/null @@ -1,154 +0,0 @@ -using Distributions -using HiddenMarkovModels -import HiddenMarkovModels as HMMs -using HiddenMarkovModels: LightDiagNormal, LightCategorical, rand_prob_vec, rand_trans_mat -using JET: @test_opt, @test_call -using LinearAlgebra -using Random: AbstractRNG, default_rng, seed! -using SimpleUnPack -using SparseArrays -using Test - -rng = default_rng() -seed!(rng, 63) - -function test_type_stability( - rng::AbstractRNG, - hmm::AbstractHMM; - control_seq::AbstractVector, - seq_ends::AbstractVector{Int}=[length(control_seq)], -) - state_seq, obs_seq = rand(rng, hmm, control_seq) - - @test_opt target_modules = (HMMs,) rand(hmm, control_seq) - @test_call target_modules = (HMMs,) rand(hmm, control_seq) - - @test_opt target_modules = (HMMs,) logdensityof(hmm, obs_seq; control_seq, seq_ends) - @test_call target_modules = (HMMs,) logdensityof(hmm, obs_seq; control_seq, seq_ends) - @test_opt target_modules = (HMMs,) logdensityof( - hmm, obs_seq, state_seq; control_seq, seq_ends - ) - @test_call target_modules = (HMMs,) logdensityof( - hmm, obs_seq, state_seq; control_seq, seq_ends - ) - - @test_opt target_modules = (HMMs,) forward(hmm, obs_seq; control_seq, seq_ends) - @test_call target_modules = (HMMs,) forward(hmm, obs_seq; control_seq, seq_ends) - - @test_opt target_modules = (HMMs,) viterbi(hmm, obs_seq; control_seq, seq_ends) - @test_call target_modules = (HMMs,) viterbi(hmm, obs_seq; control_seq, seq_ends) - - @test_opt target_modules = (HMMs,) forward_backward(hmm, obs_seq; control_seq, seq_ends) - @test_call target_modules = (HMMs,) forward_backward( - hmm, obs_seq; control_seq, seq_ends - ) - - @test_opt target_modules = (HMMs,) baum_welch( - hmm, obs_seq; control_seq, seq_ends, max_iterations=1 - ) - @test_call target_modules = (HMMs,) baum_welch( - hmm, obs_seq; control_seq, seq_ends, max_iterations=1 - ) -end - -function test_allocations( - rng::AbstractRNG, - hmm::AbstractHMM; - control_seq::AbstractVector, - seq_ends::AbstractVector{Int}, -) - obs_seq = mapreduce(vcat, eachindex(seq_ends)) do k - t1, t2 = HMMs.seq_limits(seq_ends, k) - rand(rng, hmm, control_seq[t1:t2]).obs_seq - end - - ## Forward - forward(hmm, obs_seq; control_seq, seq_ends) # compile - f_storage = HMMs.initialize_forward(hmm, obs_seq; control_seq, seq_ends) - allocs = @allocated HiddenMarkovModels.forward!( - f_storage, hmm, obs_seq; control_seq, seq_ends - ) - @test allocs == 0 - - ## Viterbi - viterbi(hmm, obs_seq; control_seq, seq_ends) # compile - v_storage = HMMs.initialize_viterbi(hmm, obs_seq; control_seq, seq_ends) - allocs = @allocated HMMs.viterbi!(v_storage, hmm, obs_seq; control_seq, seq_ends) - @test allocs == 0 - - ## Forward-backward - forward_backward(hmm, obs_seq; control_seq, seq_ends) # compile - fb_storage = HMMs.initialize_forward_backward(hmm, obs_seq; control_seq, seq_ends) - allocs = @allocated HMMs.forward_backward!( - fb_storage, hmm, obs_seq; control_seq, seq_ends - ) - @test allocs == 0 - - ## Baum-Welch - baum_welch(hmm, obs_seq; control_seq, seq_ends, max_iterations=1) # compile - fb_storage = HMMs.initialize_forward_backward(hmm, obs_seq; control_seq, seq_ends) - logL_evolution = Float64[] - sizehint!(logL_evolution, 1) - hmm_guess = deepcopy(hmm) - allocs = @allocated HMMs.baum_welch!( - fb_storage, - logL_evolution, - hmm_guess, - obs_seq; - control_seq, - seq_ends, - atol=-Inf, - max_iterations=1, - loglikelihood_increasing=false, - ) - @test allocs == 0 -end - -N, D, T, K = 3, 2, 10, 4 -R = Float64 -control_seq = fill(nothing, K * T) -seq_ends = T:T:(K * T) -init = ones(R, N) / N -trans = ones(R, N, N) / N - -## Distributions - -@testset "Normal" begin - dists = [Normal(randn(rng, R), 1) for i in 1:N] - hmm = HMM(init, trans, dists) - test_type_stability(rng, hmm; control_seq, seq_ends) - test_allocations(rng, hmm; control_seq, seq_ends) -end - -@testset "DiagNormal" begin - dists = [MvNormal(randn(rng, R, D), I) for i in 1:N] - hmm = HMM(init, trans, dists) - test_type_stability(rng, hmm; control_seq, seq_ends) -end - -## Light distributions - -@testset "LightCategorical" begin - dists = [LightCategorical(rand_prob_vec(rng, R, D)) for i in 1:N] - hmm = HMM(init, trans, dists) - test_type_stability(rng, hmm; control_seq, seq_ends) - test_allocations(rng, hmm; control_seq, seq_ends) -end - -@testset "LightDiagNormal" begin - dists = [LightDiagNormal(randn(rng, R, D), ones(R, D)) for i in 1:N] - hmm = HMM(init, trans, dists) - test_type_stability(rng, hmm; control_seq, seq_ends) - test_allocations(rng, hmm; control_seq, seq_ends) -end - -## Weird arrays - -@testset "Normal sparse" begin - trans_sparse = sparse(SymTridiagonal(rand(rng, R, N), rand(rng, R, N - 1))) - foreach(HMMs.sum_to_one!, eachrow(trans_sparse)) - dists = [Normal(randn(rng, R), 1) for i in 1:N] - hmm = HMM(init, trans_sparse, dists) - test_type_stability(rng, hmm; control_seq, seq_ends) - test_allocations(rng, hmm; control_seq, seq_ends) -end