Fill multiplication (#70)

* Add Fill multiplication special cases

* v0.7

* Drop Julia v0.7 support

* Add Eye tests

* Improved ones broadcasting

Author: dlfivefifty

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
   - osx
 julia:
-  - 0.7
   - 1.0
   - 1.1
   - 1.2

Project.toml

@@ -1,14 +1,14 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.6.4"
+version = "0.7"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-julia = "0.7, 1"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

appveyor.yml

@@ -1,6 +1,5 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1
   - julia_version: 1.1
   - julia_version: 1.2

src/fillalgebra.jl

@@ -15,23 +15,72 @@ adjoint(a::Fill{T,2}) where T = Fill{T}(adjoint(a.value), reverse(a.axes))
 
 ## Algebraic identities
 
+
+function mult_fill(a::AbstractFill, b::AbstractFill{<:Any,2})
+    axes(a, 2) ≠ axes(b, 1) &&
+        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
+    return Fill(getindex_value(a)*getindex_value(b), (axes(a, 1), axes(b, 2)))
+end
+
+function mult_fill(a::AbstractFill, b::AbstractFill{<:Any,1})
+    axes(a, 2) ≠ axes(b, 1) &&
+        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
+    return Fill(getindex_value(a)*getindex_value(b), (axes(a, 1),))
+end
+
+function mult_ones(a, b::AbstractMatrix)
+    axes(a, 2) ≠ axes(b, 1) &&
+        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
+    return Ones{promote_type(eltype(a), eltype(b))}((axes(a, 1), axes(b, 2)))
+end
+function mult_ones(a, b::AbstractVector)
+    axes(a, 2) ≠ axes(b, 1) &&
+        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
+    return Ones{promote_type(eltype(a), eltype(b))}((axes(a, 1),))
+end
+
 function mult_zeros(a, b::AbstractMatrix)
-    size(a, 2) ≠ size(b, 1) &&
+    axes(a, 2) ≠ axes(b, 1) &&
         throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
-    return Zeros{promote_type(eltype(a), eltype(b))}(size(a, 1), size(b, 2))
+    return Zeros{promote_type(eltype(a), eltype(b))}((axes(a, 1), axes(b, 2)))
 end
 function mult_zeros(a, b::AbstractVector)
-    size(a, 2) ≠ size(b, 1) &&
+    axes(a, 2) ≠ axes(b, 1) &&
         throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
-    return Zeros{promote_type(eltype(a), eltype(b))}(size(a, 1))
+    return Zeros{promote_type(eltype(a), eltype(b))}((axes(a, 1),))
 end
 
-const ZerosVecOrMat{T} = Union{Zeros{T,1}, Zeros{T,2}}
-*(a::ZerosVecOrMat, b::AbstractMatrix) = mult_zeros(a, b)
-*(a::AbstractMatrix, b::ZerosVecOrMat) = mult_zeros(a, b)
-*(a::ZerosVecOrMat, b::AbstractVector) = mult_zeros(a, b)
-*(a::AbstractVector, b::ZerosVecOrMat) = mult_zeros(a, b)
-*(a::ZerosVecOrMat, b::ZerosVecOrMat) = mult_zeros(a, b)
+*(a::AbstractFill{<:Any,1}, b::AbstractFill{<:Any,2}) = mult_fill(a,b)
+*(a::AbstractFill{<:Any,2}, b::AbstractFill{<:Any,2}) = mult_fill(a,b)
+*(a::AbstractFill{<:Any,2}, b::AbstractFill{<:Any,1}) = mult_fill(a,b)
+
+*(a::Ones{<:Any,1}, b::Ones{<:Any,2}) = mult_ones(a, b)
+*(a::Ones{<:Any,2}, b::Ones{<:Any,2}) = mult_ones(a, b)
+*(a::Ones{<:Any,2}, b::Ones{<:Any,1}) = mult_ones(a, b)
+
+*(a::Zeros{<:Any,1}, b::Zeros{<:Any,2}) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::Zeros{<:Any,2}) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::Zeros{<:Any,1}) = mult_zeros(a, b)
+
+*(a::Zeros{<:Any,1}, b::AbstractFill{<:Any,2}) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::AbstractFill{<:Any,2}) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::AbstractFill{<:Any,1}) = mult_zeros(a, b)
+*(a::AbstractFill{<:Any,1}, b::Zeros{<:Any,2}) = mult_zeros(a,b)
+*(a::AbstractFill{<:Any,2}, b::Zeros{<:Any,2}) = mult_zeros(a,b)
+*(a::AbstractFill{<:Any,2}, b::Zeros{<:Any,1}) = mult_zeros(a,b)
+
+*(a::Zeros{<:Any,1}, b::AbstractMatrix) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::AbstractMatrix) = mult_zeros(a, b)
+*(a::AbstractMatrix, b::Zeros{<:Any,1}) = mult_zeros(a, b)
+*(a::AbstractMatrix, b::Zeros{<:Any,2}) = mult_zeros(a, b)
+*(a::Zeros{<:Any,1}, b::AbstractVector) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::AbstractVector) = mult_zeros(a, b)
+*(a::AbstractVector, b::Zeros{<:Any,2}) = mult_zeros(a, b)
+
+*(a::Zeros{<:Any,1}, b::Diagonal) = mult_zeros(a, b)
+*(a::Zeros{<:Any,2}, b::Diagonal) = mult_zeros(a, b)
+*(a::Diagonal, b::Zeros{<:Any,1}) = mult_zeros(a, b)
+*(a::Diagonal, b::Zeros{<:Any,2}) = mult_zeros(a, b)
 
 
 function *(a::Adjoint{T, <:AbstractVector{T}}, b::Zeros{S, 1}) where {T, S}

src/fillbroadcast.jl

@@ -28,34 +28,49 @@ broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::Zeros) = _broadcasted
 
 broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::Ones) = _broadcasted_zeros(a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::Fill) = _broadcasted_zeros(a, b)
-function broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::AbstractRange)
+broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::AbstractRange) =
     return _broadcasted_zeros(a, b)
-end
-function broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::AbstractArray)
+broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::AbstractArray) =
     return _broadcasted_zeros(a, b)
-end
 
 broadcasted(::DefaultArrayStyle, ::typeof(*), a::Ones, b::Zeros) = _broadcasted_zeros(a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(*), a::Fill, b::Zeros) = _broadcasted_zeros(a, b)
-function broadcasted(::DefaultArrayStyle, ::typeof(*), a::AbstractRange, b::Zeros)
+broadcasted(::DefaultArrayStyle, ::typeof(*), a::AbstractRange, b::Zeros) =
     return _broadcasted_zeros(a, b)
-end
-function broadcasted(::DefaultArrayStyle, ::typeof(*), a::AbstractArray, b::Zeros)
+broadcasted(::DefaultArrayStyle, ::typeof(*), a::AbstractArray, b::Zeros) =
     return _broadcasted_zeros(a, b)
-end
+broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::Zeros, b::AbstractRange) =
+    return _broadcasted_zeros(a, b)
+broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractRange, b::Zeros) =
+    return _broadcasted_zeros(a, b)
+
 
 broadcasted(::DefaultArrayStyle, ::typeof(*), a::Ones, b::Ones) = _broadcasted_ones(a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(/), a::Ones, b::Ones) = _broadcasted_ones(a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(\), a::Ones, b::Ones) = _broadcasted_ones(a, b)
 
+# special case due to missing converts for ranges
+_range_convert(::Type{AbstractVector{T}}, a::AbstractRange{T}) where T = a
+_range_convert(::Type{AbstractVector{T}}, a::AbstractUnitRange) where T = convert(T,first(a)):convert(T,last(a))
+_range_convert(::Type{AbstractVector{T}}, a::AbstractRange) where T = convert(T,first(a)):step(a):convert(T,last(a))
+
+function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::Ones{T}, b::AbstractRange{V}) where {T,V}
+    broadcast_shape(size(a), size(b)) # Check sizes are compatible.
+    return _range_convert(AbstractVector{promote_type(T,V)}, b)
+end
+
+function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractRange{V}, b::Ones{T}) where {T,V}
+    broadcast_shape(size(a), size(b)) # Check sizes are compatible.
+    return _range_convert(AbstractVector{promote_type(T,V)}, a)
+end
 
 
-function broadcasted(::DefaultArrayStyle, ::typeof(*), a::AbstractFill, b::AbstractRange)
+function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractFill, b::AbstractRange)
     broadcast_shape(size(a), size(b)) # Check sizes are compatible.
     return broadcasted(*, getindex_value(a), b)
 end
 
-function broadcasted(::DefaultArrayStyle, ::typeof(*), a::AbstractRange, b::AbstractFill)
+function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractRange, b::AbstractFill)
     broadcast_shape(size(a), size(b)) # Check sizes are compatible.
     return broadcasted(*, a, getindex_value(b))
 end

test/runtests.jl

@@ -150,7 +150,6 @@ import FillArrays: AbstractFill, RectDiagonal
     end
 end
 
-
 @testset "RectDiagonal" begin
     data = 1:3
     expected_size = (5, 3)
@@ -189,7 +188,6 @@ end
     @test_throws ArgumentError mut[2, 1] = 9
 end
 
-
 # Check that all pair-wise combinations of + / - elements of As and Bs yield the correct
 # type, and produce numerically correct results.
 function test_addition_and_subtraction(As, Bs, Tout::Type)
@@ -313,7 +311,6 @@ end
 
 end
 
-
 @testset "IndexStyle" begin
     @test IndexStyle(Zeros(5,5)) == IndexStyle(typeof(Zeros(5,5))) == IndexLinear()
 end
@@ -326,9 +323,9 @@ end
     @test Zeros(3, 4) * randn(4) == Zeros(3, 4) * Zeros(4) == Zeros(3)
     @test Zeros(3, 4) * Zeros(4, 5) === Zeros(3, 5)
 
-    @test [1,2,3]*Zeros(1) ≡ Zeros(3)
+    @test_throws MethodError [1,2,3]*Zeros(1) # Not defined for [1,2,3]*[0] either
     @test [1,2,3]*Zeros(1,3) ≡ Zeros(3,3)
-    @test_throws DimensionMismatch [1,2,3]*Zeros(3)
+    @test_throws MethodError [1,2,3]*Zeros(3) # Not defined for [1,2,3]*[0,0,0] either
 
     # Check multiplication by Adjoint vectors works as expected.
     @test randn(4, 3)' * Zeros(4) === Zeros(3)
@@ -503,6 +500,17 @@ end
     @test Fill(1,10) .- 1 ≡ Fill(1,10) .- Ref(1) ≡ Fill(1,10) .- Ref(1I)
     @test Fill([1 2; 3 4],10) .- Ref(1I) == Fill([0 2; 3 3],10)
     @test Ref(1I) .+ Fill([1 2; 3 4],10) == Fill([2 2; 3 5],10)
+
+    @testset "Special Ones" begin
+        @test Ones{Int}(5) .* (1:5) ≡ (1:5) .* Ones{Int}(5) ≡ 1:5 
+        @test Ones(5) .* (1:5) ≡ (1:5) .* Ones(5) ≡ 1.0:5 
+        @test Ones{Int}(5) .* Ones{Int}(5) ≡ Ones{Int}(5) 
+        @test Ones{Int}(5,2) .* (1:5) == Array(Ones{Int}(5,2)) .* Array(1:5)
+        @test (1:5) .* Ones{Int}(5,2)  == Array(1:5) .* Array(Ones{Int}(5,2)) 
+        @test_throws DimensionMismatch Ones{Int}(6) .* (1:5)
+        @test_throws DimensionMismatch (1:5) .* Ones{Int}(6)
+        @test_throws DimensionMismatch Ones{Int}(5) .* Ones{Int}(6)
+    end
 end
 
 @testset "Sub-arrays" begin
@@ -743,4 +751,71 @@ end
 
     @test fill!(F,1) == F
     @test_throws ArgumentError fill!(F,2)
-end
+end
+
+@testset "mult" begin
+    @test Fill(2,10)*Fill(3,1,12) == Vector(Fill(2,10))*Matrix(Fill(3,1,12))
+    @test Fill(2,10)*Fill(3,1,12) ≡ Fill(6,10,12)
+    @test Fill(2,3,10)*Fill(3,10,12) ≡ Fill(6,3,12)
+    @test Fill(2,3,10)*Fill(3,10) ≡ Fill(6,3)
+    @test_throws DimensionMismatch Fill(2,10)*Fill(3,2,12)
+    @test_throws DimensionMismatch Fill(2,3,10)*Fill(3,2,12)
+
+    @test Ones(10)*Fill(3,1,12) ≡ Fill(3.0,10,12)
+    @test Ones(10,3)*Fill(3,3,12) ≡ Fill(3.0,10,12)
+    @test Ones(10,3)*Fill(3,3) ≡ Fill(3.0,10)
+
+    @test Fill(2,10)*Ones(1,12) ≡ Fill(2.0,10,12)
+    @test Fill(2,3,10)*Ones(10,12) ≡ Fill(2.0,3,12)
+    @test Fill(2,3,10)*Ones(10) ≡ Fill(2.0,3)
+
+    @test Ones(10)*Ones(1,12) ≡ Ones(10,12)
+    @test Ones(3,10)*Ones(10,12) ≡ Ones(3,12)
+    @test Ones(3,10)*Ones(10) ≡ Ones(3)
+
+    @test Zeros(10)*Fill(3,1,12) ≡   Zeros(10,12)
+    @test Zeros(10,3)*Fill(3,3,12) ≡ Zeros(10,12)
+    @test Zeros(10,3)*Fill(3,3) ≡    Zeros(10)
+
+    @test Fill(2,10)*  Zeros(1,12) ≡  Zeros(10,12)
+    @test Fill(2,3,10)*Zeros(10,12) ≡ Zeros(3,12)
+    @test Fill(2,3,10)*Zeros(10) ≡    Zeros(3)
+
+    @test Zeros(10)*Zeros(1,12) ≡ Zeros(10,12)
+    @test Zeros(3,10)*Zeros(10,12) ≡ Zeros(3,12)
+    @test Zeros(3,10)*Zeros(10) ≡ Zeros(3)
+
+    a = randn(3)
+    A = randn(1,4)
+
+    @test Fill(2,3)*A == Vector(Fill(2,3))*A
+    @test Fill(2,3,1)*A == Matrix(Fill(2,3,1))*A
+    @test Fill(2,3,3)*a == Matrix(Fill(2,3,3))*a
+    @test Ones(3)*A ==   Vector(Ones(3))*A
+    @test Ones(3,1)*A == Matrix(Ones(3,1))*A
+    @test Ones(3,3)*a == Matrix(Ones(3,3))*a
+    @test Zeros(3)*A  ≡ Zeros(3,4)
+    @test Zeros(3,1)*A == Zeros(3,4)
+    @test Zeros(3,3)*a == Zeros(3)
+
+    @test A*Fill(2,4) == A*Vector(Fill(2,4))
+    @test A*Fill(2,4,1) == A*Matrix(Fill(2,4,1))
+    @test a*Fill(2,1,3) == a*Matrix(Fill(2,1,3))
+    @test A*Ones(4) ==   A*Vector(Ones(4))
+    @test A*Ones(4,1) == A*Matrix(Ones(4,1))
+    @test a*Ones(1,3) == a*Matrix(Ones(1,3))
+    @test A*Zeros(4)  ≡ Zeros(1)
+    @test A*Zeros(4,1) ≡ Zeros(1,1)
+    @test a*Zeros(1,3) ≡ Zeros(3,3)
+
+    D = Diagonal(randn(1))
+    @test Zeros(1,1)*D ≡ Zeros(1,1)
+    @test Zeros(1)*D ≡ Zeros(1,1)
+    @test D*Zeros(1,1) ≡ Zeros(1,1)
+    @test D*Zeros(1) ≡ Zeros(1)
+
+    E = Eye(5)
+    @test E*(1:5) ≡ 1.0:5.0
+    @test (1:5)'E == (1.0:5)'
+    @test E*E ≡ E
+end