fixes and keep_below

Author: lkdvos

src/implementations/orthnull.jl

@@ -203,11 +203,11 @@ end
 # --------------------------------
 function null_truncation_strategy(; atol=nothing, rtol=nothing, maxnullity=nothing)
     if isnothing(maxnullity) && isnothing(atol) && isnothing(rtol)
-        return NoTruncation()
+        return notrunc()
     end
     atol = @something atol 0
     rtol = @something rtol 0
-    trunc = trunctol(; atol, rtol, rev=true)
+    trunc = trunctol(; atol, rtol, keep_below=true)
     return !isnothing(maxnullity) ? trunc & truncrank(maxnullity; rev=false) : trunc
 end
 

src/implementations/truncation.jl

@@ -41,29 +41,37 @@ function findtruncated(values::AbstractVector, strategy::TruncationByOrder)
     return partialsortperm(values, 1:howmany; strategy.by, strategy.rev)
 end
 function findtruncated_sorted(values::AbstractVector, strategy::TruncationByOrder)
-    howmany = min(strategy.howmany, length(values))
-    return strategy.rev ? (1:howmany) : ((length(values) - howmany + 1):length(values))
+    @assert strategy.by === abs
+    if strategy.by === abs
+        howmany = min(strategy.howmany, length(values))
+        return strategy.rev ? (1:howmany) : ((length(values) - howmany + 1):length(values))
+    else
+        return findtruncated(values, strategy)
+    end
 end
 
 function findtruncated(values::AbstractVector, strategy::TruncationByFilter)
-    ind = findall(strategy.filter, values)
-    return ind
+    # pre-allocate bitvector to enforce the filter function returns a Bool
+    mask = similar(BitArray, eachindex(values))
+    mask .= strategy.filter.(values)
+    return mask
 end
 
 function findtruncated(values::AbstractVector, strategy::TruncationByValue)
     atol = max(strategy.atol, strategy.rtol * norm(values, strategy.p))
-    filter = (strategy.rev ? ≤(atol) : ≥(atol)) ∘ strategy.by
-    return findall(filter, values)
+    filter = (strategy.keep_below ? ≤(atol) : ≥(atol)) ∘ strategy.by
+    return findtruncated(values, truncfilter(filter))
 end
-function findtruncated_sorted(values::AbstractVector, strategy::TruncationByValue)
+function findtruncated_svd(values::AbstractVector, strategy::TruncationByValue)
+    strategy.by === abs || return findtruncated(values, strategy)
+
     atol = max(strategy.atol, strategy.rtol * norm(values, strategy.p))
-    @assert strategy.by === abs || strategy.by === real "sorting strategy incompatible with implementation"
-    if strategy.rev
-        i = searchsortedfirst(values, atol; by=strategy.by, rev=true)
-        return i:length(values)
-    else
-        i = searchsortedlast(values, atol; by=strategy.by, rev=true)
+    if strategy.keep_above
+        i = searchsortedlast(values, atol; by=abs, rev=true)
         return 1:i
+    else
+        i = searchsortedfirst(values, atol; by=abs, rev=true)
+        return i:length(values)
     end
 end
 
@@ -97,10 +105,20 @@ function _truncerr_impl(values::AbstractVector, I; atol::Real=0, rtol::Real=0, p
 end
 
 function findtruncated(values::AbstractVector, strategy::TruncationIntersection)
-    inds = map(Base.Fix1(findtruncated, values), strategy.components)
-    return intersect(inds...)
+    return mapreduce(Base.Fix1(findtruncated, values), _ind_intersect, strategy.components;
+                     init=trues(length(values)))
 end
 function findtruncated_sorted(values::AbstractVector, strategy::TruncationIntersection)
-    inds = map(Base.Fix1(findtruncated_sorted, values), strategy.components)
-    return intersect(inds...)
+    return mapreduce(Base.Fix1(findtruncated_sorted, values), _ind_intersect,
+                     strategy.components; init=trues(length(values)))
 end
+
+# when one of the ind selections is a bitvector, have to handle differently
+function _ind_intersect(A::AbstractVector{Bool}, B::AbstractVector)
+    result = falses(length(A))
+    result[B] .= @view A[B]
+    return result
+end
+_ind_intersect(A::AbstractVector, B::AbstractVector{Bool}) = _ind_intersect(B, A)
+_ind_intersect(A::AbstractVector{Bool}, B::AbstractVector{Bool}) = A .& B
+_ind_intersect(A, B) = intersect(A, B)

src/interface/truncation.jl

@@ -41,6 +41,8 @@ Truncation strategy that does nothing, and keeps all the values.
 """
 notrunc() = NoTruncation()
 
+# TODO: Base.Ordering
+
 """
     TruncationByOrder(howmany::Int, by::Function, rev::Bool)
 
@@ -59,7 +61,10 @@ end
 
 Truncation strategy to keep the first `howmany` values when sorted according to `by` or the last `howmany` if `rev` is true.
 """
-truncrank(howmany::Integer; by=abs, rev::Bool=true) = TruncationByOrder(howmany, by, rev)
+function truncrank(howmany::Integer; by=abs, rev::Bool=true)
+    order = ...
+    return TruncationByOrder(howmany, order)
+end
 
 """
     TruncationByFilter(filter::Function)
@@ -80,31 +85,31 @@ Truncation strategy to keep the values for which `filter` returns true.
 truncfilter(f) = TruncationByFilter(f)
 
 """
-    TruncationByValue(atol::Real, rtol::Real, p::Real, by, rev::Bool=false)
+    TruncationByValue(atol::Real, rtol::Real, p::Real, by, keep_below::Bool=false)
 
-Truncation strategy to keep the values that satisfy `by(val) > max(atol, rtol * norm(values, p)`
-if `rev = false`, or discard them when `rev = true`.
+Truncation strategy to keep the values that satisfy `by(val) > max(atol, rtol * norm(values, p)`.
+If `keep_below = true`, discard these values instead.
 See also [`trunctol`](@ref)
 """
 struct TruncationByValue{T<:Real,P<:Real,F} <: TruncationStrategy
     atol::T
     rtol::T
     p::P
     by::F
-    rev::Bool
+    keep_below::Bool
 end
-function TruncationByValue(atol::Real, rtol::Real, p::Real=2, by=abs, rev::Bool=true)
-    return TruncationByValue(promote(atol, rtol)..., p, by, rev)
+function TruncationByValue(atol::Real, rtol::Real, p::Real=2, by=abs, keep_below::Bool=true)
+    return TruncationByValue(promote(atol, rtol)..., p, by, keep_below)
 end
 
 """
-    trunctol(; atol::Real=0, rtol::Real=0, p::Real=2, by=abs, rev::Bool=false)
+    trunctol(; atol::Real=0, rtol::Real=0, p::Real=2, by=abs, keep_below::Bool=false)
 
-Truncation strategy to keep the values that satisfy `by(val) > max(atol, rtol * norm(values, p)`
-if `rev = false`, or discard them when `rev = true`.
+Truncation strategy to keep the values that satisfy `by(val) > max(atol, rtol * norm(values, p)`.
+If `keep_below = true`, discard these values instead.
 """
-function trunctol(; atol::Real=0, rtol::Real=0, p::Real=2, by=abs, rev::Bool=false)
-    return TruncationByValue(atol, rtol, p, by, rev)
+function trunctol(; atol::Real=0, rtol::Real=0, p::Real=2, by=abs, keep_below::Bool=false)
+    return TruncationByValue(atol, rtol, p, by, keep_below)
 end
 
 """

test/algorithms.jl

@@ -50,7 +50,7 @@ end
                                  notrunc())
     end
 
-    alg = TruncatedAlgorithm(LAPACK_Simple(), trunctol(; atol=0.1, rev=true))
+    alg = TruncatedAlgorithm(LAPACK_Simple(), trunctol(; atol=0.1, keep_below=true))
     for f in (eig_trunc!, eigh_trunc!, svd_trunc!)
         @test @constinferred(select_algorithm(eig_trunc!, A, alg)) === alg
         @test_throws ArgumentError select_algorithm(eig_trunc!, A, alg; trunc=(; maxrank=2))

test/orthnull.jl

@@ -125,7 +125,7 @@ end
 
         rtol = eps(real(T))
         for (trunc_orth, trunc_null) in (((; rtol=rtol), (; rtol=rtol)),
-                                         (trunctol(; rtol), trunctol(; rtol, rev=true)))
+                                         (trunctol(; rtol), trunctol(; rtol, keep_below=true)))
             V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc=trunc_orth)
             N2 = @constinferred left_null!(copy!(Ac, A), N; trunc=trunc_null)
             @test V2 !== V

test/truncate.jl

@@ -53,10 +53,10 @@ using MatrixAlgebraKit: NoTruncation, TruncationIntersection, TruncationByOrder,
     @test @constinferred(findtruncated(values, strategy)) == [2, 3, 4, 5]
 
     for strategy in
-        (trunctol(; atol=0.4, rev=true), trunctol(; atol=0.2, by=identity, rev=true))
+        (trunctol(; atol=0.4, keep_below=true), trunctol(; atol=0.2, by=identity, keep_below=true))
         @test @constinferred(findtruncated(values, strategy)) == [1, 4]
     end
-    strategy = trunctol(; atol=0.2, rev=true)
+    strategy = trunctol(; atol=0.2, keep_below=true)
     @test @constinferred(findtruncated(values, strategy)) == [1]
 
     strategy = truncfilter(x -> 0.1 < x < 1)