RQIso cov function is type-unstable

[maximerischard]

using GaussianProcesses
k2=RQIso(0.0,0.0,0.0)
cov(k2, 3.0)
@code_warntype cov(k2, 3.0)

produces

Variables:
  #self#::Base.#cov
  rq::GaussianProcesses.RQIso
  r::Float64

Body:
  begin 
      SSAValue(2) = (Core.getfield)(rq::GaussianProcesses.RQIso,:σ2)::Float64
      SSAValue(1) = (Base.box)(Base.Float64,(Base.add_float)(1.0,(Base.box)(Base.Float64,(Base.div_float)(r::Float64,(Base.box)(Base.Float64,(Base.mul_float)((Base.box)(Base.Float64,(Base.mul_float)(2.0,(Core.getfield)(rq::GaussianProcesses.RQIso,:α)::Float64)),(Core.getfield)(rq::GaussianProcesses.RQIso,:ℓ2)::Float64))))))
      SSAValue(0) = (Base.box)(Base.Float64,(Base.neg_float)((Core.getfield)(rq::GaussianProcesses.RQIso,:α)::Float64))
      return (SSAValue(2) * $(Expr(:invoke, LambdaInfo for ^(::Float64, ::Float64), :(Base.^), SSAValue(1), SSAValue(0))))::ANY
  end::ANY

I think this is actually a julia bug, but I'm not sure how to reproduce it...

[maximerischard]

For reference, here's the code for cov:

cov(rq::RQIso, r::Float64) = rq.σ2*(1.0+r/(2.0*rq.α*rq.ℓ2)).^(-rq.α)

Changing it to

cov(rq::RQIso, r::Float64) = rq.σ2*(1.0+r/(2.0*rq.α*rq.ℓ2))^(-rq.α)

(.^ -> ^) doesn't fix it.

However, copy-pasting the function definition into the REPL, for some bizarre reason, does fix it:

Variables:
  #self#::Base.#cov
  rq::GaussianProcesses.RQIso
  r::Float64

Body:
  begin
      SSAValue(1) = (Core.getfield)(rq::GaussianProcesses.RQIso,:σ2)::Float64
      SSAValue(0) = $(Expr(:invoke, LambdaInfo for ^(::Float64, ::Float64), :(Main.^), :((Base.box)(Base.Float64,(Base.add_float)(1.0,(Base.box)(Base.Float64,(Base.div_float)(r,(Base.box)(Base.Float64,(Base.mul_float)((Base.box)(Base.Float64,(Base.mul_float)(2.0,(Core.getfield)(rq,:α)::Float64)),(Core.getfield)(rq,:ℓ2)::Float64))))))), :((Base.box)(Base.Float64,(Base.neg_float)((Core.getfield)(rq,:α)::Float64)))))
      return (Base.box)(Base.Float64,(Base.mul_float)(SSAValue(1),SSAValue(0)))
  end::Float64

[maximerischard]

This problem seems to be caused by importing the Optim package. I've been able to boil down the problem down to just a few lines. I created a new package, FooBar, which is just this:

module FooBar

using Optim

f(a::Float64, b::Float64) = a^b

end # module

FooBar.f turns out to be type-unstable. Commenting out using Optim makes it type-stable. Very odd...

[maximerischard]

So… it seems we're bumping into JuliaNLSolvers/Optim.jl#283

[maximerischard]

This is fixed in the development branch of julia 0.5

using GaussianProcesses
k2=RQIso(0.0,0.0,0.0)
cov(k2, 3.0)
@code_warntype cov(k2, 3.0)

Variables:
  #self#::Base.#cov
  rq::GaussianProcesses.RQIso
  r::Float64

Body:
  begin 
      SSAValue(1) = (Core.getfield)(rq::GaussianProcesses.RQIso,:σ2)::Float64
      SSAValue(0) = $(Expr(:invoke, LambdaInfo for ^(::Float64, ::Float64), :(GaussianProcesses.^), :((Base.box)(Base.Float64,(Base.add_float)(1.0,(Base.box)(Base.Float64,(Base.div_float)(r,(Base.box)(Base.Float64,(Base.mul_float)((Base.box)(Base.Float64,(Base.mul_float)(2.0,(Core.getfield)(rq,:α)::Float64)),(Core.getfield)(rq,:ℓ2)::Float64))))))), :((Base.box)(Base.Float64,(Base.neg_float)((Core.getfield)(rq,:α)::Float64)))))
      return (Base.box)(Base.Float64,(Base.mul_float)(SSAValue(1),SSAValue(0)))
  end::Float64

this should make it to the release of julia 0.5 soon, at which point we should close this issue.

[chris-nemeth]

Agreed