[julia-users] Re: General question on indexing
Seth,
thanks for that, and yes we can move over there, but I need some time to
study it before I contribute in the context. I've been looking at
GraphLayout - spring.jl, and thinking along the same lines about speed. I'm
delaying a commit because the package has a coherent coding style but which
I find hard to follow - and an interface problem with Compose. Also, I
don't realistically expect to make a real improvement except in my own
skills.
>From the documentation and small tests, I did not expect any improvements
using @simd as long as there are "if"-sentences in the middle of it. So
perhaps this format of the adjacency matrix may help? This for both
directed and non-directed graphs. The type definition AdjMatrix just a
small adaption of Maq7 above. Each node could be processed from one row and
a submatrix pointing to the positions of just the other nodes.
# using 'current node -neighbour node ' notation:
4x4 Array{Int64,2}:
0 21 31 41
12 0 32 42
13 23 0 43
14 24 34 0
julia> ama =AdjMatrix(am)
3x4 AdjMatrix:
12 21 31 41
13 23 32 42
14 24 34 43
[julia-users] Re: General question on indexing
Seth, I appreciate that. As a novice programmer I appreciate to have exchanges with you guys who are making this. I'm using commercial software every day that hasn't made progress since the nineties. I worked through something like fifteen types for the adjacency matrix (and, by the way, tuples for storage worked out, but was slower). I suppose Julia is fast so that we don't actually need to optimize like this all the time, but it's an excellent lesson. In the example I gave above you may have noted some confusion about rows-columns order. Those algorithms which were adapted from text books or other languages seem to walk along rows. In Julia, you walk faster downhill. And with @simd.
[julia-users] Re: General question on indexing
I wish I walked faster with @simd - for some reason I'm getting the same performance as described here: https://groups.google.com/d/msg/julia-users/DycY6jwDcWs/yYXHsWF9AwAJ One of the key objectives for LightGraphs is to be as fast as possible, so while Julia itself is pretty darn good, we want to eke every possible performance gain out of it so that we can scale to very large graphs. This also includes memory allocation, so fast and small. We're getting to the point where improving one results in bad performance in the other. That's both a good sign (that we've optimized away the low-hanging fruit) and a bad sign (that we may be hitting the limits of what we can do). I'm sure we've missed some stuff, though. Thanks for the discussion. Feel free to join us over at https://github.com/JuliaGraphs/LightGraphs.jl/issues if this sort of stuff interests you :) Seth. On Monday, November 16, 2015 at 11:28:15 AM UTC-8, hustf wrote: > > Seth, I appreciate that. As a novice programmer I appreciate to have > exchanges with you guys who are making this. I'm using commercial software > every day that hasn't made progress since the nineties. > > I worked through something like fifteen types for the adjacency matrix > (and, by the way, tuples for storage worked out, but was slower). I suppose > Julia is fast so that we don't actually need to optimize like this all the > time, but it's an excellent lesson. > > In the example I gave above you may have noted some confusion about > rows-columns order. Those algorithms which were adapted from text books or > other languages seem to walk along rows. In Julia, you walk faster > downhill. And with @simd. >
Re: [julia-users] Re: General question on indexing
This may also be a case where writing an iterator that just visits the
elements with i >= j might be what you want. You can study the implementation
of CartesianIndex in Base for inspiration.
--Tim
On Saturday, November 14, 2015 02:18:45 PM hustf wrote:
> I believe this work-in-progress type may be adapted? I wanted to make this
> faster by using tuples or immutable fixed vectors to store the data or
> whatnot, but I guess I'm currently stuck. Help wanted.
>
> Compared to looping unnecessarily over a full matrix, this speeds things up
> by around 35%.
>
> immutable Maq7 <:AbstractMatrix{Int}
> data::Array{Int64,2} # upper triangular except diagonal, stored as a
> horizontal vector
> width::Int64
> function Maq7(mat::Matrix{Int64})
> if issym(mat)
> n=size(mat,1)
> if n>1
> v=mat[2:n,1]
> for col = 2:(n-1)
> append!(v,mat[(col+1):n,col])
> end
> new(v', size(mat,1))
> else
> error("Ouch, not > (1x1)")
> end
> else
> error("Ouch, not square symmetric")
> end
> end
> end
> Base.size(A::Maq7)= A.width, A.width
> Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
> Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
> function Base.getindex(A::Maq7, r::Int , c::Int )
> if c > r
> return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
> elseif r > c
> return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
> else
> return zero(eltype(A.data))
> end
> end
> Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
>
>
>
> Example:
> julia> ms
> 4x4 Array{Int64,2}:
> 7 1 2 3
> 1 8 4 5
> 2 4 9 6
> 3 5 6 10
>
>
> julia> m=Maq7(ms)
> 4x4 Maq7:
> 0 1 2 3
> 1 0 4 5
> 2 4 0 6
> 3 5 6 0
>
> julia> collect(m)
> 6-element Array{Int64,1}:
> 1
> 2
> 3
> 4
> 5
> 6
>
>
> julia> length(m)
> 6
Re: [julia-users] Re: General question on indexing
Thank you, Tim - I'll definitely check it out. Right now I can't seem to
find an efficient way of doing it, but there are a few tricks in
Iterators.jl that might work.
On Sunday, November 15, 2015 at 1:33:29 PM UTC-8, Tim Holy wrote:
>
> This may also be a case where writing an iterator that just visits the
> elements with i >= j might be what you want. You can study the
> implementation
> of CartesianIndex in Base for inspiration.
>
> --Tim
>
> On Saturday, November 14, 2015 02:18:45 PM hustf wrote:
> > I believe this work-in-progress type may be adapted? I wanted to make
> this
> > faster by using tuples or immutable fixed vectors to store the data or
> > whatnot, but I guess I'm currently stuck. Help wanted.
> >
> > Compared to looping unnecessarily over a full matrix, this speeds things
> up
> > by around 35%.
> >
> > immutable Maq7 <:AbstractMatrix{Int}
> > data::Array{Int64,2} # upper triangular except diagonal, stored as a
> > horizontal vector
> > width::Int64
> > function Maq7(mat::Matrix{Int64})
> > if issym(mat)
> > n=size(mat,1)
> > if n>1
> > v=mat[2:n,1]
> > for col = 2:(n-1)
> > append!(v,mat[(col+1):n,col])
> > end
> > new(v', size(mat,1))
> > else
> > error("Ouch, not > (1x1)")
> > end
> > else
> > error("Ouch, not square symmetric")
> > end
> > end
> > end
> > Base.size(A::Maq7)= A.width, A.width
> > Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
> > Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
> > function Base.getindex(A::Maq7, r::Int , c::Int )
> > if c > r
> > return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
> > elseif r > c
> > return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
> > else
> > return zero(eltype(A.data))
> > end
> > end
> > Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
> >
> >
> >
> > Example:
> > julia> ms
> > 4x4 Array{Int64,2}:
> > 7 1 2 3
> > 1 8 4 5
> > 2 4 9 6
> > 3 5 6 10
> >
> >
> > julia> m=Maq7(ms)
> > 4x4 Maq7:
> > 0 1 2 3
> > 1 0 4 5
> > 2 4 0 6
> > 3 5 6 0
> >
> > julia> collect(m)
> > 6-element Array{Int64,1}:
> > 1
> > 2
> > 3
> > 4
> > 5
> > 6
> >
> >
> > julia> length(m)
> > 6
>
>
Re: [julia-users] Re: General question on indexing
Haven't tested timing carefully, but something like the following should be
pretty decent:
julia> immutable LowerTriangularIterator{A<:AbstractMatrix}
data::A
end
julia> Base.start(iter::LowerTriangularIterator) = (1,1)
start (generic function with 49 methods)
julia> Base.done(iter::LowerTriangularIterator, state) = state[2] >
size(iter.data, 2)
done (generic function with 51 methods)
julia> @inline function Base.next(iter::LowerTriangularIterator, state)
item = iter.data[state[1], state[2]]
newstate = ifelse(state[1] == size(iter.data, 1), (state[2]+1,
state[2]+1), (state[1]+1, state[2]))
item, newstate
end
To make it vectorize (i.e., SIMD), you might also need definitions like these:
https://github.com/JuliaLang/julia/blob/aa24ee0d0f334e8d9ecbced64bd975a6cc5f240d/base/multidimensional.jl#L134-L150
but perhaps vectorization isn't applicable to your problem anyway.
--Tim
On Sunday, November 15, 2015 05:06:23 PM Seth wrote:
> Thank you, Tim - I'll definitely check it out. Right now I can't seem to
> find an efficient way of doing it, but there are a few tricks in
> Iterators.jl that might work.
>
> On Sunday, November 15, 2015 at 1:33:29 PM UTC-8, Tim Holy wrote:
> > This may also be a case where writing an iterator that just visits the
> > elements with i >= j might be what you want. You can study the
> > implementation
> > of CartesianIndex in Base for inspiration.
> >
> > --Tim
> >
> > On Saturday, November 14, 2015 02:18:45 PM hustf wrote:
> > > I believe this work-in-progress type may be adapted? I wanted to make
> >
> > this
> >
> > > faster by using tuples or immutable fixed vectors to store the data or
> > > whatnot, but I guess I'm currently stuck. Help wanted.
> > >
> > > Compared to looping unnecessarily over a full matrix, this speeds things
> >
> > up
> >
> > > by around 35%.
> > >
> > > immutable Maq7 <:AbstractMatrix{Int}
> > > data::Array{Int64,2} # upper triangular except diagonal, stored as a
> > > horizontal vector
> > > width::Int64
> > > function Maq7(mat::Matrix{Int64})
> > > if issym(mat)
> > > n=size(mat,1)
> > > if n>1
> > > v=mat[2:n,1]
> > > for col = 2:(n-1)
> > > append!(v,mat[(col+1):n,col])
> > > end
> > > new(v', size(mat,1))
> > > else
> > > error("Ouch, not > (1x1)")
> > > end
> > > else
> > > error("Ouch, not square symmetric")
> > > end
> > > end
> > > end
> > > Base.size(A::Maq7)= A.width, A.width
> > > Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
> > > Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
> > > function Base.getindex(A::Maq7, r::Int , c::Int )
> > > if c > r
> > > return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
> > > elseif r > c
> > > return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
> > > else
> > > return zero(eltype(A.data))
> > > end
> > > end
> > > Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
> > >
> > >
> > >
> > > Example:
> > > julia> ms
> > >
> > > 4x4 Array{Int64,2}:
> > > 7 1 2 3
> > > 1 8 4 5
> > > 2 4 9 6
> > > 3 5 6 10
> > >
> > > julia> m=Maq7(ms)
> > >
> > > 4x4 Maq7:
> > > 0 1 2 3
> > > 1 0 4 5
> > > 2 4 0 6
> > > 3 5 6 0
> > >
> > > julia> collect(m)
> > >
> > > 6-element Array{Int64,1}:
> > > 1
> > > 2
> > > 3
> > > 4
> > > 5
> > > 6
> > >
> > > julia> length(m)
> > > 6
Re: [julia-users] Re: General question on indexing
That's very cool and close to what I need. Thanks. I'm using adjacency
lists (vectors of vectors) but I'm sure I can adapt this code.
On Sunday, November 15, 2015 at 6:12:19 PM UTC-8, Tim Holy wrote:
>
> Haven't tested timing carefully, but something like the following should
> be
> pretty decent:
>
> julia> immutable LowerTriangularIterator{A<:AbstractMatrix}
>data::A
>end
>
> julia> Base.start(iter::LowerTriangularIterator) = (1,1)
> start (generic function with 49 methods)
>
> julia> Base.done(iter::LowerTriangularIterator, state) = state[2] >
> size(iter.data, 2)
> done (generic function with 51 methods)
>
> julia> @inline function Base.next(iter::LowerTriangularIterator, state)
>item = iter.data[state[1], state[2]]
>newstate = ifelse(state[1] == size(iter.data, 1), (state[2]+1,
> state[2]+1), (state[1]+1, state[2]))
>item, newstate
>end
>
>
> To make it vectorize (i.e., SIMD), you might also need definitions like
> these:
>
> https://github.com/JuliaLang/julia/blob/aa24ee0d0f334e8d9ecbced64bd975a6cc5f240d/base/multidimensional.jl#L134-L150
>
>
> but perhaps vectorization isn't applicable to your problem anyway.
>
> --Tim
>
> On Sunday, November 15, 2015 05:06:23 PM Seth wrote:
> > Thank you, Tim - I'll definitely check it out. Right now I can't seem to
> > find an efficient way of doing it, but there are a few tricks in
> > Iterators.jl that might work.
> >
> > On Sunday, November 15, 2015 at 1:33:29 PM UTC-8, Tim Holy wrote:
> > > This may also be a case where writing an iterator that just visits the
> > > elements with i >= j might be what you want. You can study the
> > > implementation
> > > of CartesianIndex in Base for inspiration.
> > >
> > > --Tim
> > >
> > > On Saturday, November 14, 2015 02:18:45 PM hustf wrote:
> > > > I believe this work-in-progress type may be adapted? I wanted to
> make
> > >
> > > this
> > >
> > > > faster by using tuples or immutable fixed vectors to store the data
> or
> > > > whatnot, but I guess I'm currently stuck. Help wanted.
> > > >
> > > > Compared to looping unnecessarily over a full matrix, this speeds
> things
> > >
> > > up
> > >
> > > > by around 35%.
> > > >
> > > > immutable Maq7 <:AbstractMatrix{Int}
> > > > data::Array{Int64,2} # upper triangular except diagonal, stored as
> a
> > > > horizontal vector
> > > > width::Int64
> > > > function Maq7(mat::Matrix{Int64})
> > > > if issym(mat)
> > > > n=size(mat,1)
> > > > if n>1
> > > > v=mat[2:n,1]
> > > > for col = 2:(n-1)
> > > > append!(v,mat[(col+1):n,col])
> > > > end
> > > > new(v', size(mat,1))
> > > > else
> > > > error("Ouch, not > (1x1)")
> > > > end
> > > > else
> > > > error("Ouch, not square symmetric")
> > > > end
> > > > end
> > > > end
> > > > Base.size(A::Maq7)= A.width, A.width
> > > > Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
> > > > Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
> > > > function Base.getindex(A::Maq7, r::Int , c::Int )
> > > > if c > r
> > > > return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
> > > > elseif r > c
> > > > return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
> > > > else
> > > > return zero(eltype(A.data))
> > > > end
> > > > end
> > > > Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
> > > >
> > > >
> > > >
> > > > Example:
> > > > julia> ms
> > > >
> > > > 4x4 Array{Int64,2}:
> > > > 7 1 2 3
> > > > 1 8 4 5
> > > > 2 4 9 6
> > > > 3 5 6 10
> > > >
> > > > julia> m=Maq7(ms)
> > > >
> > > > 4x4 Maq7:
> > > > 0 1 2 3
> > > > 1 0 4 5
> > > > 2 4 0 6
> > > > 3 5 6 0
> > > >
> > > > julia> collect(m)
> > > >
> > > > 6-element Array{Int64,1}:
> > > > 1
> > > > 2
> > > > 3
> > > > 4
> > > > 5
> > > > 6
> > > >
> > > > julia> length(m)
> > > > 6
>
>
[julia-users] Re: General question on indexing
Hustf, I don't know why my reply didn't show up here yesterday, but I
wanted to thank you for the code. I'm looking at how to make that work.
Thanks for your response.
On Saturday, November 14, 2015 at 2:18:45 PM UTC-8, hustf wrote:
>
> I believe this work-in-progress type may be adapted? I wanted to make this
> faster by using tuples or immutable fixed vectors to store the data or
> whatnot, but I guess I'm currently stuck. Help wanted.
>
> Compared to looping unnecessarily over a full matrix, this speeds things
> up by around 35%.
>
> immutable Maq7 <:AbstractMatrix{Int}
> data::Array{Int64,2} # upper triangular except diagonal, stored as a
> horizontal vector
> width::Int64
> function Maq7(mat::Matrix{Int64})
> if issym(mat)
> n=size(mat,1)
> if n>1
> v=mat[2:n,1]
> for col = 2:(n-1)
> append!(v,mat[(col+1):n,col])
> end
> new(v', size(mat,1))
> else
> error("Ouch, not > (1x1)")
> end
> else
> error("Ouch, not square symmetric")
> end
> end
> end
> Base.size(A::Maq7)= A.width, A.width
> Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
> Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
> function Base.getindex(A::Maq7, r::Int , c::Int )
> if c > r
> return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
> elseif r > c
> return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
> else
> return zero(eltype(A.data))
> end
> end
> Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
>
>
>
> Example:
> julia> ms
> 4x4 Array{Int64,2}:
> 7 1 2 3
> 1 8 4 5
> 2 4 9 6
> 3 5 6 10
>
>
> julia> m=Maq7(ms)
> 4x4 Maq7:
> 0 1 2 3
> 1 0 4 5
> 2 4 0 6
> 3 5 6 0
>
> julia> collect(m)
> 6-element Array{Int64,1}:
> 1
> 2
> 3
> 4
> 5
> 6
>
>
> julia> length(m)
> 6
>
>
>
>
>
[julia-users] Re: General question on indexing
I believe this work-in-progress type may be adapted? I wanted to make this
faster by using tuples or immutable fixed vectors to store the data or
whatnot, but I guess I'm currently stuck. Help wanted.
Compared to looping unnecessarily over a full matrix, this speeds things up
by around 35%.
immutable Maq7 <:AbstractMatrix{Int}
data::Array{Int64,2} # upper triangular except diagonal, stored as a
horizontal vector
width::Int64
function Maq7(mat::Matrix{Int64})
if issym(mat)
n=size(mat,1)
if n>1
v=mat[2:n,1]
for col = 2:(n-1)
append!(v,mat[(col+1):n,col])
end
new(v', size(mat,1))
else
error("Ouch, not > (1x1)")
end
else
error("Ouch, not square symmetric")
end
end
end
Base.size(A::Maq7)= A.width, A.width
Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
function Base.getindex(A::Maq7, r::Int , c::Int )
if c > r
return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
elseif r > c
return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
else
return zero(eltype(A.data))
end
end
Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
Example:
julia> ms
4x4 Array{Int64,2}:
7 1 2 3
1 8 4 5
2 4 9 6
3 5 6 10
julia> m=Maq7(ms)
4x4 Maq7:
0 1 2 3
1 0 4 5
2 4 0 6
3 5 6 0
julia> collect(m)
6-element Array{Int64,1}:
1
2
3
4
5
6
julia> length(m)
6
