Generalize blockisequal to arrays, introduce blockequals and blockisapprox (#477)

Generalize `blockisequal` to arrays, where the definition is that the
axes have the same block structure and the array values compare equal
with `isequal`. This PR also introduces `blockequals` and
`blockisapprox`, which are like `blockisequal` but use `==` and
`isapprox` respectively to compare values.

Co-authored-by: Sheehan Olver <[email protected]>

Author: mtfishman

docs/src/lib/public.md

@@ -32,6 +32,8 @@ Block
 BlockIndex
 blockaxes
 blockisequal
+blockequals
+blockisapprox
 blocksize
 blockfirsts
 blocklasts

src/BlockArrays.jl

@@ -5,7 +5,8 @@ using LinearAlgebra, ArrayLayouts, FillArrays
 export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
 export Block, getblock, getblock!, setblock!, eachblock, blocks
 export blockaxes, blocksize, blocklength, blockcheckbounds, BlockBoundsError, BlockIndex, BlockIndexRange
-export blocksizes, blocklengths, eachblockaxes, blocklasts, blockfirsts, blockisequal
+export blocksizes, blocklengths, blocklasts, blockfirsts, blockisequal, blockequals, blockisapprox
+export eachblockaxes
 export BlockRange, blockedrange, BlockedUnitRange, BlockedOneTo
 
 export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat, mortar

src/blockaxis.jl

@@ -225,9 +225,12 @@ end
 length(a::AbstractBlockedUnitRange) = isempty(blocklasts(a)) ? zero(eltype(a)) : Integer(last(blocklasts(a))-first(a)+oneunit(eltype(a)))
 
 """
+    blockisequal(a::AbstractArray, b::AbstractArray)
     blockisequal(a::AbstractUnitRange{<:Integer}, b::AbstractUnitRange{<:Integer})
 
-Check if `a` and `b` have the same block structure.
+Check if `a` and `b` have the same values (compared with `isequal`) and block structure.
+
+See also blockequals and blockisapprox.
 
 # Examples
 ```jldoctest
@@ -251,6 +254,29 @@ true
 
 julia> blockisequal(b1, b2)
 false
+
+julia> A = reshape([1:6;], 2, 3)
+2×3 Matrix{Int64}:
+ 1  3  5
+ 2  4  6
+
+julia> B1 = BlockedMatrix(A, [1,1], [1,2])
+2×2-blocked 2×3 BlockedMatrix{Int64}:
+ 1  │  3  5
+ ───┼──────
+ 2  │  4  6
+
+julia> B2 = BlockedMatrix(A, [1,1], [2,1])
+2×2-blocked 2×3 BlockedMatrix{Int64}:
+ 1  3  │  5
+ ──────┼───
+ 2  4  │  6
+
+julia> blockisequal(B1, B1)
+true
+
+julia> blockisequal(B1, B2)
+false
 ```
 """
 blockisequal(a::AbstractUnitRange{<:Integer}, b::AbstractUnitRange{<:Integer}) = first(a) == first(b) && blocklasts(a) == blocklasts(b)
@@ -265,6 +291,92 @@ blockisequal(::Tuple{}, ::Tuple{}) = true
 blockisequal(::Tuple, ::Tuple{}) = false
 blockisequal(::Tuple{}, ::Tuple) = false
 
+blockisequal(a::AbstractArray, b::AbstractArray) =
+    blockisequal(axes(a), axes(b)) && isequal(a, b)
+
+"""
+    blockequals(a::AbstractArray, b::AbstractArray)
+
+Check if `a` and `b` have the same values (compared with `==`) and block structure.
+
+See also blockisequal and blockisapprox.
+
+# Examples
+```jldoctest
+julia> A = reshape([1:6;], 2, 3)
+2×3 Matrix{Int64}:
+ 1  3  5
+ 2  4  6
+
+julia> B1 = BlockedMatrix(A, [1,1], [1,2])
+2×2-blocked 2×3 BlockedMatrix{Int64}:
+ 1  │  3  5
+ ───┼──────
+ 2  │  4  6
+
+julia> B2 = BlockedMatrix(A, [1,1], [2,1])
+2×2-blocked 2×3 BlockedMatrix{Int64}:
+ 1  3  │  5
+ ──────┼───
+ 2  4  │  6
+
+julia> blockequals(B1, B1)
+true
+
+julia> blockequals(B1, B2)
+false
+```
+"""
+blockequals(a::AbstractArray, b::AbstractArray) =
+    blockisequal(axes(a), axes(b)) && (a == b)
+
+"""
+    blockisapprox(a::AbstractArray, b::AbstractArray; kwargs...)
+
+Check if `a` and `b` have the same block structure and approximately the same
+values, compared with `isapprox`. Accepts the same keyword arguments as `isapprox`.
+
+See also blockisequal and blockequals.
+
+# Examples
+```jldoctest
+julia> A1 = reshape([1:6;], 2, 3)
+2×3 Matrix{Int64}:
+ 1  3  5
+ 2  4  6
+
+julia> A2 = A1 .+ 1e-5
+2×3 Matrix{Float64}:
+ 1.00001  3.00001  5.00001
+ 2.00001  4.00001  6.00001
+
+julia> B1 = BlockedMatrix(A1, [1,1], [1,2])
+2×2-blocked 2×3 BlockedMatrix{Int64}:
+ 1  │  3  5
+ ───┼──────
+ 2  │  4  6
+
+julia> B2 = BlockedMatrix(A2, [1,1], [1,2])
+2×2-blocked 2×3 BlockedMatrix{Float64}:
+ 1.00001  │  3.00001  5.00001
+ ─────────┼──────────────────
+ 2.00001  │  4.00001  6.00001
+
+julia> B3 = BlockedMatrix(A2, [1,1], [2,1])
+2×2-blocked 2×3 BlockedMatrix{Float64}:
+ 1.00001  3.00001  │  5.00001
+ ──────────────────┼─────────
+ 2.00001  4.00001  │  6.00001
+
+julia> blockisapprox(B1, B2; atol=1e-4)
+true
+
+julia> blockisapprox(B1, B3; atol=1e-4)
+false
+```
+"""
+blockisapprox(a::AbstractArray, b::AbstractArray; kwargs...) =
+    blockisequal(axes(a), axes(b)) && isapprox(a, b; kwargs...)
 
 _shift_blocklengths(::AbstractBlockedUnitRange, bl, f) = bl
 _shift_blocklengths(::Any, bl, f) = bl .+ (f - 1)

test/test_blockindices.jl

@@ -937,15 +937,49 @@ end
     @test B[1,2] == 0
 end
 
-@testset "blockisequal" begin
-    B = BlockArray(rand(4,4), [1,3], [1,3])
+@testset "blockisequal, blockequals" begin
+    A = rand(4,4)
+    B = BlockArray(A, [1,3], [1,3])
+    B2 = BlockArray(A, [1,3], [2,2])
+    B3 = BlockArray(randn(4,4), [1,3], [1,3])
     v = BlockArray(rand(4), [1,3])
     axB = axes(B)
+    axB2 = axes(B2)
     axv = axes(v)
     @test blockisequal(axB, axB)
+    @test blockisequal(axB[1], axB[1])
+    @test blockequals(axB[1], axB[1])
+    @test !blockisequal(axB[2], axB2[2])
+    @test !blockequals(axB[2], axB2[2])
+    @test !blockisequal(axB2[2], axB[2])
+    @test !blockequals(axB2[2], axB[2])
     @test blockisequal(axv, axv)
     @test !blockisequal(axB, axv)
     @test !blockisequal(axv, axB)
+
+    @test blockisequal(A, A)
+    @test blockisequal(B, B)
+    @test !blockisequal(B, A)
+    @test !blockisequal(A, B)
+    @test !blockisequal(B, B2)
+    @test !blockisequal(B2, B)
+    @test !blockisequal(B, B3)
+    @test !blockisequal(B3, B)
+end
+
+@testset "blockisapprox" begin
+    A1 = reshape([1:16;], (4,4))
+    A2 = A1 .+ 1e-10
+    B1 = BlockArray(A1, [1,3], [1,3])
+    B2 = BlockArray(A2, [1,3], [1,3])
+    B12 = BlockArray(A1, [1,3], [2,2])
+    @test !blockisapprox(A1, A2; rtol=0)
+    @test blockisapprox(A1, A2; rtol=1e-9)
+    @test !blockisapprox(A1, B1; rtol=1e-9)
+    @test !blockisapprox(B1, B2; rtol=0)
+    @test blockisapprox(B1, B2; rtol=1e-9)
+    @test !blockisapprox(B1, B12; rtol=0)
+    @test !blockisapprox(B1, B12; rtol=1e-9)
 end
 
 @testset "BlockIndices" begin