From c0e73abd06a0a504354dec0047b1af86dd4e4f0c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Tue, 20 Aug 2024 18:30:51 +0200
Subject: [PATCH] Fixes for BigFloat

---
 src/extract.jl    |  4 ++++
 test/Project.toml |  1 +
 test/extract.jl   | 45 +++++++++++++++++++++++++--------------------
 3 files changed, 30 insertions(+), 20 deletions(-)

diff --git a/src/extract.jl b/src/extract.jl
index 8816848..2122e6b 100644
--- a/src/extract.jl
+++ b/src/extract.jl
@@ -30,6 +30,10 @@ function solve_weight(
     centers,
     ::MomentMatrixWeightSolver,
 ) where {T}
+    if isempty(centers)
+        # Failing for `BigFloat`
+        return T[]
+    end
     vars = MP.variables(ν)
     A = Matrix{T}(undef, length(ν.Q.Q), length(centers))
     vbasis = vectorized_basis(ν)
diff --git a/test/Project.toml b/test/Project.toml
index f1a5512..c839dd0 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e"
diff --git a/test/extract.jl b/test/extract.jl
index 2a5a241..30cd030 100644
--- a/test/extract.jl
+++ b/test/extract.jl
@@ -3,6 +3,7 @@ using MultivariatePolynomials
 import MultivariateBases as MB
 using SemialgebraicSets
 using MultivariateMoments
+import GenericLinearAlgebra # Needed for `svd` to work on `Matrix{BigFloat}`
 
 struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end
 SemialgebraicSets.promote_for(::Type{T}, ::Type{DummySolver}) where {T} = T
@@ -34,18 +35,18 @@ function _atoms(atoms, rank_check, solver)
     atoms = atomic_measure(ν, rank_check, solver)
     @test atoms !== nothing
     if !isnothing(atoms)
-        @test atoms ≈ η
+        @test atoms ≈ η rtol = 1e-8
     end
 end
 
-function atoms_1(rank_check, solver)
-    atoms = [[1, 2]]
+function atoms_1(T, rank_check, solver)
+    atoms = [T[1, 2]]
     _atoms(atoms, rank_check, solver)
     return
 end
 
-function atoms_2(rank_check, solver)
-    atoms = [[1, 1], [1, -1]]
+function atoms_2(T, rank_check, solver)
+    atoms = [T[1, 1], T[1, -1]]
     _atoms(atoms, rank_check, solver)
     return
 end
@@ -351,10 +352,10 @@ function jcg14_6_1(rank_check, ok::Bool = true)
     end
 end
 
-function large_norm(rank_check)
+function large_norm(T, rank_check)
     # If the norm of `M` is given to `rref!` instead of `√||M||`, `atomic_measure` will error.
     Mod.@polyvar x[1:2]
-    Q = [
+    Q = T[
         586.8034549325414 -800.152792847183
         -800.152792847183 2749.376669556701
     ]
@@ -390,19 +391,21 @@ _short(x) = _short(string(x))
 
 function test_extract()
     default_solver = SemialgebraicSets.default_algebraic_solver([1.0x - 1.0x])
-    @testset "$(_short(solver))" for solver in [
-        FlatExtension(),
-        FlatExtension(NewtonTypeDiagonalization()),
-        Echelon(),
-        ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), Echelon()),
-        ShiftNullspace(),
-        ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), ShiftNullspace()),
-    ]
+    @testset "$T $(_short(solver))" for T in [Float64, BigFloat],
+        solver in [
+            FlatExtension(),
+            FlatExtension(NewtonTypeDiagonalization()),
+            Echelon(),
+            ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), Echelon()),
+            ShiftNullspace(),
+            ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), ShiftNullspace()),
+        ]
+
         @testset "Atom 1" begin
-            atoms_1(1e-10, solver)
+            atoms_1(T, 1e-10, solver)
         end
         @testset "Atom 2" begin
-            atoms_2(1e-10, solver)
+            atoms_2(T, 1e-10, solver)
         end
     end
     @testset "hl05_2_3 $(_short(lrc))" for lrc in
@@ -465,11 +468,13 @@ function test_extract()
     lpj20_3_9(FixedRank(7), 0)
     jcg14_6_1(6e-3)
     jcg14_6_1(8e-4, false)
-    large_norm(1e-2)
+    @testset "$T" for T in [Float64, BigFloat]
+        large_norm(T, 1e-2)
+    end
     @testset "No comm fix $(_short(solver)) $T" for solver in
                                                     [ShiftNullspace()],
-        T in [Float64]
-        #, BigFloat]
+        T in [Float64, BigFloat]
+
         no_com(solver, T)
     end
     return