Skip to content

Commit

Permalink
Remove Array+-Number methods. Define UniformScaling and the identity …
Browse files Browse the repository at this point in the history
…operator I. Let Number/Matrix be the inverse. Update docs.
  • Loading branch information
andreasnoack committed Mar 10, 2014
1 parent b871f72 commit e4c79b0
Show file tree
Hide file tree
Showing 20 changed files with 211 additions and 69 deletions.
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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]).
Expand Down
8 changes: 1 addition & 7 deletions base/abstractarray.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
10 changes: 7 additions & 3 deletions base/array.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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))
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion base/bitarray.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)),
Expand Down
10 changes: 10 additions & 0 deletions base/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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`)
Expand Down
2 changes: 2 additions & 0 deletions base/exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ export
GeneralizedSVD,
Hermitian,
Hessenberg,
UniformScaling,
InsertionSort,
IntSet,
IO,
Expand Down Expand Up @@ -189,6 +190,7 @@ export
γ, eulergamma,
catalan,
φ, golden,
I,

# Operators
!,
Expand Down
7 changes: 6 additions & 1 deletion base/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ export
Symmetric,
Triangular,
Diagonal,
UniformScaling,

# Functions
axpy!,
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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")
Expand Down
2 changes: 1 addition & 1 deletion base/linalg/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 3 additions & 1 deletion base/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
74 changes: 74 additions & 0 deletions base/linalg/uniformscaling.jl
Original file line number Diff line number Diff line change
@@ -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.λ)
8 changes: 4 additions & 4 deletions base/sparse/sparsematrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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]
Expand Down
6 changes: 3 additions & 3 deletions base/statistics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion doc/manual/arrays.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
9 changes: 8 additions & 1 deletion doc/manual/linear-algebra.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
---------------------
Expand All @@ -74,6 +75,8 @@ Elementary operations
| ``Diagonal`` | X | X | XY | XY | ``inv``, ``det``, |
| | | | | | ``logdet``, ``/`` |
+--------------------+-------+-------+-------+-------+---------------------+
| ``UniformScaling`` | X | X | XYZ | XYZ | ``/`` |
+--------------------+-------+-------+-------+-------+---------------------+

Legend:

Expand Down Expand Up @@ -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.

15 changes: 7 additions & 8 deletions test/arrayops.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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"]
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit e4c79b0

Please sign in to comment.