`diagm` doesn't work

[tkoolen]

A = OffsetArray(rand(3), -3)
diagm(A)

results in

ERROR: size not supported for arrays with indices (-2:0,); see http://docs.julialang.org/en/latest/devdocs/offset-arrays/
Stacktrace:
 [1] errmsg(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.6/OffsetArrays/src/OffsetArrays.jl:48
 [2] diagm at ./linalg/dense.jl:250 [inlined] (repeats 2 times)

on nightly. This should probably be fixed in Base, but I wanted to get your input first.

As for the desired output type of diagm, I guess it should always be

typeof(similar(A, first(indices(A)), first(indices(A))))

?

[alsam]

Thanks for the proposal.
Unfortunately zeros - the important integral part for diagm algorithm is deprecated.
The excerpt from dense.jl is below for reference:

function diagm{T}(v::AbstractVector{T}, k::Integer=0)
    n = length(v) + abs(k)
    A = zeros(T,n,n)
    A[diagind(A,k)] = v
    A
end

Please see #12 for more details.

[tkoolen]

The other important part I think is the indexing. diagind assumes one-indexed arrays, and zeros(T) also produces a one-indexed Array. Shouldn't the indexing of A match that of v?

The current diagm implementation in Base is also a problem for JuMP, where, for T == JuMP.Variable, typeof(zero(T)) == JuMP.AffExpr != JuMP.Variable.

A sketch of something that would solve both problems:

function diagm2{T}(v::AbstractVector{T})
    S = promote_type(T, typeof(zero(T)))
    A = similar(v, S, first(indices(v)), first(indices(v)))
    for i in eachindex(v), j in eachindex(v)
        A[i, j] = i == j ? v[i] : zero(T)
    end
    A
end

with which

using OffsetArrays, JuMP
m = Model()
@variable(m, v[1 : 3])
diagm2(OffsetArray(v, -3))

produces

OffsetArrays.OffsetArray{JuMP.GenericAffExpr{Float64,JuMP.Variable},2,Array{JuMP.GenericAffExpr{Float64,JuMP.Variable},2}} with indices -2:0×-2:0:
 v[1]  0     0   
 0     v[2]  0   
 0     0     v[3]

Problems of course are that k == 0 is assumed, and that things could probably be improved by linearly indexing into A.

That brings up a different question: what should indices(A) be if k != 0?

[timholy]

That improved diagm would be a great contribution to Base, by all means please submit it.

That brings up a different question: what should indices(A) be if k != 0?

Seems like one might want an interface diagm(v, Δrow, Δcol). But for diagm(v, k), I would follow the rule that, if you snipped out a matrix that would be equal in value to diagm(v), that it should have the same indices. So for your example above, diagm(v, 0) would be a -2:0x-2:0 matrix, diagm(v, 1) would be a -2:1x-3:0 matrix, and diagm(v, -1) would be a -3:0x-2:1 matrix.

[tkoolen]

diagm(v, 1) would be a -2:1x-3:0 matrix, and diagm(v, -1) would be a -3:0x-2:1 matrix

Interesting. That does seem hard to make consistent with the current behavior for regular Arrays though, where map(first, indices(A)) will be (1, 1) no matter what k is. Also, I think one of the main applications of OffsetArrays is having all arrays you ever have to deal with be zero-based, so changing the start of the range would probably be quite surprising for users.

Another option is to change the last index for each dimension, so -2:1x-2:1 for both diagm(v, 1) and diagm(v, -1). That may be a little bit harder to implement generically in Base (but not that much harder).

[fredrikekre]

Cross-ref JuliaLang/julia#24047: I don't think we can call similar since we will now allow more than just one vector as input. Perhaps the solution to this issue is to define a method diagm(::Pair{<:OffsetVector,<:Integer}...)) in OffsetArrays.jl instead?

[jishnub]

The specific example listed in the OP seems to work now

julia> A = OffsetArray(rand(3), -3)
3-element OffsetArray(::Array{Float64,1}, -2:0) with eltype Float64 with indices -2:0:
 0.5009784688952947
 0.09987117852205074
 0.6784469095270629

julia> diagm(A)
3×3 Array{Float64,2}:
 0.500978  0.0        0.0
 0.0       0.0998712  0.0
 0.0       0.0        0.678447

Closing as old. Please reopen if further discussion on this is necessary.