different pivoting for hmatrices

Project.toml

@@ -22,6 +22,7 @@ StaticArrays = "1.9.13"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -30,4 +31,4 @@ TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 
 [targets]
-test = ["Aqua", "JET", "JuliaFormatter", "Test", "TestItemRunner", "TestItems", "Random"]
+test = ["Aqua", "JET", "JuliaFormatter", "Test", "TestItemRunner", "TestItems", "Random", "HMatrices"]

ext/ACAHMatrices/hmatrices_interface.jl

@@ -19,14 +19,8 @@ function (aca::ACA)(
         bufs.B.data[1:(length(irange) * maxrank)], (length(irange), maxrank)
     )
 
-    cur_aca = AdaptiveCrossApproximation.ACA(;
-        rowpivoting=AdaptiveCrossApproximation.MaximumValue(zeros(Bool, length(irange))),
-        columnpivoting=AdaptiveCrossApproximation.MaximumValue(zeros(Bool, length(jrange))),
-        convergence=aca.convergence,
-    )
-    npivots, U, V = cur_aca(
-        K, rowBuffer, colBuffer, maxrank; rowidcs=irange, colidcs=jrange
-    )
+    aca = aca(K, irange, jrange)
+    npivots, U, V = aca(K, rowBuffer, colBuffer, maxrank; rowidcs=irange, colidcs=jrange)
     return HMatrices.RkMatrix(
         colBuffer[:, 1:npivots], Matrix(transpose(rowBuffer[1:npivots, :]))
     )

src/aca.jl

@@ -20,13 +20,55 @@ function ACA(;
     return ACA(rowpivoting, columnpivoting, convergence)
 end
 
+function (aca::ACA)(
+    K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int}
+)
+    return ACA(aca.rowpivoting(rowidcs), aca.columnpivoting(colidcs), aca.convergence())
+end
+
+function (aca::ACA{PS,PS,CC})(
+    K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int}
+) where {PS<:PivStrat,CC<:RandomSampling}
+    convergence = aca.convergence(K, rowidcs, colidcs)
+    rowpivoting = if aca.rowpivoting isa RandomSamplingPivoting
+        aca.rowpivoting(convergence)
+    else
+        aca.rowpivoting(rowidcs)
+    end
+    columnpivoting = if aca.columnpivoting isa RandomSamplingPivoting
+        aca.columnpivoting(convergence)
+    else
+        aca.columnpivoting(colidcs)
+    end
+
+    return ACA(rowpivoting, columnpivoting, convergence)
+end
+
+function (aca::ACA{RP,CP,CC})(
+    K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int}
+) where {RP,CP,CC<:CombinedConvCrit}
+    convergence = aca.convergence(K, rowidcs, colidcs)
+    rowpivoting = if aca.rowpivoting isa CombinedPivStrat
+        aca.rowpivoting(convergence, rowidcs)
+    else
+        aca.rowpivoting(rowidcs)
+    end
+    columnpivoting = if aca.columnpivoting isa CombinedPivStrat
+        aca.columnpivoting(convergence, colidcs)
+    else
+        aca.columnpivoting(colidcs)
+    end
+
+    return ACA(rowpivoting, columnpivoting, convergence)
+end
+
 function (aca::ACA)(
     A,
     rowbuffer::AbstractMatrix{K},
     colbuffer::AbstractMatrix{K},
     maxrank::Int;
-    rowidcs=Vector(1:size(colbuffer, 1)),
-    colidcs=Vector(1:size(rowbuffer, 2)),
+    rowidcs = Vector(1:size(colbuffer, 1)),
+    colidcs = Vector(1:size(rowbuffer, 2)),
 ) where {K}
     rows = Int[1]
     cols = Int[]

src/convergence/combinedconvcrit.jl

@@ -19,3 +19,17 @@ function (convcrit::CombinedConvCrit)(
     end
     return npivot, any(convcrit.isconverged)
 end
+
+function (convcrit::CombinedConvCrit)(
+    K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int}
+)
+    curr_crits = Vector{ConvCrit}(undef, length(convcrit.crits))
+    for (i, crit) in enumerate(convcrit.crits)
+        if isa(crit, RandomSampling)
+            curr_crits[i] = crit(K, rowidcs, colidcs)
+        else
+            curr_crits[i] = crit()
+        end
+    end
+    return CombinedConvCrit(curr_crits, zeros(Bool, length(curr_crits)))
+end

src/convergence/estimation.jl

@@ -7,7 +7,10 @@ function FNormEstimator(tolerance::F) where {F}
     return FNormEstimator(F(0.0), tolerance)
 end
 
-(::FNormEstimator{F})() where {F} = FNormEstimator(0.0, 1e-4)
+function (cc::FNormEstimator{F})() where {F}
+    return FNormEstimator(0.0, cc.tol)
+end
+
 tolerance(cc::FNormEstimator) = cc.tol
 
 function (convcrit::FNormEstimator{F})(

src/convergence/extrapolation.jl

@@ -7,7 +7,9 @@ function FNormExtrapolator(tol::F) where {F}
     return FNormExtrapolator(F[], FNormEstimator(F(0.0), tol))
 end
 
-(::FNormExtrapolator{F})() where {F} = FNormExtrapolator(F[], FNormEstimator(0.0, 1e-4))
+function (cc::FNormExtrapolator{F})() where {F}
+    return FNormExtrapolator(F[], FNormEstimator(0.0, cc.estimator.tol))
+end
 tolerance(cc::FNormExtrapolator) = cc.estimator.tol
 
 # ACA

src/convergence/randomsampling.jl

@@ -12,9 +12,17 @@ end
 function RandomSampling(
     ::Type{K}; factor::F=1.0, nsamples::Int=0, tol::F=1e-4
 ) where {K,F<:Real}
-    return RandomSampling{F,K}(
-        F(0.0), nsamples, factor, zeros(Int, 0, 0), zeros(K, 0, 0), tol
-    )
+    return RandomSampling{F,K}(F(0.0), nsamples, factor, zeros(Int, 0, 0), zeros(K, 0), tol)
+end
+function (cc::RandomSampling)(
+    K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int}
+)
+    rowlen = length(rowidcs)
+    collen = length(colidcs)
+    nsamples = cc.nsamples == 0 ? cc.factor * (rowlen + collen) : cc.nsamples
+    indices = hcat(rand(1:rowlen, nsamples), rand(1:collen, nsamples))
+    rest = [K[rc[1], rc[2]][1] for rc in eachrow(indices)]
+    return RandomSampling(0.0, nsamples, cc.factor, indices, rest, cc.tol)
 end
 
 tolerance(cc::RandomSampling) = cc.tol

src/pivoting/combinedpivstrat.jl

@@ -9,7 +9,6 @@ function (pivstrat::CombinedPivStrat)()
 end
 
 function (pivstrat::CombinedPivStrat)(rc::AbstractArray)
-    #öprintln(pivstrat.convcrit.isconverged)
     length(pivstrat.strats) > length(pivstrat.convcrit.isconverged) &&
         push!(pivstrat.convcrit.isconverged, false)
     for (i, conv) in enumerate(pivstrat.convcrit.isconverged)
@@ -18,3 +17,17 @@ function (pivstrat::CombinedPivStrat)(rc::AbstractArray)
         return pivstrat.strats[i](rc)
     end
 end
+
+function (pivstrat::CombinedPivStrat)(
+    convergence::CombinedConvCrit, idcs::AbstractArray{Int}
+)
+    curr_strats = Vector{PivStrat}(undef, length(pivstrat.strats))
+    for (i, strat) in enumerate(pivstrat.strats)
+        if isa(strat, RandomSamplingPivoting)
+            curr_strats[i] = strat(convergence.crits[i])
+        else
+            curr_strats[i] = strat(idcs)
+        end
+    end
+    return CombinedPivStrat(convergence, curr_strats)
+end

src/pivoting/filldistance.jl

@@ -8,7 +8,7 @@ function FillDistance(pos::Vector{SVector{D,F}}) where {D,F<:Real}
     return FillDistance{D,F}(F[], pos)
 end
 
-function (pivstrat::FillDistance{D,F})(idcs::Vector{Int}) where {D,F}
+function (pivstrat::FillDistance{D,F})(idcs::AbstractArray{Int}) where {D,F}
     return FillDistance{D,F}(zeros(F, length(idcs)), pivstrat.pos[idcs])
 end
 

src/pivoting/lejapoints.jl

@@ -16,7 +16,7 @@ function Leja2(pos::Vector{SVector{D,F}}) where {D,F}
     return Leja2{D,F}(F[], pos)
 end
 
-function (pivstrat::Leja2{D,F})(idcs::Vector{Int}) where {D,F}
+function (pivstrat::Leja2{D,F})(idcs::AbstractArray{Int}) where {D,F}
     return Leja2{D,F}(zeros(F, length(idcs)), pivstrat.pos[idcs])
 end
 

src/pivoting/maxvalue.jl

@@ -6,6 +6,7 @@ end
 MaximumValue() = MaximumValue(Bool[])
 
 (::MaximumValue)(len::Int) = MaximumValue(zeros(Bool, len))
+(::MaximumValue)(idcs::AbstractArray{Int}) = MaximumValue(zeros(Bool, length(idcs)))
 
 function (pivstrat::MaximumValue)()
     pivstrat.usedidcs[1] = true

src/pivoting/randomsampling.jl

@@ -6,3 +6,7 @@ end
 function (pivstrat::RandomSamplingPivoting{F,K})(::AbstractArray) where {F<:Real,K}
     return pivstrat.convcrit.indices[argmax(pivstrat.convcrit.rest), pivstrat.rc]
 end
+
+function (pivstrat::RandomSamplingPivoting{F,K})(convergence::ConvCrit) where {F<:Real,K}
+    return RandomSamplingPivoting(convergence, pivstrat.rc)
+end

test/runtests.jl

@@ -5,6 +5,7 @@ using AdaptiveCrossApproximation
     include("test_aca.jl")
     include("test_convergence.jl")
     include("test_pivoting.jl")
+    include("test_ACAHMatrices.jl")
 end
 
 @testitem "Code quality (Aqua.jl)" begin

test/test_ACAHMatrices.jl

@@ -1,30 +1,87 @@
-using BEAST
 using HMatrices
-using CompScienceMeshes
+using StaticArrays
 using AdaptiveCrossApproximation
+using Random
 using Test
 using LinearAlgebra
 
-Γ = meshsphere(1.0, 0.2)
-op = Helmholtz3D.singlelayer()
-spaceX = lagrangecxd0(Γ)
-dim = length(spaceX.pos)
-xclt = ClusterTree(spaceX.pos)
+const Point3D = SVector{3,Float64}
+X = rand(Point3D, 1000)
+Y = rand(Point3D, 1000)
+const k = 2π
+function G(buf, irange, jrange)
+    for (a, i) in enumerate(irange), (b, j) in enumerate(jrange)
+        buf[a, b] = exp(-k * norm(X[i] - Y[j]))
+    end
+    return buf
+end
+
+function G(x, y)
+    return exp(-k * norm(x - y))
+end
 
-K = HMatrices.KernelMatrix(op, spaceX, spaceX);
+K = KernelMatrix(G, X, Y)
+dim = size(K, 2)
 
-rtols = [10.0^i for i in collect(-4:-1:-10)]
+rtols = [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14]
 tst_vec = rand(dim)
 
-fullmat = BEAST.assemble(op, spaceX, spaceX)
-trueResult = fullmat * tst_vec;
-##
+trueResult = K * tst_vec;
 
+# PivStrat :: MaximumValue
 for rtol in rtols
-    K = HMatrices.KernelMatrix(op, spaceX, spaceX)
     aca = ACA(; convergence=FNormEstimator(0.0, rtol))
     hmat = HMatrices.assemble_hmatrix(K; comp=aca)
-    @test size(hmat, 1) == size(fullmat, 1)
-    @test size(hmat, 2) == size(fullmat, 2)
+    @test size(hmat, 1) == size(K, 1)
+    @test size(hmat, 2) == size(K, 2)
+    @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol
+end
+# PivStrat :: leja2
+for rtol in rtols
+    aca = ACA(; rowpivoting=Leja2(X), convergence=FNormEstimator(0.0, rtol))
+    hmat = HMatrices.assemble_hmatrix(K; comp=aca)
+    @test size(hmat, 1) == size(K, 1)
+    @test size(hmat, 2) == size(K, 2)
+    @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol
+
+    aca = ACA(; columnpivoting=Leja2(Y), convergence=FNormEstimator(0.0, rtol))
+    hmat = HMatrices.assemble_hmatrix(K; comp=aca)
+    @test size(hmat, 1) == size(K, 1)
+    @test size(hmat, 2) == size(K, 2)
+    @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol
+end
+
+# PivStrat :: FillDistance
+for rtol in rtols
+    aca = ACA(; rowpivoting=FillDistance(X), convergence=FNormEstimator(0.0, rtol))
+    hmat = HMatrices.assemble_hmatrix(K; comp=aca)
+    @test size(hmat, 1) == size(K, 1)
+    @test size(hmat, 2) == size(K, 2)
+    @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol
+
+    aca = ACA(; columnpivoting=FillDistance(Y), convergence=FNormEstimator(0.0, rtol))
+    hmat = HMatrices.assemble_hmatrix(K; comp=aca)
+    @test size(hmat, 1) == size(K, 1)
+    @test size(hmat, 2) == size(K, 2)
     @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol
 end
+
+# PivStrat :: CombinedPivStrat
+for rtol in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10]#, 1e-12, 1e-14]
+    Random.seed!(1)
+    cc1 = AdaptiveCrossApproximation.FNormEstimator(rtol)
+    cc2 = AdaptiveCrossApproximation.RandomSampling(
+        Float64; nsamples=100, factor=1.0, tol=rtol
+    )
+    convergence = AdaptiveCrossApproximation.CombinedConvCrit([cc1, cc2], zeros(Bool, 2))
+
+    ps1 = MaximumValue()
+    ps2 = AdaptiveCrossApproximation.RandomSamplingPivoting(cc2, 1)
+    rp = AdaptiveCrossApproximation.CombinedPivStrat(convergence, [ps1, ps2])
+    aca = ACA(; rowpivoting=rp, convergence=convergence)
+    hmat = HMatrices.assemble_hmatrix(K; comp=aca)
+    @test size(hmat, 1) == size(K, 1)
+    @test size(hmat, 2) == size(K, 2)
+    @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol
+end
+