Manifest.toml

 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.4"
+julia_version = "1.10.5"
 manifest_format = "2.0"
 project_hash = "170a37ab5d68a18e29bd757512b09ed22a40cf6c"
 
@@ -16,39 +16,45 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[deps.CSV]]
 deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"]
-git-tree-sha1 = "44dbf560808d49041989b8a96cae4cffbeb7966a"
+git-tree-sha1 = "deddd8725e5e1cc49ee205a1964256043720a6c3"
 uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-version = "0.10.11"
+version = "0.10.15"
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "02aa26a4cf76381be7f66e020a3eddeb27b0a092"
+git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.2"
+version = "0.7.6"
 
 [[deps.ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
+git-tree-sha1 = "c7acce7a7e1078a20a285211dd73cd3941a871d6"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.11.4"
+version = "0.12.0"
+
+    [deps.ColorTypes.extensions]
+    StyledStringsExt = "StyledStrings"
+
+    [deps.ColorTypes.weakdeps]
+    StyledStrings = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
 
 [[deps.Colors]]
 deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
-git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a"
+git-tree-sha1 = "64e15186f0aa277e174aa81798f7eb8598e0157e"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.12.10"
+version = "0.13.0"
 
 [[deps.CommonSubexpressions]]
-deps = ["MacroTools", "Test"]
-git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+deps = ["MacroTools"]
+git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5"
 uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-version = "0.3.0"
+version = "0.3.1"
 
 [[deps.Compat]]
-deps = ["UUIDs"]
-git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7"
+deps = ["TOML", "UUIDs"]
+git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.9.0"
+version = "4.16.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -65,21 +71,21 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 version = "4.1.1"
 
 [[deps.DataAPI]]
-git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c"
+git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.15.0"
+version = "1.16.0"
 
 [[deps.DataFrames]]
-deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
-git-tree-sha1 = "04c738083f29f86e62c8afc341f0967d8717bdb8"
+deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
+git-tree-sha1 = "fb61b4812c49343d7ef0b533ba982c46021938a6"
 uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-version = "1.6.1"
+version = "1.7.0"
 
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d"
+git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.15"
+version = "0.18.20"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -115,41 +121,60 @@ version = "1.6.0"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "299dc33549f68299137e51e6d49a13b5b1da9673"
+git-tree-sha1 = "2dd20384bf8c6d411b5c7370865b1e9b26cb2ea3"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.1"
+version = "1.16.6"
+
+    [deps.FileIO.extensions]
+    HTTPExt = "HTTP"
+
+    [deps.FileIO.weakdeps]
+    HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 
 [[deps.FilePathsBase]]
-deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
-git-tree-sha1 = "e27c4ebe80e8699540f2d6c805cc12203b614f12"
+deps = ["Compat", "Dates"]
+git-tree-sha1 = "7878ff7172a8e6beedd1dea14bd27c3c6340d361"
 uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
-version = "0.9.20"
+version = "0.9.22"
+
+    [deps.FilePathsBase.extensions]
+    FilePathsBaseMmapExt = "Mmap"
+    FilePathsBaseTestExt = "Test"
+
+    [deps.FilePathsBase.weakdeps]
+    Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
+    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
-deps = ["LinearAlgebra", "Random"]
-git-tree-sha1 = "a20eaa3ad64254c61eeb5f230d9306e937405434"
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "6a70198746448456524cb442b8af316927ff3e1a"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.6.1"
-weakdeps = ["SparseArrays", "Statistics"]
+version = "1.13.0"
 
     [deps.FillArrays.extensions]
+    FillArraysPDMatsExt = "PDMats"
     FillArraysSparseArraysExt = "SparseArrays"
     FillArraysStatisticsExt = "Statistics"
 
+    [deps.FillArrays.weakdeps]
+    PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
-git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
-version = "0.8.4"
+version = "0.8.5"
 
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
-git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad"
+git-tree-sha1 = "a2df1b776752e3f344e5116c06d75a10436ab853"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.36"
+version = "0.10.38"
 weakdeps = ["StaticArrays"]
 
     [deps.ForwardDiff.extensions]
@@ -160,19 +185,26 @@ deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
 [[deps.InlineStrings]]
-deps = ["Parsers"]
-git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461"
+git-tree-sha1 = "45521d31238e87ee9f9732561bfee12d4eebd52d"
 uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
-version = "1.4.0"
+version = "1.4.2"
+
+    [deps.InlineStrings.extensions]
+    ArrowTypesExt = "ArrowTypes"
+    ParsersExt = "Parsers"
+
+    [deps.InlineStrings.weakdeps]
+    ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd"
+    Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.InvertedIndices]]
-git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
+git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe"
 uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
@@ -186,14 +218,14 @@ version = "1.0.0"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
+git-tree-sha1 = "a007feb38b422fbdab534406aeca1b86823cb4d6"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.5.0"
+version = "1.7.0"
 
 [[deps.LaTeXStrings]]
-git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
+git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-version = "1.3.0"
+version = "1.4.0"
 
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
@@ -224,15 +256,15 @@ uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
+git-tree-sha1 = "be484f5c92fad0bd8acfef35fe017900b0b73809"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.17.0+0"
+version = "1.18.0+0"
 
 [[deps.LightXML]]
 deps = ["Libdl", "XML2_jll"]
-git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
+git-tree-sha1 = "3a994404d3f6709610701c7dabfc03fed87a81f8"
 uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
-version = "0.9.0"
+version = "0.9.1"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
@@ -240,9 +272,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
+git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.26"
+version = "0.3.29"
 
     [deps.LogExpFunctions.extensions]
     LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -258,10 +290,9 @@ version = "0.3.26"
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[deps.MacroTools]]
-deps = ["Markdown", "Random"]
-git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
+git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.11"
+version = "0.5.15"
 
 [[deps.Markdown]]
 deps = ["Base64"]
@@ -274,9 +305,9 @@ version = "2.28.2+1"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
-git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"
+git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d"
 uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "1.1.0"
+version = "1.2.0"
 
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
@@ -287,9 +318,9 @@ version = "2023.1.10"
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
-git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
+git-tree-sha1 = "030ea22804ef91648f29b7ad3fc15fa49d0e6e71"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "1.0.2"
+version = "1.0.3"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -306,21 +337,21 @@ uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
 version = "0.8.1+2"
 
 [[deps.OpenSpecFun_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335"
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.5+0"
+version = "0.5.6+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3"
+git-tree-sha1 = "12f1439c4f986bb868acda6ea33ebc78e19b95ad"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.2"
+version = "1.7.0"
 
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
-git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851"
+git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.7.2"
+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"]
@@ -329,27 +360,27 @@ version = "1.10.0"
 
 [[deps.PooledArrays]]
 deps = ["DataAPI", "Future"]
-git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
+git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3"
 uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.0"
+version = "1.4.3"
 
 [[deps.PrettyTables]]
-deps = ["Crayons", "LaTeXStrings", "Markdown", "Printf", "Reexport", "StringManipulation", "Tables"]
-git-tree-sha1 = "ee094908d720185ddbdc58dbe0c1cbe35453ec7a"
+deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"]
+git-tree-sha1 = "1101cd475833706e4d0e7b122218257178f48f34"
 uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
-version = "2.2.7"
+version = "2.4.0"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -380,9 +411,9 @@ version = "0.7.0"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "04bdff0b09c65ff3e06a05e3eb7b120223da3d39"
+git-tree-sha1 = "712fb0231ee6f9120e005ccd56297abbc053e7e0"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.0"
+version = "1.4.8"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -392,9 +423,9 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
 [[deps.SortingAlgorithms]]
 deps = ["DataStructures"]
-git-tree-sha1 = "c60ec5c62180f27efea3ba2908480f8055e17cee"
+git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085"
 uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
-version = "1.1.1"
+version = "1.2.1"
 
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
@@ -403,9 +434,9 @@ version = "1.10.0"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
+git-tree-sha1 = "64cca0c26b4f31ba18f13f6c12af7c85f478cfde"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.3.1"
+version = "2.5.0"
 
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -414,19 +445,23 @@ version = "2.3.1"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 
 [[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
-git-tree-sha1 = "51621cca8651d9e334a659443a74ce50a3b6dfab"
+deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
+git-tree-sha1 = "47091a0340a675c738b1304b58161f3b0839d454"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.6.3"
-weakdeps = ["Statistics"]
+version = "1.9.10"
 
     [deps.StaticArrays.extensions]
+    StaticArraysChainRulesCoreExt = "ChainRulesCore"
     StaticArraysStatisticsExt = "Statistics"
 
+    [deps.StaticArrays.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -435,9 +470,9 @@ version = "1.10.0"
 
 [[deps.StringManipulation]]
 deps = ["PrecompileTools"]
-git-tree-sha1 = "a04cabe79c5f01f4d723cc6704070ada0b9d46d5"
+git-tree-sha1 = "a6b1675a536c5ad1a60e5a5153e1fee12eb146e3"
 uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e"
-version = "0.3.4"
+version = "0.4.0"
 
 [[deps.SuiteSparse_jll]]
 deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
@@ -456,25 +491,20 @@ uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
 version = "1.0.1"
 
 [[deps.Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
-git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec"
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"]
+git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.10.1"
+version = "1.12.0"
 
 [[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]]
-deps = ["Random", "Test"]
-git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769"
+git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.9.13"
+version = "0.11.3"
 
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
@@ -501,15 +531,15 @@ version = "1.6.1"
 
 [[deps.WriteVTK]]
 deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams", "VTKBase"]
-git-tree-sha1 = "7b46936613e41cfe1c6a5897d243ddcab8feabec"
+git-tree-sha1 = "1d8042d58334ab7947ce505709df7009da6f3375"
 uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192"
-version = "1.18.0"
+version = "1.21.1"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "04a51d15436a572301b5abbb9d099713327e9fc4"
+git-tree-sha1 = "a2fccc6559132927d4c5dc183e3e01048c6dcbd6"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.10.4+0"
+version = "2.13.5+0"
 
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
@@ -519,7 +549,7 @@ version = "1.2.13+1"
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.11.0+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]

scripts/run_experiments.jl

@@ -609,21 +609,45 @@ end
 
 huber(x, gamma) = abs(x) < gamma ? x^2 / (2 * gamma) : abs(x) - gamma / 2
 
-# TODO: finish!
-function refine_and_estimate_pd(st::L1L2TVState)
-    # globally refine
+function globally_refine!(hmesh::HMesh, st::L1L2TVState)
     marked_cells = Set(axes(st.mesh.cells, 2))
-    mesh_new, fs_new = refine(st.mesh, marked_cells;
-        st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
+    mesh_new, fs_new = refine!(hmesh, marked_cells;
+        st.tdata, st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
     st_new = L1L2TVState(st; mesh = mesh_new,
+        fs_new.tdata,
         fs_new.est, fs_new.g, fs_new.u, fs_new.p1, fs_new.p2,
         fs_new.du, fs_new.dp1, fs_new.dp2)
+    return st_new
+end
+
+function refine_and_estimate_pd!(st, ctx, i, interpolate_image_data!, f0 = nothing, f1 = nothing)
+    hmesh = HMesh(st.mesh)
+    hmesh_coarse_cells = cells(hmesh)
+
+    st_new = globally_refine!(hmesh, st)
 
     # compute primal dual error indicators
+    if isnothing(f0) && isnothing(f1)
+        interpolate_image_data!(st_new)
+    else
+        # TODO: we should maybe create a distinct state object for optflow and
+        # have f0, f1 be part of it to avoid this special casing.
+        interpolate_image_data!(f0, f1)
+    end
     estimate_pd!(st_new)
 
+    # Output optional fine data
+    if i == 1
+        # Hack to get numbering from 0 as the other vtu data.
+        output(st_new, joinpath(ctx.outdir, "output_pd_fine_$(lpad(0, 5, '0')).vtu"), st_new.g, st_new.est)
+    end
+    output(st_new, joinpath(ctx.outdir, "output_pd_fine_$(lpad(i, 5, '0')).vtu"), st_new.g, st_new.est)
+
     # transfer data to old state
-    #
+    mesh_coarse, fs_coarse = coarsen(hmesh, hmesh_coarse_cells; st_new.est)
+    @assert cells(st.mesh) == cells(mesh_coarse)
+
+    st.est.data .= fs_coarse.est.data
 end
 
 # this computes the primal-dual error indicator which is not really useful
@@ -1359,7 +1383,7 @@ function inpaint(ctx)
         ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
 
-    function interpolate_image_data!()
+    function interpolate_image_data!(st)
         println("interpolate image data ...")
         if ctx.params.n_refine == 0
             # if we use a grid at image resolution, we wan't to avoid
@@ -1381,7 +1405,7 @@ function inpaint(ctx)
 
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
 
-        interpolate_image_data!()
+        interpolate_image_data!(st)
 
         pvd[0] = save_step(0)
 
@@ -1403,7 +1427,11 @@ function inpaint(ctx)
             # plot
             i += 1
             display(plot(grayclamp.(to_img(sample(st.u)))))
-            estimate_res!(st)
+            if hasproperty(ctx.params, :estimator) && ctx.params.estimator == :primal_dual
+                refine_and_estimate_pd!(st, ctx, i, interpolate_image_data!)
+            else
+                estimate_res!(st)
+            end
             pvd[i] = save_step(i)
             #break
 
@@ -1424,11 +1452,15 @@ function inpaint(ctx)
 
             # manually mark all cell within inpainting domain, since the
             # estimator is not reliable there
-            for cell in cells(mesh)
-                bind!(st.tdata, cell)
-                maskv = evaluate(st.tdata, SA[1/3, 1/3])
-                if abs(maskv[begin] - 1.) > 1e-8
-                    push!(marked_cells, cell)
+            if !hasproperty(ctx.params, :estimator) || ctx.params.estimator != :primal_dual
+                # The primal-dual estimator seems reliable inside the
+                # inpainting domain
+                for cell in cells(mesh)
+                    bind!(st.tdata, cell)
+                    maskv = evaluate(st.tdata, SA[1/3, 1/3])
+                    if abs(maskv[begin] - 1.) > 1e-8
+                        push!(marked_cells, cell)
+                    end
                 end
             end
             #marked_cells = Set(axes(mesh.cells, 2))
@@ -1436,7 +1468,7 @@ function inpaint(ctx)
                 st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2, st.tdata)
             st = L1L2TVState(st; mesh, fs.tdata,
                 fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
-            interpolate_image_data!()
+            interpolate_image_data!(st)
         end
     end
     end # @elapsed
@@ -1462,6 +1494,7 @@ function experiment_inpaint_adaptive(ctx)
         params = (
             name = "test",
             n_refine = coarsen,
+            #estimator = :primal_dual,
             g_arr, mask_arr, mesh,
             #alpha1 = 0.2, alpha2 = 8., lambda = 1.,
             alpha1 = 0., alpha2 = 50., lambda = 1.,
@@ -1582,7 +1615,7 @@ function optflow(ctx)
         interpolate!(st.g, g_optflow; u0 = st.u, f0, fw, st.tdata)
     end
 
-    function interpolate_image_data!()
+    function interpolate_image_data!(f0, f1)
         println("interpolate image data ...")
         project_image!(f0, imgf0)
         project_image!(f1, imgf1)
@@ -1596,7 +1629,7 @@ function optflow(ctx)
 
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
 
-        interpolate_image_data!()
+        interpolate_image_data!(f0, f1)
         warp!() # in the first step just to fill st.g
         pvd[0] = save_step(0)
 
@@ -1640,7 +1673,11 @@ function optflow(ctx)
             k_refine += 1
             k_refine > ctx.params.n_refine && break
             println("refine ...")
-            estimate_res!(st)
+            if hasproperty(ctx.params, :estimator) && ctx.params.estimator == :primal_dual
+                refine_and_estimate_pd!(st, ctx, i, interpolate_image_data!, f0, f1)
+            else
+                estimate_res!(st)
+            end
             marked_cells = mark(st; theta = 0.5)
             #marked_cells = Set(axes(mesh.cells, 2))
             mesh, fs = refine(mesh, marked_cells;
@@ -1649,7 +1686,7 @@ function optflow(ctx)
             st = L1L2TVState(st; mesh, fs.tdata,
                 fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
             f0, f1, fw = (fs.f0, fs.f1, fs.fw)
-            interpolate_image_data!()
+            interpolate_image_data!(f0, f1)
         end
     end
     #CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
@@ -1684,7 +1721,7 @@ function experiment_optflow_middlebury(ctx)
     imgf0 = loadimg(joinpath(ctx.indir, "frame10.png"))
     imgf1 = loadimg(joinpath(ctx.indir, "frame11.png"))
     gtflow = FileIO.load(joinpath(ctx.indir, "flow10.flo"))
-    maxflow = OpticalFlowUtils._maxflow(gtflow)
+    maxflow = OpticalFlowUtils.maxflow(gtflow)
 
     ctx = Util.Context(ctx; imgf0, imgf1, maxflow, gtflow)
     saveimg(joinpath(ctx.outdir, "ground_truth.png"), colorflow(gtflow; maxflow))
@@ -1700,12 +1737,6 @@ function experiment_optflow_middlebury_all_benchmarks(ctx)
     end
 end
 
-# FIXME: legacy
-function experiment_optflow_middlebury_all(ctx)
-    experiment_optflow_middlebury_all_benchmarks(
-        Util.Context(ctx; warp = true, refine = true, n_refine = 6))
-end
-
 function experiment_optflow_middlebury_warping_comparison(ctx)
     ctx(experiment_optflow_middlebury_all_benchmarks, "vanilla";
         warp = false, refine = false, n_refine = 0)
@@ -1720,6 +1751,8 @@ function experiment_optflow_middlebury_warping_comparison_adaptive(ctx)
         warp = true, refine = false, n_refine = 0)
     ctx(experiment_optflow_middlebury_all_benchmarks, "adaptive-warping";
         warp = true, refine = true, n_refine = 6)
+    #ctx(experiment_optflow_middlebury_all_benchmarks, "adaptive-warping";
+    #    warp = true, refine = true, n_refine = 6, estimator = :residual)
 end
 
 function experiment_optflow_schoice(ctx)

src/function.jl

@@ -24,7 +24,7 @@ evaluate_basis(::DP1, x) = evaluate_basis(P1(), x)
 struct FeSpace{M, Fe, S}
     mesh::M
     element::Fe
-    dofmap::Array{Int, 3} # (rdim, eldof, cell) -> gdof
+    dofmap::Array{Int, 3} # (rangedim, eldof, cell) -> gdof
     ndofs::Int # = maximum(dofmap)
 end
 
@@ -64,7 +64,7 @@ function FeSpace(mesh, el::DP1, size_=(1,))
 end
 
 Base.show(io::IO, ::MIME"text/plain", x::FeSpace) =
-    print("$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(ndofs(x)) dofs")
+    print(io, "$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(ndofs(x)) dofs")
 
 function Base.getproperty(obj::FeSpace{<:Any, <:Any, S}, sym::Symbol) where S
     if sym === :size
@@ -104,7 +104,7 @@ function FeFunction(space; name=string(gensym("f")))
 end
 
 Base.show(io::IO, ::MIME"text/plain", f::FeFunction) =
-    print("$(nameof(typeof(f))), size $(f.space.size) with $(length(f.data)) dofs")
+    print(io, "$(nameof(typeof(f))), size $(f.space.size) with $(length(f.data)) dofs")
 
 interpolate!(dst::FeFunction, expr::Function; params...) =
     interpolate!(dst, dst.space.element, expr; params...)
@@ -256,7 +256,7 @@ struct Derivative{F, M}
 end
 
 Base.show(io::IO, ::MIME"text/plain", x::Derivative) =
-    print("$(nameof(typeof(x))) of $(typeof(x.f))")
+    print(io, "$(nameof(typeof(x))) of $(typeof(x.f))")
 
 # TODO: should be called "derivative"
 function nabla(f)
@@ -380,6 +380,9 @@ function prolong!(new_f, old_f, old_cell, new_cells, ::DP0)
 end
 
 
+# IO
+
+
 vtk_append!(vtkfile, f::FeFunction, fs::FeFunction...) =
     (vtk_append!(vtkfile, f); vtk_append!(vtkfile, fs...))
 

src/mesh.jl

-export init_grid, init_hgrid, save, refine!, refine, cells, vertices,
-    ndims_domain, ndims_space, nvertices, nvertices_cell, ncells
+export HMesh, init_grid, init_hgrid, save, refine!, refine, cells, vertices,
+    ndims_domain, ndims_space, nvertices, nvertices_cell, ncells, coarsen
 
 using LinearAlgebra: norm
 
@@ -9,12 +9,14 @@ using StaticArrays: SVector
 struct HMesh
     vertices::Vector{SVector{2, Float64}}
     cells::Vector{NTuple{3, Int}}
+    "direct descendents of a cell"
+    children::Vector{Union{NTuple{2, Int}, Nothing}}
     "refinement depth below each cell"
     levels::Vector{Int}
 end
 
 Base.show(io::IO, ::MIME"text/plain", x::HMesh) =
-    print("$(nameof(typeof(x))), $(ncells(x)) cells")
+    print(io, "$(nameof(typeof(x))), $(ncells(x)) cells")
 
 ndims_domain(::HMesh) = 2
 ndims_space(::HMesh) = 2
@@ -44,8 +46,10 @@ function init_hgrid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
 	push!(cells, (vidx[I], vidx[I + e1], vidx[I + e1 + e2]))
 	push!(cells, (vidx[I], vidx[I + e1 + e2], vidx[I + e2]))
     end
+    children = fill(nothing, axes(cells))
+    levels = zeros(Int, axes(cells))
 
-    return HMesh(vertices, cells, zeros(Int, axes(cells)))
+    return HMesh(vertices, cells, children, levels)
 end
 
 function init_hgrid(img::Array{<:Any, 2}; type=:vertex)
@@ -79,13 +83,14 @@ function HMesh(mesh::Mesh)
     for c in axes(mesh.cells, 2)
 	push!(cells, NTuple{3}(mesh.cells[:, c]))
     end
+    children = fill(nothing, axes(cells))
     levels = zeros(Int, axes(cells))
 
-    return HMesh(vertices, cells, levels)
+    return HMesh(vertices, cells, children, levels)
 end
 
 Base.show(io::IO, ::MIME"text/plain", x::Mesh) =
-    print("$(nameof(typeof(x))), $(ncells(x)) cells")
+    print(io, "$(nameof(typeof(x))), $(ncells(x)) cells")
 
 ndims_domain(::Mesh) = 2
 ndims_space(::Mesh) = 2
@@ -156,10 +161,10 @@ function init_grid(a::Array{<:Any, 2}, m::Int, n::Int = m; type=:vertex)
         init_grid((m, n)..., (0.5, 0.5), s .- (0.5, 0.5))
 end
 
-# horribly implemented, please don't curse me
+# TODO: cleanup
 function bisect!(mesh::HMesh, marked_cells::Set)
     refined_cells = Pair{Int, NTuple{2, Int}}[]
-    # assemble edge -> cells map
+    # assemble mapping: edge -> cells
     edgemap = Dict{NTuple{2, Int}, Vector{Int}}()
     for cell in cells(mesh)
 	vs = sort(SVector(vertices(mesh, cell)))
@@ -173,8 +178,6 @@ function bisect!(mesh::HMesh, marked_cells::Set)
 	edgemap[e3] = push!(get!(edgemap, e3, []), cell)
     end
 
-    #return edgemap
-
     function refine_cell(c1)
 	c2 = -1
 	# c1 -> c11 + c12
@@ -208,6 +211,7 @@ function bisect!(mesh::HMesh, marked_cells::Set)
 	# take care to produce positively oriented cells
 	push!(mesh.cells, Tuple(replace(SVector(vertices(mesh, c1)), c1_vs[1] => vbisect)))
 	push!(mesh.levels, 0)
+	push!(mesh.children, nothing)
 	c3 = lastindex(mesh.cells)
 	@assert(length(setdiff(c1_vs, c1_vs[1])) == 2)
 	replace!(edgemap[NTuple{2}(setdiff(c1_vs, c1_vs[1]))], c1 => c3)
@@ -215,6 +219,7 @@ function bisect!(mesh::HMesh, marked_cells::Set)
 	# take care to produce positively oriented cells
 	push!(mesh.cells, Tuple(replace(SVector(vertices(mesh, c1)), c1_vs[2] => vbisect)))
 	push!(mesh.levels, 0)
+	push!(mesh.children, nothing)
 	c4 = lastindex(mesh.cells)
 	@assert(length(setdiff(c1_vs, c1_vs[2])) == 2)
 	replace!(edgemap[NTuple{2}(setdiff(c1_vs, c1_vs[2]))], c1 => c4)
@@ -227,11 +232,13 @@ function bisect!(mesh::HMesh, marked_cells::Set)
 	mesh.levels[c1] += 1
 	delete!(marked_cells, c1)
 	push!(refined_cells, c1 => (c3, c4))
+	mesh.children[c1] = (c3, c4)
 
 	if c2 > 0
 	    # take care to produce positively oriented cells
 	    push!(mesh.cells, Tuple(replace(SVector(vertices(mesh, c2)), c1_vs[1] => vbisect))) # c1_vs is correct
 	    push!(mesh.levels, 0)
+	    push!(mesh.children, nothing)
 	    c5 = lastindex(mesh.cells)
 	    @assert(length(setdiff(c2_vs, c1_vs[1])) == 2)
 	    replace!(edgemap[NTuple{2}(setdiff(c2_vs, c1_vs[1]))], c2 => c5)
@@ -239,6 +246,7 @@ function bisect!(mesh::HMesh, marked_cells::Set)
 	    # take care to produce positively oriented cells
 	    push!(mesh.cells, Tuple(replace(SVector(vertices(mesh, c2)), c1_vs[2] => vbisect))) # c1_vs is correct
 	    push!(mesh.levels, 0)
+	    push!(mesh.children, nothing)
 	    c6 = lastindex(mesh.cells)
 	    @assert(length(setdiff(c2_vs, c1_vs[2])) == 2)
 	    replace!(edgemap[NTuple{2}(setdiff(c2_vs, c1_vs[2]))], c2 => c6)
@@ -251,6 +259,7 @@ function bisect!(mesh::HMesh, marked_cells::Set)
 	    mesh.levels[c2] += 1
 	    delete!(marked_cells, c2)
 	    push!(refined_cells, c2 => (c5, c6))
+	    mesh.children[c2] = (c5, c6)
 	end
     end
 
@@ -265,46 +274,64 @@ end
 # TODO: cleanup and optimize
 # FIXME: this assumes hmesh was created from f.mesh
 function refine!(hmesh::HMesh, marked_cells::Set; fs...)
+    fs = NamedTuple(fs)
+
     old_mesh = sub_mesh(hmesh)
+    hmesh_old_leaf_cells = cells(hmesh)
+    @assert allequal(x -> x.space.mesh, fs) "Functions living on different meshes."
+    for f in fs
+        # Note this is only a superficial check. We don't spend the effort to
+        # check that meshes match structurally as well.
+        @assert length(cells(f.space.mesh)) == length(hmesh_old_leaf_cells) "Mismatching number of cells for mesh function and lowest layer of hierarchical mesh: $(length(cells(f.space.mesh))) vs $(length(hmesh_leaf_cells))."
+    end
+
     cell_refinements = bisect!(hmesh, marked_cells)
-    # we use a extended mesh containing all cells (both fine and
-    # coarse) since a cell might have been refinemened multiple times and the
-    # intermediate cells will not be present in the final mesh for direct
-    # prolongation of data
+    hmesh_new_leaf_cells = cells(hmesh)
+
+    # We use an extended mesh containing all cells and build functions on top
+    # of that since prolongation across multiple levels becomes more convenient
+    # when having a single set of cells to index with `cell_refinements`.
     extended_mesh = sub_mesh(hmesh, axes(hmesh.cells, 1))
+    # Final leaf mesh, reindexed.
     new_mesh = sub_mesh(hmesh)
+
     refined_cells = map(x -> first(x), cell_refinements)
 
-    # extended functions onto newly created cells
-    extended_fs = map(NamedTuple(fs)) do f
+    # create extended functions on all cells of extended mesh
+    extended_fs = map(fs) do f
 	space = FeSpace(extended_mesh, f.space.element, f.space.size)
 	return FeFunction(space; f.name)
     end
-    # copy over previous data for unmodified cells
-    for (f, extended_f) in zip(NamedTuple(fs), extended_fs)
-	copyto!(extended_f.data, f.data)
+    # copy over data for unmodified cells
+    for (cell, h_cell) in enumerate(hmesh_old_leaf_cells)
+	for (ext_f, f) in zip(extended_fs, fs)
+	    gdofs = f.space.dofmap[:, :, cell]
+	    ext_gdofs = ext_f.space.dofmap[:, :, h_cell]
+	    ext_f.data[ext_gdofs] .= f.data[gdofs]
+	end
     end
     # prolong data for refined cells
+    # TODO: use hmesh.children instead of `cell_refinements`.
     for (old_cell, extended_cells) in cell_refinements
 	for extended_f in extended_fs
 	    prolong!(extended_f, extended_f, old_cell, extended_cells)
 	end
     end
 
-    # retain only non-refined cells
-    new_fs = map(NamedTuple(extended_fs)) do f
+    # create final functions on the new leaf mesh
+    new_fs = map(extended_fs) do f
 	space = FeSpace(new_mesh, f.space.element, f.space.size)
 	return FeFunction(space; f.name)
     end
-    retained_cells = setdiff(cells(extended_mesh), refined_cells)
-    @assert(retained_cells == cells(hmesh))
-    for (new_cell, old_cell) in enumerate(retained_cells)
+    # copy over all leaf cell data into newly indexed mesh
+    for (new_cell, old_cell) in enumerate(hmesh_new_leaf_cells)
 	for (f, new_f) in zip(extended_fs, new_fs)
 	    gdofs = f.space.dofmap[:, :, old_cell]
 	    new_gdofs = new_f.space.dofmap[:, :, new_cell]
 	    new_f.data[new_gdofs] .= f.data[gdofs]
 	end
     end
+
     return new_mesh, new_fs
 end
 
@@ -314,6 +341,47 @@ function refine(mesh::Mesh, marked_cells; fs...)
     return refine!(hmesh, Set(marked_cells); fs...)
 end
 
+"""
+Coarsen functions to a new coarse mesh.
+
+Assumptions:
+  - `hmesh` contains the coarse mesh with the given cell indices,
+  - The meshes of `fs...` are exactly the finest level of `hmesh`.
+"""
+function coarsen(hmesh::HMesh, hmesh_coarse_cells; fs...)
+    fs = NamedTuple(fs)
+    coarse_mesh = sub_mesh(hmesh, hmesh_coarse_cells)
+    hmesh_fine_cells = cells(hmesh)
+
+    coarse_fs = map(fs) do f
+        @assert length(cells(f.space.mesh)) == count(==(0), hmesh.levels)
+	space = FeSpace(coarse_mesh, f.space.element, f.space.size)
+	return FeFunction(space; f.name)
+    end
+
+    function coarse_dofs(hcell, f)
+        if hmesh.levels[hcell] == 0
+            # reverse cell index map
+            fine_cell = findfirst(==(hcell), hmesh_fine_cells)
+            @assert f.space.element isa DP0 "only DP0 implemented"
+            gdofs = f.space.dofmap[:, 1, fine_cell]
+            return f.data[gdofs]
+        else
+            return sum(coarse_dofs(subhcell, f) for subhcell in hmesh.children[hcell]) / length(hmesh.children[hcell])
+        end
+    end
+
+    for coarse_cell in cells(coarse_mesh)
+        for (f, coarse_f) in zip(fs, coarse_fs)
+            # need to map to `hmesh` indexing
+            coarse_gdofs = coarse_f.space.dofmap[:, 1, coarse_cell]
+            coarse_f.data[coarse_gdofs] .= coarse_dofs(hmesh_coarse_cells[coarse_cell], f)
+        end
+    end
+
+    return coarse_mesh, coarse_fs
+end
+
 function geo_tolocal(A, v)
     J = jacobian(x -> A * [1 - x[1] - x[2], x[1], x[2]], [0., 0.])
     return J \ (v - A[:, 1])

test/Manifest.toml

 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.0"
+julia_version = "1.10.5"
 manifest_format = "2.0"
+project_hash = "a87db4a96b9f8108c69768cc4f6470d5533e4e6e"
 
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
@@ -18,11 +19,12 @@ deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[deps.Random]]
-deps = ["SHA", "Serialization"]
+deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

test/runtests.jl

@@ -53,23 +53,38 @@ end
     @test img == img_sampled
 end
 
-@testset "mesh refinement" begin
+@testset "hierarchical mesh" begin
     # simple hashing
     val(f) = integrate(f.space.mesh, (x; f) -> sum(f .* f); f)
 
     mesh = init_grid(zeros(5, 5))
+    hmesh = HMesh(mesh)
 
-    f = FeFunction(FeSpace(mesh, P1(), (2, 3)))
+    f = FeFunction(FeSpace(mesh, DP0(), (2, 3)))
     f.data .= rand(size(f.data)...)
 
-    for i = 1:2
-        # refine roughly half of all cells
-        cells_ref = Set(findall(_ -> rand(Bool), cells(mesh)))
+    hmesh_cells_per_refinement = Vector{Int}[]
+    n = 5
 
-        mesh_new, (f_new,) = refine(mesh, cells_ref; f)
+    @testset "refinement $i" for i in 1:n
+        push!(hmesh_cells_per_refinement, cells(hmesh))
+        # Random set of cells on the finest level.
+        cells_ref = Set(filter(_ -> rand(Bool), cells(hmesh)))
+        mesh_new, (f_new,) = refine!(hmesh, cells_ref; f)
 
-        @test isapprox(val(f), val(f_new))
+        @test all(isnothing, hmesh.children[iszero.(hmesh.levels)])
+        @test all(!isnothing, hmesh.children[hmesh.levels .> 0])
+
+        @test isapprox(val(f_new), val(f))
 
         (mesh, f) = (mesh_new, f_new)
     end
+
+    @testset "coarsening" for i in n:-1:1
+        (mesh_coarse, (f_coarse,)) = coarsen(hmesh, hmesh_cells_per_refinement[i]; f)
+
+        @test isapprox(val(f_coarse), val(f))
+
+        #(mesh, f) = (mesh_coarse, f_coarse)
+    end
 end