use MKLDirect for G₂

IntervalMatrices use Rump algorithm to matrix multiplication
This brings time to Qint'*Qint down to ~40s which is
5-8 × slower than Q'*Q (for size n=2^13).

The naive version is ~100 × slower than Q'*Q even for n = 2^10.

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "5452e4eb8836d03e65c342a444826b894c3e14dc"
+project_hash = "46bad61f58c1a6b455db53aa91fae7351cb279ab"
 
 [[deps.AbstractAlgebra]]
 deps = ["GroupsCore", "InteractiveUtils", "LinearAlgebra", "MacroTools", "Markdown", "Random", "RandomExtensions", "SparseArrays", "Test"]
@@ -10,23 +10,6 @@ git-tree-sha1 = "29e65c331f97db9189ef00a4c7aed8127c2fd2d4"
 uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
 version = "0.27.10"
 
-[[deps.Accessors]]
-deps = ["Compat", "CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Requires", "StaticArrays", "Test"]
-git-tree-sha1 = "beabc31fa319f9de4d16372bff31b4801e43d32c"
-uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.28"
-
-[[deps.Adapt]]
-deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef"
-uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.6.1"
-
-[[deps.ArgCheck]]
-git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
-uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
-version = "2.3.0"
-
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 version = "1.1.1"
@@ -34,20 +17,9 @@ version = "1.1.1"
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
-[[deps.BangBang]]
-deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
-git-tree-sha1 = "7fe6d92c4f281cf4ca6f2fba0ce7b299742da7ca"
-uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
-version = "0.3.37"
-
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
-[[deps.Baselet]]
-git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
-uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
-version = "0.1.1"
-
 [[deps.BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
 git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8"
@@ -113,48 +85,22 @@ deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 version = "1.0.1+0"
 
-[[deps.CompositionsBase]]
-git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
-uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
-version = "0.1.1"
-
-[[deps.ConstructionBase]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "89a9db8d28102b094992472d333674bd1a83ce2a"
-uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.5.1"
-
 [[deps.Cyclotomics]]
 deps = ["LRUCache", "Memoize", "Primes", "SparseArrays", "Test"]
 git-tree-sha1 = "dc2e5fd64c188399434e83fa5c10c6fa4eff962a"
 uuid = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 version = "0.3.2"
 
-[[deps.DataAPI]]
-git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630"
-uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.14.0"
-
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
 git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 version = "0.18.13"
 
-[[deps.DataValueInterfaces]]
-git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
-uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
-version = "1.0.0"
-
 [[deps.Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
-[[deps.DefineSingletons]]
-git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
-uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
-version = "0.1.2"
-
 [[deps.DiffResults]]
 deps = ["StaticArraysCore"]
 git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
@@ -187,11 +133,6 @@ git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
 uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
 version = "0.5.2"
 
-[[deps.ExternalDocstrings]]
-git-tree-sha1 = "1224740fc4d07c989949e1c1b508ebd49a65a5f6"
-uuid = "e189563c-0753-4f5e-ad5c-be4293c83fb4"
-version = "0.1.1"
-
 [[deps.FastRounding]]
 deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
 git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
@@ -201,27 +142,17 @@ version = "0.3.1"
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
-[[deps.Folds]]
-deps = ["Accessors", "BangBang", "Baselet", "DefineSingletons", "Distributed", "ExternalDocstrings", "InitialValues", "MicroCollections", "Referenceables", "Requires", "Test", "ThreadedScans", "Transducers"]
-git-tree-sha1 = "638109532de382a1f99b1aae1ca8b5d08515d85a"
-uuid = "41a02a25-b8f0-4f67-bc48-60067656b558"
-version = "0.2.8"
-
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
 git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
 version = "0.10.35"
 
-[[deps.Future]]
-deps = ["Random"]
-uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
-
 [[deps.Groups]]
-deps = ["Folds", "GroupsCore", "KnuthBendix", "LinearAlgebra", "Logging", "OrderedCollections", "PermutationGroups", "StaticArrays"]
-git-tree-sha1 = "5dbf642ee0048e6ad5f0bda11af17e40b8e8dd2f"
+deps = ["GroupsCore", "KnuthBendix", "LinearAlgebra", "Logging", "OrderedCollections", "PermutationGroups", "Random", "StaticArrays"]
+git-tree-sha1 = "47967c88672431b6d3624d1cb200abdc15aa0b47"
 uuid = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
-version = "0.7.5"
+version = "0.7.6"
 
 [[deps.GroupsCore]]
 deps = ["Markdown", "Random"]
@@ -229,11 +160,6 @@ git-tree-sha1 = "9e1a5e9f3b81ad6a5c613d181664a0efc6fe6dd7"
 uuid = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 version = "0.4.0"
 
-[[deps.InitialValues]]
-git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3"
-uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c"
-version = "0.3.1"
-
 [[deps.IntegerMathUtils]]
 git-tree-sha1 = "f366daebdfb079fd1fe4e3d560f99a0c892e15bc"
 uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
@@ -249,6 +175,12 @@ git-tree-sha1 = "c1c88395d09366dae431556bcb598ad08fa1392b"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 version = "0.20.8"
 
+[[deps.IntervalMatrices]]
+deps = ["IntervalArithmetic", "LinearAlgebra", "Random", "Reexport", "SparseArrays"]
+git-tree-sha1 = "b5a981918e38fb40aeea659ec6511187fa488c14"
+uuid = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
+version = "0.8.3"
+
 [[deps.InverseFunctions]]
 deps = ["Test"]
 git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f"
@@ -260,11 +192,6 @@ git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
 version = "0.2.2"
 
-[[deps.IteratorInterfaceExtensions]]
-git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
-uuid = "82899510-4779-5014-852e-03e436cf321d"
-version = "1.0.0"
-
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
 git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
@@ -279,9 +206,9 @@ version = "0.21.3"
 
 [[deps.JuMP]]
 deps = ["LinearAlgebra", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays"]
-git-tree-sha1 = "611b9f12f02c587d860c813743e6cec6264e94d8"
+git-tree-sha1 = "4ec0e68fecbbe1b78db2ddf1ac573963ed5adebc"
 uuid = "4076af6c-e467-56ae-b986-b466b2749572"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.KnuthBendix]]
 deps = ["MacroTools", "ProgressMeter"]
@@ -341,9 +268,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[deps.MathOptInterface]]
 deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
-git-tree-sha1 = "f219b62e601c2f2e8adb7b6c48db8a9caf381c82"
+git-tree-sha1 = "88551cb56065a7bf9276af1d1970d2bd0d4b6ed4"
 uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-version = "1.13.1"
+version = "1.14.0"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -356,12 +283,6 @@ git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
 uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
 version = "0.4.4"
 
-[[deps.MicroCollections]]
-deps = ["BangBang", "InitialValues", "Setfield"]
-git-tree-sha1 = "629afd7d10dbc6935ec59b32daeb33bc4460a42e"
-uuid = "128add7d-3638-4c79-886c-908ea0c25c34"
-version = "0.1.4"
-
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
@@ -402,9 +323,9 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.4.1"
+version = "1.6.0"
 
 [[deps.Parsers]]
 deps = ["Dates", "SnoopPrecompile"]
@@ -469,17 +390,10 @@ git-tree-sha1 = "261dddd3b862bd2c940cf6ca4d1c8fe593e457c8"
 uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 version = "1.3.3"
 
-[[deps.Referenceables]]
-deps = ["Adapt"]
-git-tree-sha1 = "e681d3bfa49cd46c3c161505caddf20f0e62aaa9"
-uuid = "42d2dcc6-99eb-4e98-b66c-637b7d73030e"
-version = "0.1.2"
-
-[[deps.Requires]]
-deps = ["UUIDs"]
-git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
-uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.3.0"
+[[deps.Reexport]]
+git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.2.2"
 
 [[deps.RoundingEmulator]]
 git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
@@ -498,12 +412,6 @@ git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
 uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
 version = "0.2.1"
 
-[[deps.Setfield]]
-deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
-git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
-uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-version = "1.1.1"
-
 [[deps.SnoopPrecompile]]
 deps = ["Preferences"]
 git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
@@ -523,12 +431,6 @@ git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.2.0"
 
-[[deps.SplittablesBase]]
-deps = ["Setfield", "Test"]
-git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5"
-uuid = "171d559e-b47b-412a-8079-5efa626c420e"
-version = "0.1.15"
-
 [[deps.StarAlgebras]]
 deps = ["LinearAlgebra", "SparseArrays"]
 git-tree-sha1 = "265b89a5dfb38fe94ad48b997a253b2393fce6f1"
@@ -537,9 +439,9 @@ version = "0.2.0"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
-git-tree-sha1 = "6aa098ef1012364f2ede6b17bf358c7f1fbe90d4"
+git-tree-sha1 = "b8d897fe7fa688e93aef573711cb207c08c9e11e"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.5.17"
+version = "1.5.19"
 
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
@@ -552,7 +454,7 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.SymbolicWedderburn]]
 deps = ["Cyclotomics", "GroupsCore", "LinearAlgebra", "PermutationGroups", "Primes", "SparseArrays", "StarAlgebras"]
-git-tree-sha1 = "3a098aae261c4c57d6ee41dab62973065f0cb68b"
+git-tree-sha1 = "7afdc5576b8fc49032c1f765c1514e5369fb77fe"
 repo-rev = "enh/towards_simple_projections"
 repo-url = "https://github.com/kalmarek/SymbolicWedderburn.jl.git"
 uuid = "858aa9a9-4c7c-4c62-b466-2421203962a2"
@@ -563,18 +465,6 @@ deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 version = "1.0.0"
 
-[[deps.TableTraits]]
-deps = ["IteratorInterfaceExtensions"]
-git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
-uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
-version = "1.0.1"
-
-[[deps.Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
-git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec"
-uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.10.1"
-
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
@@ -584,24 +474,12 @@ version = "1.10.1"
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[[deps.ThreadedScans]]
-deps = ["ArgCheck"]
-git-tree-sha1 = "ca1ba3000289eacba571aaa4efcefb642e7a1de6"
-uuid = "24d252fe-5d94-4a69-83ea-56a14333d47a"
-version = "0.1.0"
-
 [[deps.TranscodingStreams]]
 deps = ["Random", "Test"]
 git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.9.11"
 
-[[deps.Transducers]]
-deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
-git-tree-sha1 = "c42fa452a60f022e9e087823b47e5a5f8adc53d5"
-uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
-version = "0.4.75"
-
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -614,12 +492,6 @@ deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
 version = "1.2.12+3"
 
-[[deps.ZygoteRules]]
-deps = ["MacroTools"]
-git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
-uuid = "700de1a5-db45-46bc-99cf-38207098b444"
-version = "0.2.2"
-
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

Project.toml

@@ -6,6 +6,7 @@ version = "0.4.0"
 [deps]
 Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -17,6 +18,7 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 COSMO = "0.8"
 Groups = "0.7"
 IntervalArithmetic = "0.20"
+IntervalMatrices = "0.8"
 JuMP = "1.3"
 ProgressMeter = "1.7"
 SCS = "1.1"

scripts/G₂_Adj.jl

@@ -1,6 +1,6 @@
 using LinearAlgebra
 BLAS.set_num_threads(8)
-
+using MKL_jll
 ENV["OMP_NUM_THREADS"] = 4
 
 using Groups
@@ -84,16 +84,59 @@ unit = Δ
 
 warm = nothing
 
-solve_in_loop(
+let status = JuMP.OPTIMIZE_NOT_CALLED, warm = warm, eps = 1e-9
+    certified, λ = false, 0.0
+    while status ≠ JuMP.OPTIMAL
+        @time status, warm = PropertyT.solve(
             model,
-    wd,
-    varP;
-    logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
-    optimizer = scs_optimizer(;
-        eps = 1e-10,
-        max_iters = 50_000,
+            scs_optimizer(;
+                linear_solver = SCS.MKLDirectSolver,
+                eps = eps,
+                max_iters = 100_000,
                 accel = 50,
                 alpha = 1.95,
             ),
-    data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
+            warm,
         )
+
+        @info "reconstructing the solution"
+        Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP], eps = eps
+            PropertyT.__droptol!.(Ps, 100eps)
+            Qs = real.(sqrt.(Ps))
+            PropertyT.__droptol!.(Qs, eps)
+
+            PropertyT.reconstruct(Qs, wd)
+        end
+
+        @info "certifying the solution"
+        @time certified, λ = PropertyT.certify_solution(
+            elt,
+            unit,
+            JuMP.objective_value(model),
+            Q;
+            halfradius = HALFRADIUS,
+            augmented = true,
+        )
+    end
+
+    if certified && λ > 0
+        Κ(λ, S) = round(sqrt(2λ / length(S)), Base.RoundDown; digits = 5)
+        @info "Certified result: $G has property (T):" N λ Κ(λ, S)
+    else
+        @info "Could NOT certify the result:" certified λ
+    end
+end
+
+# solve_in_loop(
+#     model,
+#     wd,
+#     varP;
+#     logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
+#     optimizer = scs_optimizer(;
+#         eps = 1e-10,
+#         max_iters = 50_000,
+#         accel = 50,
+#         alpha = 1.95,
+#     ),
+#     data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
+# )

src/PropertyT.jl

@@ -3,7 +3,6 @@ module PropertyT
 using LinearAlgebra
 using SparseArrays
 
-using IntervalArithmetic
 using JuMP
 
 using Groups

src/certify.jl

@@ -1,3 +1,6 @@
+import IntervalArithmetic
+import IntervalMatrices
+
 function augment_columns!(Q::AbstractMatrix)
     for c in eachcol(Q)
         c .-= sum(c) ./ length(c)
@@ -63,7 +66,7 @@ function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
     suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
 
     eq_sign = let T = eltype(residual)
-        if T <: Interval
+        if T <: IntervalArithmetic.Interval
             "∈"
         elseif T <: Union{Rational,Integer}
             "="
@@ -119,8 +122,10 @@ function certify_solution(
         return false, λ_flpoint
     end
 
-    λ_int = @interval(λ)
-    Q_int = [@interval(q) for q in Q]
+    λ_int = IntervalArithmetic.@interval(λ)
+    Q_int = IntervalMatrices.IntervalMatrix([
+        IntervalArithmetic.@interval(q) for q in Q
+    ])
 
     check, sos_int = @time if should_we_augment
         @info("Projecting columns of Q to the augmentation ideal...")
@@ -141,5 +146,5 @@ function certify_solution(
     λ_certified =
         sufficient_λ(elt, orderunit, λ_int, sos_int; halfradius = halfradius)
 
-    return check && inf(λ_certified) > 0.0, λ_certified
+    return check && IntervalArithmetic.inf(λ_certified) > 0.0, λ_certified
 end

src/roots.jl

@@ -7,73 +7,48 @@ export Root, isproportional, isorthogonal, ~, ⟂
 
 abstract type AbstractRoot{N,T} end
 
-struct Root{N,T} <: AbstractRoot{N,T}
-    coord::SVector{N,T}
-end
-
-Root(a) = Root(SVector(a...))
-
-function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
-    M == N || return false
-    r.coord == s.coord || return false
-    return true
-end
-
-Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
-
-Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord)
-Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord)
-Base.:-(r::Root{N}) where {N} = Root(-r.coord)
-
-Base.:*(a::Number, r::Root) = Root(a * r.coord)
-Base.:*(r::Root, a::Number) = a * r
-
-Base.length(r::AbstractRoot) = norm(r, 2)
-
-LinearAlgebra.norm(r::Root, p::Real = 2) = norm(r.coord, p)
-LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
+ℓ₂length(r::AbstractRoot) = norm(r, 2)
+ambient_dim(r::AbstractRoot) = length(r)
+Base.:*(r::AbstractRoot, a::Number) = a * r
 
 cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
 
-function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
-    N == M || return false
+function isproportional(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || return false
     val = abs(cos_angle(α, β))
     return isapprox(val, one(val); atol = eps(one(val)))
 end
 
-function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
-    N == M || return false
+function isorthogonal(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || return false
     val = cos_angle(α, β)
     return isapprox(val, zero(val); atol = eps(one(val)))
 end
 
-function _positive_direction(α::Root{N}) where {N}
-    v = α.coord + 1 / (N * 100) * rand(N)
-    return Root{N,Float64}(v / norm(v, 2))
-end
-
-function positive(roots::AbstractVector{<:Root{N}}) where {N}
+function positive(roots::AbstractVector{<:AbstractRoot})
+    isempty(roots) && return empty(roots)
     pd = _positive_direction(first(roots))
     return filter(α -> dot(α, pd) > 0.0, roots)
 end
 
-function Base.show(io::IO, r::Root)
+function Base.show(io::IO, r::AbstractRoot)
     return print(io, "Root $(r.coord)")
 end
 
-function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
-    lngth² = sum(x -> x^2, r.coord)
-    l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
+function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
+    l₂l = ℓ₂length(r)
+    l = isinteger(l₂l) ? "$(l₂l)" : "√$(l₂l^2)"
     return print(io, "Root in ℝ^$N of length $l\n", r.coord)
 end
 
-𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
-𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
-
-reflection(α::Root, β::Root) = β - Int(2dot(α, β) / dot(α, α)) * α
-function cartan(α, β)
+function reflection(α::AbstractRoot, β::AbstractRoot)
+    return β - Int(2dot(α, β) // dot(α, α)) * α
+end
+function cartan(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || throw("incompatible ambient dimensions")
     return [
-        length(reflection(a, b) - b) / length(a) for a in (α, β), b in (α, β)
+        ℓ₂length(reflection(a, b) - b) / ℓ₂length(a) for a in (α, β),
+        b in (α, β)
     ]
 end
 
@@ -124,7 +99,10 @@ function classify_root_system(
     end
 end
 
-function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
+function proportional_root_from_system(
+    Ω::AbstractVector{<:AbstractRoot},
+    α::AbstractRoot,
+)
     k = findfirst(v -> isproportional(α, v), Ω)
     if isnothing(k)
         error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
@@ -132,31 +110,31 @@ function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
     return Ω[k]
 end
 
-struct Plane{R<:Root}
+struct Plane{R<:AbstractRoot}
     v1::R
     v2::R
     vectors::Vector{R}
 end
 
-function Plane(α::Root, β::Root)
+function Plane(α::AbstractRoot, β::AbstractRoot)
     return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
 end
 
-function Base.in(r::Root, plane::Plane)
+function Base.in(r::AbstractRoot, plane::Plane)
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
-function _islong(α::Root, Ω)
-    lα = length(α)
-    return any(r -> lα - length(r) > eps(lα), Ω)
+function _islong(α::AbstractRoot, Ω)
+    lα = ℓ₂length(α)
+    return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
 end
 
 function classify_sub_root_system(
-    Ω::AbstractVector{<:Root{N}},
-    α::Root{N},
-    β::Root{N},
+    Ω::AbstractVector{<:AbstractRoot{N}},
+    α::AbstractRoot{N},
+    β::AbstractRoot{N},
 ) where {N}
-    @assert 1 ≤ length(unique(length, Ω)) ≤ 2
+    @assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
     v = proportional_root_from_system(Ω, α)
     w = proportional_root_from_system(Ω, β)
 
@@ -197,4 +175,45 @@ function classify_sub_root_system(
     throw("Unknown root system: $subsystem")
 end
 
+## concrete implementation:
+struct Root{N,T} <: AbstractRoot{N,T}
+    coord::SVector{N,T}
+end
+
+Root(a) = Root(SVector(a...))
+
+# convienience constructors
+𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
+𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
+
+function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
+    M == N || return false
+    r.coord == s.coord || return false
+    return true
+end
+
+Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
+
+function Base.:+(r::Root, s::Root)
+    ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
+    return Root(r.coord + s.coord)
+end
+
+function Base.:-(r::Root, s::Root)
+    ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
+    return Root(r.coord - s.coord)
+end
+Base.:-(r::Root) = Root(-r.coord)
+
+Base.:*(a::Number, r::Root) = Root(a * r.coord)
+
+Base.length(r::Root) = length(r.coord)
+
+LinearAlgebra.norm(r::Root, p::Real = 2) = norm(r.coord, p)
+LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
+
+function _positive_direction(α::Root{N}) where {N}
+    v = α.coord + 1 / (N * 100) * rand(N)
+    return Root{N,Float64}(v / norm(v, 2))
+end
 end # of module Roots

test/Chevalley.jl

@@ -22,7 +22,7 @@ end
 @testset "Exceptional root systems" begin
     @testset "F4" begin
         F4 = let Σ = PermutationGroups.PermGroup(perm"(1,2,3,4)", perm"(1,2)")
-            long = let x = (1.0, 1.0, 0.0, 0.0)
+            long = let x = (1, 1, 0, 0) .// 1
                 PropertyT.Roots.Root.(
                     union(
                         (x^g for g in Σ),
@@ -32,14 +32,14 @@ end
                 )
             end
 
-            short = let x = (1.0, 0.0, 0.0, 0.0)
+            short = let x = (1, 0, 0, 0) .// 1
                 PropertyT.Roots.Root.(
                     union((x^g for g in Σ), ((-1 .* x)^g for g in Σ))
                 )
             end
 
             signs = collect(Iterators.product(fill([-1, +1], 4)...))
-            halfs = let x = 1 / 2 .* (1.0, 1.0, 1.0, 1.0)
+            halfs = let x = (1, 1, 1, 1) .// 2
                 PropertyT.Roots.Root.(union(x .* sgn for sgn in signs))
             end
 
@@ -49,15 +49,15 @@ end
         @test length(F4) == 48
 
         a = F4[1]
-        @test isapprox(length(a), sqrt(2))
+        @test isapprox(PropertyT.Roots.ℓ₂length(a), sqrt(2))
         b = F4[6]
-        @test isapprox(length(b), sqrt(2))
+        @test isapprox(PropertyT.Roots.ℓ₂length(b), sqrt(2))
         c = a + b
-        @test isapprox(length(c), 2.0)
+        @test isapprox(PropertyT.Roots.ℓ₂length(c), 2.0)
         @test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
 
-        long = F4[findfirst(r -> length(r) == sqrt(2), F4)]
-        short = F4[findfirst(r -> length(r) == 1.0, F4)]
+        long = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == sqrt(2), F4)]
+        short = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == 1.0, F4)]
 
         subtypes = Set([:C₂, :A₂, Symbol("A₁×C₁")])
 
@@ -94,7 +94,7 @@ end
                     perm"(1,2,3,4,5,6,7,8)",
                     perm"(1,2)",
                 )
-                long = let x = (1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+                long = let x = (1, 1, 0, 0, 0, 0, 0, 0) .// 1
                     PropertyT.Roots.Root.(
                         union(
                             (x^g for g in Σ),
@@ -108,7 +108,7 @@ end
                     p for p in Iterators.product(fill([-1, +1], 8)...) if
                     iseven(count(==(-1), p))
                 )
-                halfs = let x = 1 / 2 .* ntuple(i -> 1.0, 8)
+                halfs = let x = (1, 1, 1, 1, 1, 1, 1, 1) .// 2
                     rts = unique(PropertyT.Roots.Root(x .* sgn) for sgn in signs)
                 end
 
@@ -119,7 +119,7 @@ end
 
         @testset "E8" begin
             @test length(E8) == 240
-            @test all(r -> length(r) ≈ sqrt(2), E8)
+            @test all(r -> PropertyT.Roots.ℓ₂length(r) ≈ sqrt(2), E8)
 
             let Ω = E8, α = first(Ω)
                 counts = countmap([

test/optimizers.jl

@@ -44,6 +44,8 @@ function cosmo_optimizer(;
         "decompose" => decompose,
         "eps_abs" => eps,
         "eps_rel" => eps,
+        "eps_prim_inf" => eps,
+        "eps_dual_inf" => eps,
         "max_iter" => max_iters,
         "verbose" => verbose,
         "verbose_timing" => verbose_timing,