Merge pull request #5810 from JuliaLang/anj/arraydiv

Restrict +,- and / to algebraic operations (for matrices). Use . versions for elementwise operations. Add UniformScaling matrix.

NEWS.md

@@ -156,6 +156,10 @@ Library improvements
         the same length.  This generalizes and replaces `normfro` ([#6057]),
         and `norm` is now type-stable ([#6056]).
 
+      * + and - now only works when the sizes of the arrays are the same, i.e. the 
+        operations no longer do broadcasting. New `UniformScaling` type and identity
+        `I` constant (#5810).
+
     * Sparse linear algebra
 
       * Faster sparse `kron` ([#4958]).

base/abstractarray.jl

@@ -326,11 +326,9 @@ imag{T<:Real}(x::AbstractArray{T}) = zero(x)
 *(A::Number, B::AbstractArray) = A .* B
 *(A::AbstractArray, B::Number) = A .* B
 
-/(A::Number, B::AbstractArray) = A ./ B
 /(A::AbstractArray, B::Number) = A ./ B
 
 \(A::Number, B::AbstractArray) = B ./ A
-\(A::AbstractArray, B::Number) = B ./ A
 
 ./(x::Number,y::AbstractArray )         = throw(MethodError(./, (x,y)))
 ./(x::AbstractArray, y::Number)         = throw(MethodError(./, (x,y)))
@@ -398,7 +396,7 @@ function circshift(a, shiftamts)
     for i=1:n
         s = size(a,i)
         d = i<=length(shiftamts) ? shiftamts[i] : 0
-        I[i] = d==0 ? (1:s) : mod([-d:s-1-d], s)+1
+        I[i] = d==0 ? (1:s) : mod([-d:s-1-d], s).+1
     end
     a[I...]::typeof(a)
 end
@@ -1184,10 +1182,6 @@ function mapslices(f::Function, A::AbstractArray, dims::AbstractVector)
     ndimsA = ndims(A)
     alldims = [1:ndimsA]
 
-    if dims == alldims
-        return f(A)
-    end
-
     otherdims = setdiff(alldims, dims)
 
     idx = cell(ndimsA)

base/array.jl

@@ -782,7 +782,7 @@ for f in (:+, :-, :div, :mod, :&, :|, :$)
         end
     end
 end
-for f in (:+, :-, :.*, :./, :.%, :div, :mod, :rem, :&, :|, :$)
+for f in (:.+, :.-, :.*, :./, :.%, :div, :mod, :rem, :&, :|, :$)
     @eval begin
         function ($f){T}(A::Number, B::StridedArray{T})
             F = similar(B, promote_array_type(typeof(A),T))
@@ -802,7 +802,7 @@ for f in (:+, :-, :.*, :./, :.%, :div, :mod, :rem, :&, :|, :$)
 end
 
 # functions that should give an Int result for Bool arrays
-for f in (:+, :-)
+for f in (:.+, :.-)
     @eval begin
         function ($f)(A::Bool, B::StridedArray{Bool})
             F = Array(Int, size(B))
@@ -818,12 +818,16 @@ for f in (:+, :-)
             end
             return F
         end
+    end
+end
+for f in (:+, :-)
+    @eval begin
         function ($f)(A::StridedArray{Bool}, B::StridedArray{Bool})
             F = Array(Int, promote_shape(size(A), size(B)))
             for i=1:length(A)
                 @inbounds F[i] = ($f)(A[i], B[i])
             end
-            return F
+            return F        
         end
     end
 end

base/bitarray.jl

@@ -898,7 +898,7 @@ for f in (:+, :-)
         return r
     end
 end
-for f in (:+, :-),
+for f in (:.+, :.-),
     (arg1, arg2, T, fargs) in ((:(B::BitArray), :(x::Bool)    , Int                                   , :(b, x)),
                                (:(B::BitArray), :(x::Number)  , :(promote_array_type(typeof(x), Bool)), :(b, x)),
                                (:(x::Bool)    , :(B::BitArray), Int                                   , :(x, b)),

base/deprecated.jl

@@ -189,6 +189,16 @@ export PipeString
 @deprecate svdfact(A,thin)      svdfact(A,thin=thin)
 @deprecate svdfact!(A,thin)     svdfact(A,thin=thin)
 @deprecate svd(A,thin)          svd(A,thin=thin)
+@deprecate (+)(A::Array{Bool},x::Bool)      A .+ x
+@deprecate (+)(x::Bool,A::Array{Bool})      x .+ A 
+@deprecate (-)(A::Array{Bool},x::Bool)      A .- x
+@deprecate (-)(x::Bool,A::Array{Bool})      x .- A
+@deprecate (+)(A::Array,x::Number)          A .+ x
+@deprecate (+)(x::Number,A::Array)          x .+ A
+@deprecate (-)(A::Array,x::Number)          A .- x
+@deprecate (-)(x::Number,A::Array)          x .- A
+@deprecate (/)(x::Number,A::Array)          x ./ A
+@deprecate (\)(A::Array,x::Number)          A .\ x
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)

base/exports.jl

@@ -51,6 +51,7 @@ export
     GeneralizedSVD,
     Hermitian,
     Hessenberg,
+    UniformScaling,
     InsertionSort,
     IntSet,
     IO,
@@ -189,6 +190,7 @@ export
     γ, eulergamma,
     catalan,
     φ, golden,
+    I,
 
 # Operators
     !,

base/linalg.jl

@@ -33,6 +33,7 @@ export
     Symmetric,
     Triangular,
     Diagonal,
+    UniformScaling,
 
 # Functions
     axpy!,
@@ -145,7 +146,10 @@ export
     At_mul_Bt,
     At_mul_Bt!,
     At_rdiv_B,
-    At_rdiv_Bt
+    At_rdiv_Bt,
+
+# Constants
+    I
 
 typealias BlasFloat Union(Float64,Float32,Complex128,Complex64)
 typealias BlasReal Union(Float64,Float32)
@@ -201,6 +205,7 @@ include("linalg/woodbury.jl")
 include("linalg/tridiag.jl")
 include("linalg/diagonal.jl")
 include("linalg/bidiag.jl")
+include("linalg/uniformscaling.jl")
 include("linalg/rectfullpacked.jl")
 include("linalg/givens.jl")
 include("linalg/special.jl")

base/linalg/factorization.jl

@@ -248,7 +248,7 @@ Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular
 
 /{T}(B::Matrix{T},A::LU{T}) = At_ldiv_Bt(A,B).'
 
-inv{T<:BlasFloat}(A::LU{T})=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
+inv{T<:BlasFloat}(A::LU{T}) = @assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
 
 cond{T<:BlasFloat}(A::LU{T}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
 cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)

base/linalg/generic.jl

@@ -212,12 +212,14 @@ trace(x::Number) = x
 inv(a::AbstractVector) = error("argument must be a square matrix")
 inv{T}(A::AbstractMatrix{T}) = A_ldiv_B!(A,eye(T, chksquare(A)))
 
-function \{TA<:Number,TB<:Number}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB})
+function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB})
     TC = typeof(one(TA)/one(TB))
     A_ldiv_B!(convert(typeof(A).name.primary{TC}, A), TB == TC ? copy(B) : convert(typeof(B).name.primary{TC}, B))
 end
 \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b
 /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'
+# \(A::StridedMatrix,x::Number) = inv(A)*x Should be added at some point when the old elementwise version has been deprecated long enough
+# /(x::Number,A::StridedMatrix) = x*inv(A)
 
 cond(x::Number) = x == 0 ? Inf : 1.0
 cond(x::Number, p) = cond(x)

base/linalg/uniformscaling.jl

@@ -0,0 +1,74 @@
+import Base: +, -, *, /, copy, ctranspose, getindex, showarray, transpose
+import Base.LinAlg: SingularException
+immutable UniformScaling{T<:Number} <: AbstractMatrix{T}
+		λ::T
+end
+
+const I = UniformScaling(1)
+
+getindex(J::UniformScaling, i::Integer,j::Integer) = ifelse(i==j,J.λ,zero(J.λ))
+
+showarray(io::IO,J::UniformScaling;kw...) = print(io,"$(typeof(J))\n$(J.λ)*I")
+copy(J::UniformScaling) = UniformScaling(J.λ)
+
+transpose(J::UniformScaling) = J
+ctranspose(J::UniformScaling) = UniformScaling(conj(J.λ))
+
++(J1::UniformScaling,J2::UniformScaling)  = UniformScaling(J1.λ+J2.λ)
++{T}(B::BitArray{2},J::UniformScaling{T}) = bitunpack(B) + J
++(J::UniformScaling,B::BitArray{2}) 	  = J + bitunpack(B)
+function +{TA,TJ}(A::AbstractMatrix{TA},J::UniformScaling{TJ})
+    n = chksquare(A)
+    B = similar(A,promote_type(TA,TJ))
+    copy!(B,A)
+    @inbounds for i = 1:n
+        B[i,i] += J.λ
+    end
+    B
+end
++(J::UniformScaling,A::AbstractMatrix) = A + J
+
+-(J1::UniformScaling,J2::UniformScaling) = UniformScaling(J1.λ-J2.λ)
+-(B::BitArray{2},J::UniformScaling) 	 = bitunpack(B) - J
+-(J::UniformScaling,B::BitArray{2}) 	 = J - bitunpack(B)
+function -{TA,TJ<:Number}(A::AbstractMatrix{TA},J::UniformScaling{TJ})
+    n = chksquare(A)
+    B = similar(A,promote_type(TA,TJ))
+    copy!(B,A)
+    @inbounds for i = 1:n
+        B[i,i] -= J.λ
+    end
+    B
+end
+function -{TA,TJ<:Number}(J::UniformScaling{TJ},A::AbstractMatrix{TA})
+	n = chksquare(A)
+	B = -A
+	@inbounds for i = 1:n
+		B[i,i] += J.λ
+	end
+	B
+end
+
+*(J1::UniformScaling,J2::UniformScaling) = UniformScaling(J1.λ*J2.λ)
+*(B::BitArray{2},J::UniformScaling) = *(bitunpack(B),J::UniformScaling)
+*(J::UniformScaling,B::BitArray{2}) = *(J::UniformScaling,bitunpack(B))
+*(S::SparseMatrixCSC,J::UniformScaling) = J.λ == 1 ? S : J.λ*S
+*{Tv,Ti}(J::UniformScaling,S::SparseMatrixCSC{Tv,Ti}) = J.λ == 1 ? S : S*J.λ
+*(A::AbstractMatrix,J::UniformScaling) = J.λ == 1 ? A : J.λ*A
+*(J::UniformScaling,A::AbstractVecOrMat) = J.λ == 1 ? A : J.λ*A
+
+*(x::Number,J::UniformScaling) = UniformScaling(x*J.λ)
+*(J::UniformScaling,x::Number) = UniformScaling(J.λ*x)
+
+/(J1::UniformScaling,J2::UniformScaling) = J2.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ/J2.λ)
+/(J::UniformScaling,A::AbstractMatrix) = J.λ*inv(A)
+/(A::AbstractMatrix,J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : A/J.λ
+
+/(J::UniformScaling,x::Number) = UniformScaling(J.λ/x)
+
+\(J1::UniformScaling,J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ)
+\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}),J::UniformScaling) = inv(A)*J.λ
+\(J::UniformScaling,A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A
+\(A::AbstractMatrix,J::UniformScaling) = inv(A)*J.λ
+
+\(x::Number,J::UniformScaling) = UniformScaling(x\J.λ)

base/sparse/sparsematrix.jl

@@ -1229,7 +1229,7 @@ function vcat(X::SparseMatrixCSC...)
             rX2 = X[i].colptr[c + 1] - 1
             rr2 = rr1 + (rX2 - rX1)
 
-            rowval[rr1 : rr2] = X[i].rowval[rX1 : rX2] + mX_sofar
+            rowval[rr1 : rr2] = X[i].rowval[rX1 : rX2] .+ mX_sofar
             nzval[rr1 : rr2] = X[i].nzval[rX1 : rX2]
             mX_sofar += mX[i]
             rr1 = rr2 + 1
@@ -1261,7 +1261,7 @@ function hcat(X::SparseMatrixCSC...)
     nnz_sofar = 0
     nX_sofar = 0
     for i = 1 : num
-        colptr[(1 : nX[i] + 1) + nX_sofar] = X[i].colptr + nnz_sofar
+        colptr[(1 : nX[i] + 1) + nX_sofar] = X[i].colptr .+ nnz_sofar
         rowval[(1 : nnzX[i]) + nnz_sofar] = X[i].rowval
         nzval[(1 : nnzX[i]) + nnz_sofar] = X[i].nzval
         nnz_sofar += nnzX[i]
@@ -1303,8 +1303,8 @@ function blkdiag(X::SparseMatrixCSC...)
     nX_sofar = 0
     mX_sofar = 0
     for i = 1 : num
-        colptr[(1 : nX[i] + 1) + nX_sofar] = X[i].colptr + nnz_sofar
-        rowval[(1 : nnzX[i]) + nnz_sofar] = X[i].rowval + mX_sofar
+        colptr[(1 : nX[i] + 1) + nX_sofar] = X[i].colptr .+ nnz_sofar
+        rowval[(1 : nnzX[i]) + nnz_sofar] = X[i].rowval .+ mX_sofar
         nzval[(1 : nnzX[i]) + nnz_sofar] = X[i].nzval
         nnz_sofar += nnzX[i]
         nX_sofar += nX[i]

base/statistics.jl

@@ -274,15 +274,15 @@ function quantile!(v::AbstractVector, q::AbstractVector)
     lv = length(v)
     lq = length(q)
 
-    index = 1 + (lv-1)*q
+    index = 1 .+ (lv-1)*q
     lo = ifloor(index)
     hi = iceil(index)
     sort!(v)
     isnan(v[end]) && error("quantiles are undefined in presence of NaNs")
     i = find(index .> lo)
     r = float(v[lo])
-    h = (index-lo)[i]
-    r[i] = (1-h).*r[i] + h.*v[hi[i]]
+    h = (index.-lo)[i]
+    r[i] = (1.-h).*r[i] + h.*v[hi[i]]
     return r
 end
 quantile(v::AbstractVector, q::AbstractVector) = quantile!(copy(v),q)

doc/manual/arrays.rst

@@ -286,7 +286,7 @@ operator should be used for elementwise operations.
 5.  Binary Boolean or bitwise — ``&``, ``|``, ``$``
 
 Some operators without dots operate elementwise anyway when one argument is a
-scalar. These operators are ``+``, ``-``, ``*``, ``/``, ``\``, and the bitwise
+scalar. These operators are ``*``, ``/``, ``\``, and the bitwise
 operators.
 
 Note that comparisons such as ``==`` operate on whole arrays, giving a single

doc/manual/linear-algebra.rst

@@ -51,7 +51,8 @@ for them in LAPACK are available.
 +--------------------+-----------------------------------------------------------------------------------+
 | ``Diagonal``       | `Diagonal matrix <http://en.wikipedia.org/wiki/Diagonal_matrix>`_                 |
 +--------------------+-----------------------------------------------------------------------------------+
-
+| ``UniformScaling`` | `Uniform scaling operator <http://en.wikipedia.org/wiki/Uniform_scaling>`_        |
++--------------------+-----------------------------------------------------------------------------------+
 
 Elementary operations
 ---------------------
@@ -74,6 +75,8 @@ Elementary operations
 | ``Diagonal``       |   X   |   X   |  XY   |   XY  | ``inv``, ``det``,   |
 |                    |       |       |       |       | ``logdet``, ``/``   |
 +--------------------+-------+-------+-------+-------+---------------------+
+| ``UniformScaling`` |   X   |   X   |  XYZ  |  XYZ  | ``/``               |
++--------------------+-------+-------+-------+-------+---------------------+
 
 Legend:
 
@@ -116,3 +119,7 @@ Legend:
 | D | An optimized method to find the characteristic vectors corresponding to the characteristic values ``x=[x1, x2,...]`` is available | ``eigvecs(M, x)``      |
 +---+-----------------------------------------------------------------------------------------------------------------------------------+------------------------+
 
+The uniform scaling operator
+--------------------------
+A ``UniformScaling`` operator represents a scalar times the identity operator, ``λ*I``. The identity operator ``I`` is defined as a constant and is an instance of ``UniformScaling``. The size of these operators are generic and match the other matrix in the binary operations ``+``,``-``,``*`` and ``\``. For ``A+I`` and ``A-I`` this means that ``A`` must be square. Multiplication with the identity operator ``I`` is a noop (except for checking that the scaling factor is one) and therefore almost without overhead. 
+

test/arrayops.jl

@@ -12,14 +12,14 @@ b = a+a
 @test length((1,)) == 1
 @test length((1,2)) == 2
 
-@test isequal(1+[1,2,3], [2,3,4])
-@test isequal([1,2,3]+1, [2,3,4])
-@test isequal(1-[1,2,3], [0,-1,-2])
-@test isequal([1,2,3]-1, [0,1,2])
+@test isequal(1.+[1,2,3], [2,3,4])
+@test isequal([1,2,3].+1, [2,3,4])
+@test isequal(1.-[1,2,3], [0,-1,-2])
+@test isequal([1,2,3].-1, [0,1,2])
 
 @test isequal(5*[1,2,3], [5,10,15])
 @test isequal([1,2,3]*5, [5,10,15])
-@test isequal(1/[1,2,5], [1.0,0.5,0.2])
+@test isequal(1./[1,2,5], [1.0,0.5,0.2])
 @test isequal([1,2,3]/5, [0.2,0.4,0.6])
 
 a = ones(2,2)
@@ -703,10 +703,9 @@ begin
     @test all(b.==6)
 
     # issue #5141
+    ## Update Removed the version that removes the dimensions when dims==1:ndims(A) 
     c1 = mapslices(x-> maximum(-x), a, [])
     @test c1 == -a
-    c2 = mapslices(x-> maximum(-x), a, [1,2])
-    @test c2 == maximum(-a)
 
     # other types than Number
     @test mapslices(prod,["1" "2"; "3" "4"],1) == ["13" "24"]
@@ -859,7 +858,7 @@ end
 @test isequal(flipdim([2,3,1], 2), [2,3,1])
 @test isequal(flipdim([2 3 1], 1), [2 3 1])
 @test isequal(flipdim([2 3 1], 2), [1 3 2])
-@test_throws flipdim([2,3,1] -1)
+@test_throws flipdim([2,3,1], -1)
 @test isequal(flipdim(1:10, 1), 10:-1:1)
 @test isequal(flipdim(1:10, 2), 1:10)
 @test_throws flipdim(1:10, -1)