What if you add @simd after @inbounds?
On Friday, September 16, 2016 at 5:48:00 PM UTC+2, Ben Ward wrote: > > The code for that function not shown is rather similar to the version that > is shown: > > function > distance{T<:MutationType,A<:NucleotideAlphabet}(::Type{Proportion{T}}, > seqs::Vector{BioSequence{A}}) > d, l = distance(Count{T}, seqs) > D = Vector{Float64}(length(d)) > @inbounds for i in 1:length(D) > D[i] = d[i] / l[i] > end > return D, l > end > > Which has a @code_warntype which seems ok: > > *julia> **@code_warntype distance(Proportion{AnyMutation}, dnas2)* > > Variables: > > #self#::Bio.Var.#distance > > #unused#::Type{Bio.Var.Proportion{Bio.Var.AnyMutation}} > > seqs@_3::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1} > > d::Array{Int64,1} > > l::Array{Int64,1} > > #temp#@_6::Int64 > > D::Array{Float64,1} > > #temp#@_8::Int64 > > i::Int64 > > seqs@_10::Array{Bio.Seq.DNANucleotide,2} > > > Body: > > begin > > $(Expr(:inbounds, false)) > > # meta: location /Users/bward/.julia/v0.5/Bio/src/var/distances.jl > distance 133 > > # meta: location > /Users/bward/.julia/v0.5/Bio/src/var/mutation_counting.jl count_mutations > 263 > > seqs@_10::Array{Bio.Seq.DNANucleotide,2} = $(Expr(:invoke, > LambdaInfo for > seqmatrix(::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}, > ::Symbol), :(Bio.Var.seqmatrix), :(seqs@_3), :(:seq))) > > # meta: pop location > > # meta: pop location > > $(Expr(:inbounds, :pop)) > > SSAValue(0) = $(Expr(:invoke, LambdaInfo for > count_mutations(::Type{Bio.Var.AnyMutation}, > ::Array{Bio.Seq.DNANucleotide,2}), :(Bio.Var.count_mutations), > Bio.Var.AnyMutation, :(seqs@_10))) > > #temp#@_6::Int64 = $(QuoteNode(1)) > > SSAValue(9) = (Base.getfield)(SSAValue(0),1)::Array{Int64,1} > > SSAValue(10) = (Base.box)(Int64,(Base.add_int)(1,1)) > > d::Array{Int64,1} = SSAValue(9) > > #temp#@_6::Int64 = SSAValue(10) > > SSAValue(11) = (Base.getfield)(SSAValue(0),2)::Array{Int64,1} > > SSAValue(12) = (Base.box)(Int64,(Base.add_int)(2,1)) > > l::Array{Int64,1} = SSAValue(11) > > #temp#@_6::Int64 = SSAValue(12) # line 256: > > SSAValue(6) = (Base.arraylen)(d::Array{Int64,1})::Int64 > > D::Array{Float64,1} = > (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Float64,1)::Type{Array{Float64,1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,Array{Float64,1},0,SSAValue(6),0)::Array{Float64,1} > > # line 257: > > $(Expr(:inbounds, true)) > > SSAValue(8) = (Base.arraylen)(D::Array{Float64,1})::Int64 > > SSAValue(13) = > (Base.select_value)((Base.sle_int)(1,SSAValue(8))::Bool,SSAValue(8),(Base.box)(Int64,(Base.sub_int)(1,1)))::Int64 > > #temp#@_8::Int64 = 1 > > 26: > > unless (Base.box)(Base.Bool,(Base.not_int)((#temp#@_8::Int64 === > (Base.box)(Int64,(Base.add_int)(SSAValue(13),1)))::Bool)) goto 37 > > SSAValue(14) = #temp#@_8::Int64 > > SSAValue(15) = (Base.box)(Int64,(Base.add_int)(#temp#@_8::Int64,1)) > > i::Int64 = SSAValue(14) > > #temp#@_8::Int64 = SSAValue(15) # line 258: > > SSAValue(5) = > (Base.box)(Base.Float64,(Base.div_float)((Base.box)(Float64,(Base.sitofp)(Float64,(Base.arrayref)(d::Array{Int64,1},i::Int64)::Int64)),(Base.box)(Float64,(Base.sitofp)(Float64,(Base.arrayref)(l::Array{Int64,1},i::Int64)::Int64)))) > > > (Base.arrayset)(D::Array{Float64,1},SSAValue(5),i::Int64)::Array{Float64,1} > > 35: > > goto 26 > > 37: > > $(Expr(:inbounds, :pop)) # line 260: > > return > (Core.tuple)(D::Array{Float64,1},l::Array{Int64,1})::Tuple{Array{Float64,1},Array{Int64,1}} > > end::Tuple{Array{Float64,1},Array{Int64,1}} > > But then again the loop in this function is also not vectorised which I > struggle to see why... unless division can't be. > > On Friday, September 16, 2016 at 3:27:47 AM UTC+1, Ralph Smith wrote: >> >> SSAValue(15) = (Base.getfield)(SSAValue(0),1) >> *::Union{Array{Float64,1},Array{Int64,1}}* >> >> >> indicates that the first element of SSAValue(0) is ambiguous. Earlier it >> shows that this means p from >> >> p, l = distance(Proportion{AnyMutation}, seqs) >> >> which we can't analyze from what you show here. >> >> On Thursday, September 15, 2016 at 10:08:16 AM UTC-4, Ben Ward wrote: >>> >>> Hi I have two functions and a function which calls them: >>> >>> @inline function expected_distance(::Type{JukesCantor69}, p::Float64) >>> return -0.75 * log(1 - 4 * p / 3) >>> end >>> >>> @inline function variance(::Type{JukesCantor69}, p::Float64, l::Int64) >>> return p * (1 - p) / (((1 - 4 * p / 3) ^ 2) * l) >>> end >>> >>> function distance{A<:NucleotideAlphabet}(::Type{JukesCantor69}, >>> seqs::Vector{BioSequence{A}}) >>> p, l = distance(Proportion{AnyMutation}, seqs) >>> D = Vector{Float64}(length(p)) >>> V = Vector{Float64}(length(p)) >>> @inbounds for i in 1:length(p) >>> D[i] = expected_distance(JukesCantor69, p[i]) >>> V[i] = variance(JukesCantor69, p[i], l[i]) >>> end >>> return D, V >>> end >>> >>> But I'm seeing type uncertainty: >>> >>> *@code_warntype distance(JukesCantor69, dnas)* >>> >>> Variables: >>> >>> #self#::Bio.Var.#distance >>> >>> #unused#::Type{Bio.Var.JukesCantor69} >>> >>> seqs::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1} >>> >>> p::Array{Float64,1} >>> >>> l::Array{Int64,1} >>> >>> #temp#@_6::Int64 >>> >>> D::Array{Float64,1} >>> >>> V::Array{Float64,1} >>> >>> #temp#@_9::Int64 >>> >>> i::Int64 >>> >>> >>> Body: >>> >>> begin >>> >>> SSAValue(0) = $(Expr(:invoke, LambdaInfo for >>> distance(::Type{Bio.Var.Proportion{Bio.Var.AnyMutation}}, >>> ::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}), >>> :(Bio.Var.distance), Bio.Var.Proportion{Bio.Var.AnyMutation}, :(seqs))) >>> >>> #temp#@_6::Int64 = $(QuoteNode(1)) >>> >>> SSAValue(15) = (Base.getfield)(SSAValue(0),1) >>> *::Union{Array{Float64,1},Array{Int64,1}}* >>> >>> SSAValue(16) = (Base.box)(Int64,(Base.add_int)(1,1)) >>> >>> p::Array{Float64,1} = SSAValue(15) >>> >>> #temp#@_6::Int64 = SSAValue(16) >>> >>> SSAValue(17) = (Base.getfield)(SSAValue(0),2) >>> *::Union{Array{Float64,1},Array{Int64,1}}* >>> >>> SSAValue(18) = (Base.box)(Int64,(Base.add_int)(2,1)) >>> >>> l::Array{Int64,1} = SSAValue(17) >>> >>> #temp#@_6::Int64 = SSAValue(18) # line 314: >>> >>> SSAValue(7) = (Base.arraylen)(p::Array{Float64,1})::Int64 >>> >>> D::Array{Float64,1} = >>> (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Float64,1)::Type{Array{Float64,1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,Array{Float64,1},0,SSAValue(7),0)::Array{Float64,1} >>> >>> # line 315: >>> >>> SSAValue(9) = (Base.arraylen)(p::Array{Float64,1})::Int64 >>> >>> V::Array{Float64,1} = >>> (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Float64,1)::Type{Array{Float64,1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,Array{Float64,1},0,SSAValue(9),0)::Array{Float64,1} >>> >>> # line 316: >>> >>> $(Expr(:inbounds, true)) >>> >>> SSAValue(11) = (Base.arraylen)(p::Array{Float64,1})::Int64 >>> >>> SSAValue(19) = >>> (Base.select_value)((Base.sle_int)(1,SSAValue(11))::Bool,SSAValue(11),(Base.box)(Int64,(Base.sub_int)(1,1)))::Int64 >>> >>> #temp#@_9::Int64 = 1 >>> >>> 22: >>> >>> unless (Base.box)(Base.Bool,(Base.not_int)((#temp#@_9::Int64 === >>> (Base.box)(Int64,(Base.add_int)(SSAValue(19),1)))::Bool)) goto 43 >>> >>> SSAValue(20) = #temp#@_9::Int64 >>> >>> SSAValue(21) = (Base.box)(Int64,(Base.add_int)(#temp#@_9::Int64,1)) >>> >>> i::Int64 = SSAValue(20) >>> >>> #temp#@_9::Int64 = SSAValue(21) # line 317: >>> >>> SSAValue(12) = >>> (Base.arrayref)(p::Array{Float64,1},i::Int64)::Float64 >>> >>> $(Expr(:inbounds, false)) >>> >>> # meta: location /Users/bward/.julia/v0.5/Bio/src/var/distances.jl >>> expected_distance 69 >>> >>> SSAValue(13) = $(Expr(:invoke, LambdaInfo for log(::Float64), >>> :(Bio.Var.log), >>> :((Base.box)(Base.Float64,(Base.sub_float)((Base.box)(Float64,(Base.sitofp)(Float64,1)),(Base.box)(Base.Float64,(Base.div_float)((Base.box)(Base.Float64,(Base.mul_float)((Base.box)(Float64,(Base.sitofp)(Float64,4)),SSAValue(12))),(Base.box)(Float64,(Base.sitofp)(Float64,3))))))))) >>> >>> # meta: pop location >>> >>> $(Expr(:inbounds, :pop)) >>> >>> SSAValue(5) = >>> (Base.box)(Base.Float64,(Base.mul_float)(-0.75,SSAValue(13))) >>> >>> >>> (Base.arrayset)(D::Array{Float64,1},SSAValue(5),i::Int64)::Array{Float64,1} >>> # line 318: >>> >>> SSAValue(14) = >>> (Base.arrayref)(p::Array{Float64,1},i::Int64)::Float64 >>> >>> SSAValue(6) = >>> (Base.box)(Base.Float64,(Base.div_float)((Base.box)(Base.Float64,(Base.mul_float)(SSAValue(14),(Base.box)(Base.Float64,(Base.sub_float)((Base.box)(Float64,(Base.sitofp)(Float64,1)),SSAValue(14))))),(Base.box)(Base.Float64,(Base.mul_float)((Base.Math.box)(Base.Math.Float64,(Base.Math.powi_llvm)((Base.box)(Base.Float64,(Base.sub_float)((Base.box)(Float64,(Base.sitofp)(Float64,1)),(Base.box)(Base.Float64,(Base.div_float)((Base.box)(Base.Float64,(Base.mul_float)((Base.box)(Float64,(Base.sitofp)(Float64,4)),SSAValue(14))),(Base.box)(Float64,(Base.sitofp)(Float64,3)))))),(Base.box)(Int32,(Base.checked_trunc_sint)(Int32,2))))::Float64,(Base.box)(Float64,(Base.sitofp)(Float64,(Base.arrayref)(l::Array{Int64,1},i::Int64)::Int64)))))) >>> >>> >>> (Base.arrayset)(V::Array{Float64,1},SSAValue(6),i::Int64)::Array{Float64,1} >>> >>> 41: >>> >>> goto 22 >>> >>> 43: >>> >>> $(Expr(:inbounds, :pop)) # line 320: >>> >>> return >>> (Core.tuple)(D::Array{Float64,1},V::Array{Float64,1})::Tuple{Array{Float64,1},Array{Float64,1}} >>> >>> end::Tuple{Array{Float64,1},Array{Float64,1}} >>> >>> But I'm not sure which those lines correspond to in my code, as they're >>> temporary values. I think at some point some code either results in an >>> integer or a float. I wondered if it was inside the smaller function called >>> by the larger one. >>> >>> Thanks, >>> Ben. >>> >>