将 Julia 回调函数传递给 C

Julia 的一大优势是，它可以非常轻松地调用 C 代码，无需任何特殊的“胶水”例程或开销来编组参数和转换返回值。例如，如果你想调用 GNU GSL 来计算一个特殊的函数，例如德拜积分，它就像这样简单

debye_1(x) = ccall((:gsl_sf_debye_1,:libgsl), Cdouble, (Cdouble,), x)

此时，你可以计算 debye_1(2)、debye_1(3.7) 等。(更简单的方法是使用 Jiahao Chen 的 GSL 包，它已经为你创建了这样的包装器。) 这使得大量现有的 C 库在 Julia 中可用 (以及 Fortran 库和其他具有 C 可访问调用约定的语言)。

事实上，你甚至可以反过来，将 Julia 例程传递给 C，以便 C 代码以回调函数的形式调用 Julia 代码。例如，一个用于数值积分的 C 库可能会期望你将被积函数作为函数参数传递，该库将调用该函数来尽可能多次地评估被积函数以估计积分。回调函数在优化、求根以及许多其他数值任务中也很自然，在许多非数值问题中也是如此。这篇博客文章的目的是说明将 Julia 函数作为回调传递给 C 例程的技术，这很直观且高效，但需要对函数和其他值作为参数传递方式的低级理解。

本文中的代码需要 Julia 0.2 (或其最新的 git 仿制品)；回调函数所需的关键功能 (尤其是 unsafe_pointer_to_objref) 在 Julia 0.1 中不可用。

使用 `qsort` 进行排序

也许最著名的回调参数示例是由 qsort 函数提供的，它是 ANSI C 标准库的一部分，在 C 中声明为

void qsort(void *base, size_t nmemb, size_t size,
           int(*compare)(const void *a, const void *b));

base 参数是指向长度为 nmemb 的数组的指针，每个元素的大小为 size 字节。compare 是一个回调函数，它接收指向两个元素 a 和 b 的指针，如果 a 应该出现在 b 之前，则返回一个小于零的整数；如果 a 应该出现在 b 之后，则返回一个大于零的整数 (如果任何顺序都允许，则返回零)。现在，假设我们在 Julia 中有一个 1 维数组 A，我们想使用 qsort 函数对它进行排序 (而不是使用 Julia 的内置 sort 函数)。在我们担心调用 qsort 和传递参数之前，我们需要编写一个适用于任意类型 T 的比较函数，例如

function mycompare{T}(a_::Ptr{T}, b_::Ptr{T})
    a = unsafe_load(a_)
    b = unsafe_load(b_)
    return a < b ? cint(-1) : a > b ? cint(+1) : cint(0)
end
cint(n) = convert(Cint, n)

请注意，我们使用内置函数 unsafe_load 来获取参数 a_ 和 b_ 指向的值 (这是“不安全的”，因为它将在这些不是有效指针时崩溃，但 qsort 将始终传递有效指针)。此外，我们必须对返回值小心一点：qsort 期望一个返回 C int 的函数，因此我们必须确保通过调用 convert 返回 Cint (Julia 中的对应类型)。

现在，我们如何将它传递给 C？C 中的函数指针本质上只是指向实现该函数的机器代码内存位置的指针，而 Julia 中的函数值 mycompare (类型为 Function) 则完全不同。由于 Julia 的 JIT 编译方法，Julia 函数可能甚至直到第一次调用时才会编译，一般来说，相同的 Julia 函数可能会被编译成多个机器代码实例，这些实例针对不同类型的参数进行了专门化 (例如，在这种情况下不同的 T)。因此，你可以想象 mycompare 必须在内部指向一个相当复杂的数据结构 (一个 jl_function_t，如果你感兴趣的话，在 julia.h 中)，它包含有关参数类型、已编译版本 (如果有) 等的信息。一般来说，它必须存储一个闭包，其中包含有关定义函数的环境的信息；我们将在下面进一步讨论这一点。无论如何，它是一个与指向一组参数类型机器代码的简单指针完全不同的对象。幸运的是，我们可以简单地通过调用一个内置的 Julia 函数来获取后者，该函数称为 cfunction

const mycompare_c = cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}))

在这里，我们将三个参数传递给 cfunction：函数 mycompare、返回类型 Cint 以及参数类型的元组，在这种情况下，是对 Cdouble (Float64) 元素数组进行排序。Julia 编译了一个针对这些参数类型专门化的 mycompare 版本 (如果它还没有这样做)，并返回一个 Ptr{Void}，它保存机器代码的地址，这正是我们需要传递给 qsort 的。现在我们已准备好对一些示例数据调用 qsort

A = [1.3, -2.7, 4.4, 3.1]
ccall(:qsort, Void, (Ptr{Cdouble}, Csize_t, Csize_t, Ptr{Void}),
      A, length(A), sizeof(eltype(A)), mycompare_c)

执行完此操作后，A 将更改为已排序的数组 [ -2.7, 1.3, 3.1, 4.4]。请注意，Julia 知道如何将数组 A::Vector{Cdouble} 转换为 Ptr{Cdouble}，如何计算类型的大小 (以字节为单位) (与 C 的 sizeof 运算符相同) 等。为了好玩，尝试在 mycompare 中插入一行 println("mycompare($a,$b)")，它将允许你看到 qsort 正在执行的比较 (并验证它是否真的在调用你传递给它的 Julia 函数)。

闭包问题

然而，我们还没有完成。如果你开始将回调函数传递给 C 例程，那么你很快就会发现 cfunction 并不总是有效。例如，假设我们尝试通过以下方式内联声明我们的比较函数：

mycomp = cfunction((a_,b_) -> unsafe_load(a_) < unsafe_load(b_) ?
                              cint(-1) : cint(+1),
                   Cint, (Ptr{Cdouble}, Ptr{Cdouble}))

Julia 会对此抛出异常，并打印 ERROR: function is not yet c-callable。一般来说，cfunction 仅适用于“顶级”函数：在顶级 (全局或模块) 范围内定义的命名函数，但不包括匿名 (args -> value) 函数，也不包括在其他函数内定义的函数 (“嵌套”函数)。造成这种情况的原因源于计算机科学中的一个重要概念：闭包。

为了理解对闭包的需求以及它们对回调函数带来的困难，假设我们想为 qsort 提供一个更好的接口，该接口允许用户只传递一个返回 true 或 false 的 lessthan 函数，同时隐藏与指针、Cint 等相关的所有底层操作。我们可能想做一些类似于以下形式的操作

function qsort!{T}(A::Vector{T}, lessthan::Function)
    function mycompare(a_::Ptr{T}, b_::Ptr{T})
        a = unsafe_load(a_)
        b = unsafe_load(b_)
        return lessthan(a, b) ? cint(-1) : cint(+1)
    end
    mycompare_c = cfunction(mycompare, Cint, (Ptr{T}, Ptr{T}))
    ccall(:qsort, Void, (Ptr{T}, Csize_t, Csize_t, Ptr{Void}),
          A, length(A), sizeof(T), mycompare_c)
    A
end

然后，我们可以简单地调用 qsort!([1.3, -2.7, 4.4, 3.1], <) 来使用内置的 < 比较按升序进行排序，或者使用我们想要的任何其他比较函数。不幸的是，当你尝试调用 qsort! 时，cfunction 会再次抛出异常，现在我们不难理解为什么。请注意，嵌套的 mycompare 函数不再是自包含的：它使用了来自周围作用域的变量 lessthan。这是嵌套函数和匿名函数的常见模式：它们通常由定义函数的环境中的局部变量参数化。从技术上讲，这种依赖关系的能力是由编程语言 (如 Julia) 中的词法作用域提供的，并且是任何将函数作为 "一等" 对象的语言的典型特征。为了支持词法作用域，Julia Function 对象需要在内部携带一个指向封闭环境中变量的指针，这种封装称为闭包。

相反，C 函数指针不是闭包。它不会将指针封闭到定义函数的环境中，或者任何其他东西；它只是指向一串指令的地址。这使得在 C 中很难编写函数来转换其他函数 (高阶函数) 或通过局部变量参数化函数。这似乎给我们留下了两种选择，这两种选择都不是特别吸引人

我们可以将 lessthan 存储在一个全局变量中，并从顶级的 mycompare 函数中引用它。(这是 C 程序员使用参数化比较函数调用 qsort 的传统解决方案。) 此策略的问题是它不是可重入的：它阻止我们递归调用 qsort! (例如，如果比较函数本身需要对某些复杂的数据结构进行排序)，或者从多个线程调用 qsort! (当未来的 Julia 版本支持共享内存并行时)。尽管如此，它总比没有好。
每次调用 qsort! 时，Julia 都会 JIT 编译 mycompare 的一个新版本，该版本会对该调用传递的 lessthan 参数的引用进行硬编码。从技术上讲，这是可能的，并且在某些语言中已经实现了 (例如，据报道 GNU Guile 和 Lua 确实做了类似的事情)。但是，这种策略是有代价的：它要求每次参数发生变化时都重新编译回调，而全局变量策略则没有这种要求。无论如何，它还没有在 Julia 中实现。

幸运的是，通常还有第三种选择，因为 C 程序员很久以前就认识到了函数指针的这些局限性，并想出了一种解决方法：大多数现代 C 回调接口允许通过一个“传递”(或“thunk”) 指针参数将任意数据传递给回调。正如下一节所解释的那样，我们可以利用这种技术在 Julia 中传递一个“真正”的闭包作为回调。

通过传递指针传递闭包

qsort 接口如今被认为是相当过时的。几年前，它在 BSD-Unix 系统上得到了补充，并在 GNU libc 中最终添加了一个名为 qsort_r 的函数，该函数解决了以可重入的方式将参数传递给回调的问题。这是 BSD (例如 MacOS) qsort_r 函数的定义

void qsort_r(void *base, size_t nmemb, size_t size, void *thunk,
             int (*compare)(void *thunk, const void *a, const void *b));

与 qsort 相比，多了一个 thunk 参数，并且该参数会传递给 compare 函数作为其第一个参数。这样，你可以将一个指向任意数据的指针传递给你的回调，我们可以利用它将闭包传递给任意 Julia 回调。

我们只需要一种方法将 Julia 的 `Function` 转换为不透明的 `Ptr{Void}`，这样我们就可以将其传递给回调函数，然后将不透明指针转换回 `Function`。前者如果我们简单地将 `ccall` 参数声明为类型 `Any`（这将以不透明的 Julia 对象指针传递参数）会自动完成，而后者由内置函数 `unsafe_pointer_to_objref` 完成。（从技术上讲，我们可以使用类型 `Function` 或显式调用 `pointer_from_objref` 而不是 `Any`。）使用这些，我们现在可以定义一个工作的高级 `qsort!` 函数，它接受一个任意的 `lessthan` 比较函数参数

function qsort!_compare{T}(lessthan_::Ptr{Void}, a_::Ptr{T}, b_::Ptr{T})
    a = unsafe_load(a_)
    b = unsafe_load(b_)
    lessthan = unsafe_pointer_to_objref(lessthan_)::Function
    return lessthan(a, b) ? cint(-1) : cint(+1)
end

function qsort!{T}(A::Vector{T}, lessthan::Function=<)
    compare_c = cfunction(qsort!_compare, Cint, (Ptr{Void}, Ptr{T}, Ptr{T}))
    ccall(:qsort_r, Void, (Ptr{T}, Csize_t, Csize_t, Any, Ptr{Void}),
          A, length(A), sizeof(T), lessthan, compare_c)
    return A
end

qsort!_compare 是一个顶级函数，因此 `cfunction` 对它没有问题，并且它将只为要排序的类型 `T` 编译一次（而不是为每次调用 `qsort!` 或每次 `lessthan` 函数编译一次）。我们使用显式的 `::Function` 断言来告诉编译器我们只会在 `lessthan_` 中传递 `Function` 指针。请注意，我们给 `lessthan` 参数设置了一个默认值为 `<`（默认参数是 Julia 中最近添加的功能）。

现在我们可以执行 `qsort!([1.3, -2.7, 4.4, 3.1])`，它将返回按升序排序的数组，或者 `qsort!([1.3, -2.7, 4.4, 3.1], >)` 来按降序排序。

警告：`qsort_r` 不便携

上面的例子有一个与 Julia 无关的主要问题：qsort_r 函数不具有可移植性。上面的例子在 Windows 上无法工作，因为 Windows C 库没有定义 `qsort_r`（相反，它有一个名为 qsort_s 的函数，它当然使用与 BSD 和 GNU qsort_r 函数都不兼容 的参数顺序）。更糟糕的是，它会在 GNU/Linux 系统上崩溃，这些系统确实提供了 qsort_r，但具有不兼容的调用约定。因此，很难以一种既不会在 GNU/Linux 也不会在 BSD（例如 MacOS）系统上崩溃的方式使用 `qsort_r`。这就是 glibc 的 `qsort_r` 的定义方式

void qsort_r(void *base, size_t nmemb, size_t size,
             int (*compare)(const void *a, const void *b, void *thunk),
              void *thunk);

请注意，`thunk` 参数的位置已更改，无论是在 `qsort_r` 本身还是在比较函数中。因此，GNU/Linux 系统上的相应 `qsort!` Julia 代码应该是

function qsort!_compare{T}(a_::Ptr{T}, b_::Ptr{T}, lessthan_::Ptr{Void})
    a = unsafe_load(a_)
    b = unsafe_load(b_)
    lessthan = unsafe_pointer_to_objref(lessthan_)::Function
    return lessthan(a, b) ? cint(-1) : cint(+1)
end

function qsort!{T}(A::Vector{T}, lessthan::Function=<)
    compare_c = cfunction(qsort!_compare, Cint, (Ptr{T}, Ptr{T}, Ptr{Void}))
    ccall(:qsort_r, Void, (Ptr{T}, Csize_t, Csize_t, Ptr{Void}, Any),
          A, length(A), sizeof(T), compare_c, lessthan)
    return A
end

如果您确实需要从 Julia 调用 `qsort_r`，您可以使用上面的定义，如果 `OS_NAME == :Linux`，则使用 BSD 定义，否则，在 Windows 上使用第三个版本使用 `qsort_s`，但幸运的是，Julia 自带了它自己的完全合适的 `sort` 和 `sort!` 例程，因此没有太多必要。

在数据结构中传递闭包

作为另一个更面向数值计算的例子，我们将研究如何调用 GNU Scientific Library (GSL) 中的数值积分例程。已经有一个 GSL 包可以为您处理下面的包装器工作，但是查看它是如何实现的很有帮助，因为 GSL 以稍微不同的方式模拟闭包，使用数据结构。

像大多数现代接受回调的 C 库一样，GSL 使用 `void*` 传递参数，以允许将任意数据传递到回调例程中，我们可以用它来支持 Julia 中的任意闭包。但是，与 `qsort_r` 不同，GSL 将 C 函数指针和传递指针都包装在一个名为 `gsl_function` 的数据结构中

struct {
    double (*function)(double x, void *params);
    void *params;
} gsl_function;

使用上面的技术，我们可以轻松地在 Julia 中声明一个 `GSL_Function` 类型，它与这种 C 类型相匹配，并且带有一个构造函数 `GSL_Function(f::Function)`，它在任意 Julia 函数 `f` 周围创建一个包装器

function gsl_function_wrap(x::Cdouble, params::Ptr{Void})
    f = unsafe_pointer_to_objref(params)::Function
    convert(Cdouble, f(x))::Cdouble
end
const gsl_function_wrap_c = cfunction(gsl_function_wrap,
                                      Cdouble, (Cdouble, Ptr{Void}))

type GSL_Function
    func::Ptr{Void}
    params::Any
    GSL_Function(f::Function) = new(gsl_function_wrap_c, f)
end

上面代码的一个细微之处是，我们需要显式地将 `f` 的返回值转换为 `Cdouble`（以防调用者的代码为某些 `x` 返回其他数值类型，例如 `Int`）。此外，我们需要显式断言（`::Cdouble`） `convert` 的结果是 `Cdouble`。与 `qsort` 示例一样，这是因为 `cfunction` 只有在 Julia 可以保证 `gsl_function_wrap` 返回指定的 `Cdouble` 类型时才起作用，并且 Julia 无法推断 `convert` 的返回类型，因为它不知道 `f(x)` 的返回类型。

有了上面的定义，将它传递给 GSL 自适应积分例程中的包装函数 `gsl_integration_qag` 只是一个简单的事情

function gsl_integration_qag(f::Function, a::Real, b::Real, epsrel::Real=1e-12,
                             maxintervals::Integer=10^7)
    s = ccall((:gsl_integration_workspace_alloc,:libgsl), Ptr{Void}, (Csize_t,),
              maxintervals)
    result = Array(Cdouble,1)
    abserr = Array(Cdouble,1)
    ccall((:gsl_integration_qag,:libgsl), Cint,
          (Ptr{GSL_Function}, Cdouble,Cdouble, Cdouble, Csize_t, Cint, Ptr{Void},
           Ptr{Cdouble}, Ptr{Cdouble}),
          &GSL_Function(f), a, b, epsrel, maxintervals, 1, s, result, abserr)
    ccall((:gsl_integration_workspace_free,:libgsl), Void, (Ptr{Void},), s)
    return (result[1], abserr[1])
end

请注意，`&GSL_Function(f)` 将指针传递给一个包含指向 `gsl_function_wrap_c` 和 `f` 的指针的 `GSL_Function` “结构体”，对应于 C 中的 `gsl_function*` 参数。返回值是估计积分和估计误差的元组。

例如，`gsl_integration_qag(cos, 0, 1)` 返回 `(0.8414709848078965,9.34220461887732e-15)`，它以机器精度计算出正确的积分 `sin(1)`。

清理垃圾（或不清理）

在上面的例子中，我们将一个 Julia `Function` 的不透明指针（对象引用）传递给 C。无论何时将指向 Julia 数据的指针传递给 C 代码，都必须确保 Julia 数据不会在 C 代码完成之前被垃圾回收，函数也不例外。不再被任何 Julia 变量引用的匿名函数可能会被垃圾回收，此时指向它的任何 C 指针都会变为无效。

这听起来很可怕，但实际上您不必经常担心，因为 Julia 保证 `ccall` 参数在 `ccall` 退出之前不会被垃圾回收。因此，在所有上面的例子中，我们都是安全的：`Function` 只需要在 `ccall` 的生命周期内存在即可。

唯一的危险出现在您将函数指针传递给 C 并且 C 代码保存指针 在某些数据结构中，它将在以后的 `ccall` 中使用它时。在这种情况下，您有责任确保 `Function` 变量存活（被某些 Julia 变量引用）只要 C 代码可能需要它。

例如，在 GSL 的一维最小化接口中，您不是简单地将目标函数传递给最小化例程并等待它被最小化。相反，您调用一个 GSL 例程来创建一个“最小化器对象”，将您的函数指针存储在这个对象中，调用例程来迭代最小化，然后在您完成时释放最小化器。Julia 函数在该过程完成之前不得被垃圾回收。确保这一点的最简单方法是在最小化器对象周围创建一个 Julia 包装器类型，该类型存储对 Julia 函数的显式引用，如下所示

type GSL_Minimizer
    m::Ptr{Void} # the gsl_min_fminimizer pointer
    f::Any  # explicit reference to objective, to prevent garbage-collection
    function GSL_Minimizer(t)
       m = ccall((:gsl_min_fminimizer_alloc,:libgsl), Ptr{Void}, (Ptr{Void},), t)
       p = new(m, nothing)
       finalizer(p, p -> ccall((:gsl_min_fminimizer_free,:libgsl),
                               Void, (Ptr{Void},), p.m))
       p
    end
end

它包装一个类型为 `t` 的 `gsl_min_fminimizer` 对象，带有一个占位符 `f` 来存储对目标函数的引用（一旦它在下面设置），包括一个 `finalizer` 来释放 `GSL_Minimizer` 被垃圾回收时 GSL 对象。参数 `t` 用于指定最小化算法，它可以默认通过以下方式使用 Brent 算法：

const gsl_brent = unsafe_load(cglobal((:gsl_min_fminimizer_brent,:libgsl), Ptr{Void}))
GSL_Minimizer() = GSL_Minimizer(gsl_brent)

（对 `cglobal` 的调用会产生指向 GSL 中的 `gsl_min_fminimizer_brent` 全局变量的指针，然后我们通过 `unsafe_load` 对其进行解引用以获取实际指针。）

然后，当我们设置要最小化的函数（“目标”）时，我们在 `GSL_Minimizer` 中存储对它的额外引用，以防止 `GSL_Minimizer` 的生命周期内发生垃圾回收，同样使用上面定义的 `GSL_Function` 类型来包装回调

function gsl_minimizer_set!(m::GSL_Minimizer, f, x0, xmin, xmax)
    ccall((:gsl_min_fminimizer_set,:libgsl), Cint,
          (Ptr{Void}, Ptr{GSL_Function}, Cdouble, Cdouble, Cdouble),
          m.m, &GSL_Function(f), x0, xmin, xmax)
    m.f = f
    m
end

然后，有一些 GSL 例程可以迭代最小化器并检查当前的 `x`、目标值或最小值的边界，这些例程便于包装

gsl_minimizer_iterate!(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_iterate,:libgsl), Cint, (Ptr{Void},), m.m)

gsl_minimizer_x(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_x_minimum,:libgsl), Cdouble, (Ptr{Void},), m.m)

gsl_minimizer_f(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_f_minimum,:libgsl), Cdouble, (Ptr{Void},), m.m)

gsl_minimizer_xmin(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_x_lower,:libgsl), Cdouble, (Ptr{Void},), m.m)
gsl_minimizer_xmax(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_x_upper,:libgsl), Cdouble, (Ptr{Void},), m.m)

将所有这些放在一起，我们可以通过以下方式在区间 [-3,1] 中最小化一个简单函数 `sin(x)`，初始猜测为 -1：

m = GSL_Minimizer()
gsl_minimizer_set!(m, sin, -1, -3, 1)
while gsl_minimizer_xmax(m) - gsl_minimizer_xmin(m) > 1e-6
    println("iterating at x = $(gsl_minimizer_x(m))")
    gsl_minimizer_iterate!(m)
end
println("found minimum $(gsl_minimizer_f(m)) at x = $(gsl_minimizer_x(m))")

经过几次迭代后，它会打印 `found minimum -1.0 at x = -1.5707963269964016`，这是大约 10 位数的正确最小值（−π/2）。

此时，我将厚颜无耻地推广一下我自己的 Julia NLopt 包，它包装了我的免费/开源 NLopt 库，以提供比 GSL 多得多的优化算法，并且可能具有更好的接口。但是，传递给 NLopt 的回调函数所使用的技术实际上与 GSL 所使用的技术非常相似。

在 PyCall 包中，可以找到这些技术的更复杂版本，用于从 Julia 调用 Python。为了将 Julia 函数传递给 Python，我们再次在处理类型转换等的包装函数上使用 `cfunction`，并通过传递指针传递实际的 Julia 闭包。但在这种情况下，传递指针包含一个 Python 对象，该对象已使用允许它包装 Julia 对象的新类型创建，并且通过将 Julia 对象存储在一个保存对象的全局字典中来推迟垃圾回收（通过新类型的 Python 析构函数将其删除）。所有这些都非常棘手，超出了这篇文章的范围；我只是提到它来说明这样一个事实，即通过构建上述技术，可以在 Julia 中实现相当复杂的跨语言调用行为。