Add more SubstitutionMatrix operations (#60)

SubstitutionMatrix is not a subtype of AbstractArray.
This may or may not be intentional, but currently, the number of supported operations is fairly minimal.
This provides a few more conveniences for users who might want to manipulate them beyond mere lookup tables.

Author: timholy

docs/src/pairalign.md

@@ -5,6 +5,7 @@ CurrentModule = BioAlignments
 DocTestSetup = quote
     using BioSequences
     using BioAlignments
+    using BioAlignments.BioSymbols
 end
 ```
 
@@ -237,6 +238,15 @@ Make sure to create a copy of the original matrix when you create a matrix from
 a predefined matrix. In the above case, `BLOSUM62` is shared in the whole
 program and modification on it will affect any result that uses `BLOSUM62`.
 
+Other supported operations include `similar`, `isassigned`, and `BioSymbols.alphabet`:
+```jldoctest
+julia> reshape(alphabet(BLOSUM62), 1, :)
+1×27 Matrix{AminoAcid}:
+ AA_A  AA_R  AA_N  AA_D  AA_C  AA_Q  AA_E  …  AA_B  AA_J  AA_Z  AA_X  AA_Term
+```
+
+As shown, there is no entry corresponding to `AA_Gap`.
+
 `DichotomousSubstitutionMatrix` is a specialized matrix for matching or
 mismatching substitution.  This is a preferable choice when performance is more
 important than flexibility because looking up score is faster than

src/submat.jl

@@ -13,7 +13,7 @@ The required method:
 
 * `Base.getindex(submat, x, y)`: substitution score/cost from `x` to `y`
 """
-abstract type AbstractSubstitutionMatrix{S<:Real} end
+abstract type AbstractSubstitutionMatrix{S<:Real} end # is it deliberate that this is not <:AbstractMatrix{S}?
 
 """
 Substitution matrix.
@@ -82,18 +82,37 @@ Base.convert(::Type{Matrix}, submat::SubstitutionMatrix) = copy(submat.data)
     return submat.data[i,j]
 end
 
+Base.getindex(submat::SubstitutionMatrix, mask::AbstractMatrix{Bool}) = submat.data[mask]
+Base.getindex(submat::SubstitutionMatrix, mask1::AbstractVector{Bool}, mask2::AbstractVector{Bool}) = submat.data[mask1, mask2]
+
 function Base.setindex!(submat::SubstitutionMatrix{T}, val, x, y) where T
     i = index(convert(T, x))
     j = index(convert(T, y))
     submat.data[i,j] = val
     submat.defined[i,j] = true
-    return submat
+    return val
+end
+
+Base.eltype(::SubstitutionMatrix{T,S}) where {T,S} = S
+
+BioSymbols.alphabet(::SubstitutionMatrix{T}) where {T} = alphabet_without_gap(T)
+
+Base.isassigned(submat::SubstitutionMatrix) = submat.defined
+function Base.isassigned(submat::SubstitutionMatrix{T}, x, y) where T
+    i = index(convert(T, x))
+    j = index(convert(T, y))
+    return submat.defined[i,j]
 end
 
 function Base.copy(submat::SubstitutionMatrix{T,S}) where {T,S}
     return SubstitutionMatrix{T,S}(copy(submat.data), copy(submat.defined))
 end
 
+function Base.similar(submat::SubstitutionMatrix{T}, ::Type{S}=eltype(submat)) where {T,S<:Real}
+    n = length(BioSymbols.alphabet(T)) - 1
+    return SubstitutionMatrix{T,S}(Matrix{S}(undef, n, n), falses(n, n))
+end
+
 Base.minimum(submat::SubstitutionMatrix) = minimum(submat.data)
 Base.maximum(submat::SubstitutionMatrix) = maximum(submat.data)
 

test/runtests.jl

@@ -376,6 +376,20 @@ end
 
         @test convert(Matrix, BioAlignments.load_submat(AminoAcid, "BLOSUM62")) == convert(Matrix, BLOSUM62)
 
+        # Various operations
+        @test eltype(BLOSUM62) === Int
+        @test eltype(similar(BLOSUM62, Float64)) === Float64
+        @test !any(isassigned(similar(BLOSUM62)))
+        @test isassigned(BLOSUM62, AA_R, AA_K)
+        @test !isassigned(BLOSUM62, AA_R, AA_O)
+        @test size(BLOSUM62[isassigned(BLOSUM62)]) == (sum(BLOSUM62.defined),)
+        keep1, keep2 = vec(all(isassigned(BLOSUM62); dims=2)), vec(all(isassigned(BLOSUM62); dims=1))
+        @test size(BLOSUM62[keep1, keep2]) == (sum(keep1), sum(keep2))
+        @test BioSymbols.alphabet(BLOSUM62) == BioAlignments.alphabet_without_gap(AminoAcid)
+        myblossum = similar(BLOSUM62, Float64)
+        myblossum[AA_K,AA_R] = myblossum[AA_R,AA_K] = 1.11
+        @test myblossum[AA_K,AA_R] == myblossum[AA_R,AA_K] == 1.11
+
         submat = SubstitutionMatrix(DNA, rand(Float64, 15, 15))
         @test isa(submat, SubstitutionMatrix{DNA,Float64})