(Юлия). + Оператор, похоже, не использует мою функцию широковещания, почему? - PullRequest
2 голосов
/ 04 ноября 2019

Я реализую пользовательскую матрицу, которая имеет только одно ненулевое значение, и независимо от того, какую операцию вы выполняете, это единственная ячейка в матрице, которая может быть ненулевой. Я назвал это SVMatrix (однозначная матрица). Пока у меня есть код

struct SVMatrix{T} <: Base.AbstractMatrix{T}
    value::T
    index::Tuple{Int,Int}
    size::Tuple{Int,Int}
end

function Base.broadcast(+, A::SVMatrix, B::AbstractArray)
    SVMatrix(A.value+B[A.index...], A.index, A.size) 
end

function Base.getindex(A::SVMatrix{T}, i::Int) where {T}
    if i == A.index[1] + A.index[2]*A.size[1]
        A.value
    else
        0
    end
end


function Base.getindex(A::SVMatrix{T}, i::Vararg{Int,2}) where {T}
    if i == A.index
        return A.value
    else
        0
    end
end

function Base.size(A::SVMatrix)
    A.size
end

Затем я синхронизировал функцию широковещания вместе с оператором. + Следующим образом

function time(n::Int)
    A = SVMatrix(1.0, (3,4), (n, n))
    B = rand(n,n)
    @time broadcast(+, A, B)
    @time A .+ B
end

time(1000)
println()
time(1000)

и получил результаты

 0.000000 seconds
 0.008207 seconds (2 allocations: 7.629 MiB, 47.51% gc time)

 0.000000 seconds
 0.008258 seconds (2 allocations: 7.629 MiB)

Похоже, что. + Не использует мою функцию широковещания, хотя в документации говорится, что

Фактически, f. (Args ... ) эквивалентен широковещанию (f, args ...), предоставляя удобный синтаксис для широковещательной передачи любой функции (точечный синтаксис).

Почему я получаю эти результаты?

1 Ответ

5 голосов
/ 04 ноября 2019

На самом деле это отличный пример того, как вы не должны расширять трансляцию.

julia> struct SVMatrix{T} <: Base.AbstractMatrix{T}
           value::T
           index::Tuple{Int,Int}
           size::Tuple{Int,Int}
       end

julia> @inline function Base.getindex(A::SVMatrix{T}, i::Vararg{Int,2}) where {T}
           @boundscheck checkbounds(A, i...)
           if i == A.index
               return A.value
           else
               return zero(T)
           end
       end

julia> Base.size(A::SVMatrix) = A.size

julia> SVMatrix(1.0, (1,1), (2, 2)) .+ ones(2, 2)
2×2 Array{Float64,2}:
 2.0  1.0
 1.0  1.0

Результат .+ должен , а не быть [2 0; 0 0]! Если бы мы использовали вашу реализацию широковещания (исправлено для отправки на ::typeof(+), как DNF отметил ), ваш массив был бы удивительно сломан, когда другие использовали его и ожидали, что он будет вести себя как все другие AbstractArray s.

Теперь операция, в которой вы могли бы возвратить умно пересчитанные значения SVMatrix, равна .*:

julia> SVMatrix(2.5, (1,1), (2, 2)) .* ones(2, 2)
2×2 Array{Float64,2}:
 2.5  0.0
 0.0  0.0

Мы можем выполнить эту операцию в O (1) пространство и время, но реализация по умолчанию зацикливается на всех значениях и возвращает плотный Array. Вот где сияет многократная диспетчеризация Джулии:

julia> Base.broadcasted(::typeof(*), A::SVMatrix, B::AbstractArray) = SVMatrix(A.value*B[A.index...], A.index, A.size)

julia> SVMatrix(2.5, (1,1), (2, 2)) .* ones(2, 2)
2×2 SVMatrix{Float64}:
 2.5  0.0
 0.0  0.0

Поскольку это операция O (1) и является огромным выигрышем , мы можем отказаться от широковещательного слияния и немедленнопересчитать новый SVMatrix - даже внутри "слитого" выражения. Здесь вы еще не закончили!

  • Необходимо реализовать проверку ошибок для совместимых фигур.
  • Необходимо разрешить трансляцию таких вещей, как SVMatrix(2.5, (1,1), (2, 2)) .* rand(2).
  • В идеале вы должны реализовать BroadcastStyle, чтобы разрешить отправку для "хотя бы одного SVMatrix в списке аргументов". Затем вы реализуете Base.broadcasted(::ArrayStyle{SVMatrix}, ::typeof(*), args...), что позволит SVMatrix появляться по обе стороны от .*, но это сложная тема.
...