AutoDiff?

[cscherrer]

Hi,

I'd like to use ApproxFun for Bayesian modeling. So I have this function:

function smoothness(f::Fun)
	s = f.space
	∫ = DefiniteIntegral(s)
	∂² = Derivative(s,2)
	-∫ * (∂² * f)^2
end

This seems to work well, for example I can do

s = Chebyshev()
ℓ(nt::NamedTuple) = smoothness(Fun(s,nt.x))

julia> ℓ((x=randn(20),))
-6.885525125826199e6 on ApproxFunBase.AnyDomain()

But most samplers work in terms of autodiff, usually ForwardDiff by default. And I get this error:

julia> chain = initialize!(DynamicHMCChain, ℓ, tr)
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:x,), Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity, 1}}}}, typeof(ℓ)}}, Float64}, Float64, 10})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:x,), Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity, 1}}}}, typeof(ℓ)}}, Float64}, Float64, 10})
    @ Base ./number.jl:7
  [2] convert(#unused#::Type{Float64}, f::Fun{ConstantSpace{ApproxFunBase.AnyDomain, Float64}, ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:x,), Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity, 1}}}}, typeof(ℓ)}}, Float64}, Float64, 10}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:x,), Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity, 1}}}}, typeof(ℓ)}}, Float64}, Float64, 10}}})
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/WeoCn/src/Spaces/ConstantSpace.jl:70
  [3] setindex!(A::Vector{Float64}, x::Function, i1::Int64)
    @ Base ./array.jl:839

Do you know of a way to make this work with autodiff? Could the Float64 type constraint possibly be removed, or would this cause trouble elsewhere?

[dlfivefifty]

your stacktrace is cut off so I have no idea what's calling setindex!. Likely the array type is wrong here.

[cscherrer]

Sorry, I didn't think earlier calls were relevant. I think I found the problem, though. I expected smoothness to return a float, but instead it returns a Fun that happens to be constant. For example, this works:

julia> using ApproxFun

julia> function smoothness(f::Fun)
               s = f.space
               ∫ = DefiniteIntegral(s)
               ∂² = Derivative(s,2)
               result = -∫ * (∂² * f)^2
               result.coefficients[1]
       end
smoothness (generic function with 1 method)

julia> s = Chebyshev()
Chebyshev()

julia> ℓ(x) = smoothness(Fun(s,x))
ℓ (generic function with 2 methods)

julia> using ForwardDiff

julia> ForwardDiff.gradient(ℓ, randn(10))
10-element Vector{Float64}:
      0.0
      0.0
    -37.724699060965875
  11184.992849921739
  -1902.2097186821084
  87554.55322301129
  -5927.948750101238
 268579.1170432344
 -21127.338515818305
 630881.0238039208

But result.coefficients[1] is probably very brittle. What's the right way to turn a constant Fun into a <: Real?

[dlfivefifty]

i think Number(result)