Cannot construct a parameters object

[rvignolo]

Hi!

First of all, thanks for the package, it is great.

Let me be as precise as posible. This would seem like a lot to process but it is actually quite simple to understand.

I am facing an interesting issue when defining parameters. As you may know, named tuples are nice because they are inmutable and type stable, so they can be used as a parameters container for a model (of course structs are also useful, but for this particular example I need to work with named tuples). However, if you try to do the following:

julia> using Parameters
julia> nt = (a = 1, f = t -> a + t)
(a = 1, f = var"#9#10"())
julia> nt.f(1)
ERROR: UndefVarError: a not defined
Stacktrace:
 [1] (::var"#9#10")(::Int64) at .\none:1
 [2] top-level scope at none:0

you get an error related to the scope of variable a.

To overcome this issue, you can use the Parameters.jl package, specifically by means of the following macro:

julia> nt2 = @with_kw (a = 1, f = t -> a + t)
##NamedTuple_kw#511 (generic function with 2 methods)
julia> nt2()
(a = 1, f = var"#388#390"{Int64}(1))
julia> nt2().f
#388 (generic function with 1 method)
julia> nt2().f(1)
2

Also, it can handle more complex examples:

julia> nt2 = @with_kw (a = 1, f = t -> a + t, g = t -> f(t) + sin(t))
##NamedTuple_kw#255 (generic function with 2 methods)

julia> nt2().g(1.)
2.8414709848078967

This is nice. What is happening in the background is that the @with_kw macro produces a constructor function that receives named arguments so you can use them to define other parameters. You can inspect this using @macroexpand or @expand (from MacroTools.jl).

Now I want to make use of these features in a more complicated example. Let's say I would like to include the numerical derivative of f in my parameters by using GSL.jl package (I have already tested other packages but this one seems to be the fastest at the moment). The first thing I have to do, is define a gsl_function object:

julia> using GSL

julia> const h = 9.765625e-4
0.0009765625

julia> const abserr = Cdouble[0]
const 1-element Array{Float64,1}:
 0.0

julia> const result = Cdouble[0]
1-element Array{Float64,1}:
 0.0

julia> nt2 = @with_kw (a = 1, f = t -> a + t, f′ = t -> deriv_central(@gsl_function(f), t, h, result, abserr))##NamedTuple_kw#262 (generic function with 2 methods)

julia> nt2()
(a = 1, f = var"#80#84"{Int64}(1), f′ = var"#81#85"())

julia> nt2().f′(0.1)
ERROR: UndefVarError: f not defined
Stacktrace:
 [1] (::var"#88#89")(::Float64, ::Ptr{Nothing}) at /Users/ramirovignolo/.julia/packages/GSL/yvp3l/src/manual_wrappers.jl:45
 [2] deriv_central(::gsl_function, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/ramirovignolo/.julia/packages/GSL/yvp3l/src/gen/direct_wrappers/gsl_deriv_h.jl:36
 [3] (::var"#81#85")(::Float64) at ./none:1
 [4] top-level scope at none:1

It seems that @gsl_function expands and uses @cfunction (from Base). When this macro is expanded, f is not captured. Does anyone knows how to avoid this problem at all? Thanks!!

[rvignolo]

I have also tried using a let block:

let a = 1, f = t -> a + t, f′ = t -> deriv_central(@gsl_function(f), t, h, result, abserr)
    (a= a, f = f, f′ = f′)
end

And also replace the @gsl_function macro call with a function call:

julia> nt = let a = 1, f = t -> a + t, f′ = t -> deriv_central(gsl_function(@cfunction( (x,p) -> f(x), Cdouble, (Cdouble, Ptr{Cvoid})), 0), t, h, result, abserr)
           (a= a, f = f, f′ = f′)
       end
(a = 1, f = var"#106#108"{Int64}(1), f′ = var"#107#109"())
julia> nt.f(1)
2
julia> nt.f′(1.)
ERROR: UndefVarError: f not defined
Stacktrace:
 [1] (::var"#110#111")(::Float64, ::Ptr{Nothing}) at ./none:1
 [2] deriv_central(::gsl_function, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/ramirovignolo/.julia/packages/GSL/yvp3l/src/gen/direct_wrappers/gsl_deriv_h.jl:36
 [3] (::var"#107#109")(::Float64) at ./none:1
 [4] top-level scope at none:1

without success.

Thanks in advance!

[ludvigak]

Hi!
I must admit I don't know enough about Julia macros to really explain what is going on, and looks like you're in pretty deep with the macros here. However, I know that the @cfunction macro can be called so that it uses a runtime closure instead of compile-time evaluation, since we're already using this in our @gsl_function macro, see

GSL.jl/src/manual_wrappers.jl

Line 43 in 0d53507

macro gsl_function(f)

Perhaps you could get somewhere by creating the gsl_function struct yourself instead of using our macro?

As a side note, this looks like a terribly complicated way of evaluating a central difference, but perhaps you're planning to something more complicated later on?

[rvignolo]

Hi @ludvigak,

First of all, thanks for your answer!

I must admit I don't know enough about Julia macros to really explain what is going on, and looks like you're in pretty deep with the macros here. However, I know that the @cfunction macro can be called so that it uses a runtime closure instead of compile-time evaluation, since we're already using this in our @gsl_function macro,

I believe that the following example may be useful for you because it does not use @with_kw and it just uses @gsl_function:

function testing()

  a = 1
  g = t -> a + t

  h = 9.765625e-4
  result = Cdouble[0]
  abserr = Cdouble[0]

  g′ = t -> deriv_forward(@gsl_function(g), t, h, result, abserr)

  g′(1.)

  return result[], abserr[]
end

Running this function returns:

julia> testing()
ERROR: UndefVarError: g not defined
Stacktrace:
 [1] (::var"#112#113")(::Float64, ::Ptr{Nothing}) at /Users/ramirovignolo/.julia/packages/GSL/yvp3l/src/manual_wrappers.jl:45
 [2] deriv_forward(::gsl_function, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/ramirovignolo/.julia/packages/GSL/yvp3l/src/gen/direct_wrappers/gsl_deriv_h.jl:93
 [3] #109 at /Users/ramirovignolo/codigos/finance/UniversalMonteCarlo/examples/interestrate/shortrate/multifactor/FongVasicek_revisited.jl:328 [inlined]
 [4] testing() at /Users/ramirovignolo/codigos/finance/UniversalMonteCarlo/examples/interestrate/shortrate/multifactor/FongVasicek_revisited.jl:330
 [5] top-level scope at none:1

So it seems that the runtime closure g is not captured.

Perhaps you could get somewhere by creating the gsl_function struct yourself instead of using our macro?

Yes, this is great advice and I have also done it. Let's apply it to the previous case:

function testing2()

  a = 1
  g = t -> a + t

  h = 9.765625e-4
  result = Cdouble[0]
  abserr = Cdouble[0]

  # g′ = t -> deriv_forward(@gsl_function(g), t, h, result, abserr)
  g′ = t -> deriv_forward(gsl_function(@cfunction((x, p) -> g(x), Cdouble, (Cdouble, Ptr{Cvoid})), 0), t, h, result, abserr)

  g′(1.)

  return result[], abserr[]
end

I get the same error as before when running it. However, maybe there are more things that we can try in this context.

As a side note, this looks like a terribly complicated way of evaluating a central difference, but perhaps you're planning to something more complicated later on?

I am interested in this point. Could you elaborate? For my part, this was a simple and silly example. Are you referring to the fact that the derivative of a function like g can be computed by other means, like automatic differentiation? Please, do not hesitate to share your thoughts about this.

Thanks and regards!

[giordano]

Looks like there is a missing esc

[rvignolo]

in @cfunction macro, right? For this particular case, this is its expansion:

$(Expr(:cfunction, Ptr{Nothing}, :($(QuoteNode(:((x, p)->begin
          g(x)
      end)))), :Cdouble, :(Core.svec(Cdouble, Ptr{Cvoid})), :(:ccall)))

[giordano]

No, in the definition of @gsl_function in the first line of your stacktrace

[rvignolo]

But you get the same error when you remove the @gsl_function macro...

[giordano]

julia> using GSL

julia> function testing2()
           a = 1
           g = t -> a + t

           h = 9.765625e-4
           result = Cdouble[0]
           abserr = Cdouble[0]

           # g′ = t -> deriv_forward(@gsl_function(g), t, h, result, abserr)
           ptr = @cfunction($((x, p) -> g(x)), Cdouble, (Cdouble, Ptr{Cvoid}))
           GC.@preserve ptr g′ = t -> deriv_forward(gsl_function(Base.unsafe_convert(Ptr{Cvoid}, ptr), 0), t, h, result, abserr)

           g′(1.)

           return result[], abserr[]
       end
testing2 (generic function with 1 method)

julia> testing2()
(0.9999999999999991, 7.544823210636764e-11)

[ludvigak]

Whoa, nice! Ugly, but nice ;)

[rvignolo]

I have no idea how you got to this point but this is amazing. Thanks for that.

[ludvigak]

Regarding my other comment, if it's just derivatives you want, I'm not sure it's worth using GSL for it. The algorithm is pretty simple, and there's also FiniteDifferences.jl (but I guess you've seen that already).

Or, you could just use the good old (f(x+h)-f(x-h))/(2h) with a suitably small h. Something like 1e-6 will take you a long way.

[giordano]

Whoa, nice! Ugly, but nice ;)

I have no idea how you got to this point but this is amazing. Thanks for that.

I've had my fair share of troubles in the past: https://discourse.julialang.org/t/using-a-base-cfunction-in-ccall/13621.

Should @gsl_function be changed accordingly?

I'm pretty sure all gsl_function does miss esc. All arguments of macros should be escaped for hygiene as appropriate

[rvignolo]

Regarding my other comment, if it's just derivatives you want, I'm not sure it's worth using GSL for it. The algorithm is pretty simple, and there's also FiniteDifferences.jl (but I guess you've seen that already).

Exactly, GSL just for derivatives. At the moment this is way faster than FinniteDifferences.jl. Also, the algorithm is "simple" but it is way better than just using a constant h. GSL modifies this value several times. Check it out here.

I've had my fair share of troubles in the past: https://discourse.julialang.org/t/using-a-base-cfunction-in-ccall/13621.

Thank you for sharing this!

[rvignolo]

@giordano and @ludvigak do you know how to improve this code:

using GSL

macro gsl_derivative(f, x, at)
  p = gensym(:p)
  return quote
    ptr = @cfunction($(Expr(:$, :(($(esc(x)), $p) -> $(esc(f))))), Cdouble, (Cdouble, Ptr{Cvoid}))
    GC.@preserve ptr fgsl = gsl_function(Base.unsafe_convert(Ptr{Cvoid}, ptr), 0)
    gsl_derivative(fgsl, $(esc(at)))
  end
end

const gsl_h = 9.765625e-4         # (1/2)^-10
const gsl_result = Cdouble[0]
const gsl_abserr = Cdouble[0]

function gsl_derivative(f::F, at::T) where {F<:gsl_function,T<:Number}
  deriv_central(f, at, gsl_h, gsl_result, gsl_abserr)
  return gsl_result[]
end

This allows to construct a parameters container such as:

using Parameters
nt = @with_kw (a = 1, g = t -> 2 * t^2 + a, g′ = t -> @gsl_derivative(g(x) + a * x + sin(x), x, t))

That can be instantiated and evaluated:

params = nt()
params.g′(1.)

but I think that each function call for g′ requires to create the ptr pointer and the gsl_function object, which is not desirable.