在r：dyn.load中实现非线性数据拟合的GSL算法

Question

我正在尝试使用r使用rcpp在r中实现曲线拟合的GSL非线性最小二乘算法。 这个问题接近我在这里提出的先前问题：在非线性最小二乘GSL中拟合函数的固定参数

我尝试实现基于GSL的非线性最小二乘如果我的目标是估计用于拟合数据的所有参数，算法就成功了。当我尝试在在非线性最小二乘GSL中的拟合函数的固定参数用于固定函数的某些参数。

当我sourcecpp我的代码（按照我的上一个问题）进行了调整时，我确实会收到以下错误：

Error in dyn.load("/private/var/folders/pq/hxwd9my563q_qpy4rbrlgkmw0000gn/T/RtmpRKdn9f/sourceCpp-x86_64-apple-darwin17.0-1.0.8.3/sourcecpp_a60fe49985e/sourceCpp_80.so") : 
  unable to load shared object '/private/var/folders/pq/hxwd9my563q_qpy4rbrlgkmw0000gn/T/RtmpRKdn9f/sourceCpp-x86_64-apple-darwin17.0-1.0.8.3/sourcecpp_a60fe49985e/sourceCpp_80.so':
  dlopen(/private/var/folders/pq/hxwd9my563q_qpy4rbrlgkmw0000gn/T/RtmpRKdn9f/sourceCpp-x86_64-apple-darwin17.0-1.0.8.3/sourcecpp_a60fe49985e/sourceCpp_80.so, 0x0006): symbol not found in flat namespace '__Z28internal_make_gsl_vector_ptrILm3EEP10gsl_vectorRKNSt3__15arrayIdXT_EEE' 

这是用于执行非线性最小二乘数据拟合的C ++包装器代码：

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
// [[Rcpp::depends(RcppGSL)]]
// [[Rcpp::depends(BH)]]
#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <RcppNumerical.h>
#include <RcppGSL.h>
#include <array>
#include <Rcpp.h>
#include <iostream>
#include <vector>
#include <cassert>
#include <functional>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlinear.h>

using namespace std;
using namespace Rcpp;
using namespace Numer;

template <class R, class... ARGS>
struct function_ripper {
    static constexpr size_t n_args = sizeof...(ARGS);
};

template <class R, class... ARGS>
auto constexpr n_params(R (ARGS...) )
{
    return function_ripper<R, ARGS...>();
}

template <typename F, size_t... Is>
auto gen_tuple_impl(F func, std::index_sequence<Is...> )
{
    return std::make_tuple(func(Is)...);
}

template <size_t N, typename F>
auto gen_tuple(F func)
{
    return gen_tuple_impl(func, std::make_index_sequence<N>{} );
}

auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
                           gsl_multifit_nlinear_parameters *params) -> std::vector<double>
{
    // This specifies a trust region method
    const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
    const size_t max_iter = 200;
    const double xtol = 1.0e-8;
    const double gtol = 1.0e-8;
    const double ftol = 1.0e-8;

    auto *work = gsl_multifit_nlinear_alloc(T, params, fdf->n, fdf->p);
    int info;

    // initialize solver
    gsl_multifit_nlinear_init(initial_params, fdf, work);
    //iterate until convergence
    gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, nullptr, nullptr, &info, work);

    // result will be stored here
    gsl_vector * y    = gsl_multifit_nlinear_position(work);
    auto result = std::vector<double>(initial_params->size);

    for(int i = 0; i < result.size(); i++)
    {
        result[i] = gsl_vector_get(y, i);
    }

    auto niter = gsl_multifit_nlinear_niter(work);
    auto nfev  = fdf->nevalf;
    auto njev  = fdf->nevaldf;
    auto naev  = fdf->nevalfvv;
    gsl_multifit_nlinear_free(work);
    gsl_vector_free(initial_params);
    return result;
}

auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*
{
    auto* result = gsl_vector_alloc(vec.size());
    int i = 0;
    for(const auto e: vec)
    {
        gsl_vector_set(result, i, e);
        i++;
    }
    return result;
}

template<typename C1>
struct fit_data
{
    const std::vector<double>& t;
    const std::vector<double>& y;
    // the actual function to be fitted
    C1 f;
};

template<typename FitData, int n_params>
int internal_f(const gsl_vector* x, void* params, gsl_vector *f)
{
    auto* d  = static_cast<FitData*>(params);
    // Convert the parameter values from gsl_vector (in x) into std::tuple
    auto init_args = [x](int index)
    {
        return gsl_vector_get(x, index);
    };
    auto parameters = gen_tuple<n_params>(init_args);

    // Calculate the error for each...
    for (size_t i = 0; i < d->t.size(); ++i)
    {
        double ti = d->t[i];
        double yi = d->y[i];
        auto func = [ti, &d](auto ...xs)
        {
            // call the actual function to be fitted
            return d->f(ti, xs...);
        };
        auto y = std::apply(func, parameters);
        gsl_vector_set(f, i, yi - y);
    }
    return GSL_SUCCESS;
}

using func_f_type   = int (*) (const gsl_vector*, void*, gsl_vector*);
using func_df_type  = int (*) (const gsl_vector*, void*, gsl_matrix*);
using func_fvv_type = int (*) (const gsl_vector*, const gsl_vector *, void *, gsl_vector *);

template<auto n>
auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*;


auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params) -> std::vector<double>;

template<typename C1>
auto curve_fit_impl(func_f_type f, func_df_type df, func_fvv_type fvv, gsl_vector* initial_params, fit_data<C1>& fd) -> std::vector<double>
{
    assert(fd.t.size() == fd.y.size());

    auto fdf = gsl_multifit_nlinear_fdf();
    auto fdf_params = gsl_multifit_nlinear_default_parameters();

    fdf.f   = f;
    fdf.df  = df;
    fdf.fvv = fvv;
    fdf.n   = fd.t.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    // "This selects the Levenberg-Marquardt algorithm with geodesic acceleration."
    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    return internal_solve_system(initial_params, &fdf, &fdf_params);
}

template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x, const std::vector<double>& y) -> std::vector<double>
{
    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x, y, f};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

这些是我用来适合的功能数据：

// [[Rcpp::export]]
double gaussian(double x, double a, double b, double c){
    const double z = (x - b) / c;
    return a * std::exp(-0.5 * z * z);
}

struct gaussian_fixed_a{
    double a;
    gaussian_fixed_a(double a) : a{a} {}
    double operator()(double x, double b, double c) const { return gaussian(x, a, b, c); }
};

// [[Rcpp::export]]
Rcpp::List fittingTest(const std::vector<double> xs,const std::vector<double> ys, const double a){

      gaussian_fixed_a g(a);
      auto r = curve_fit(g, std::array{0.444, 0.11}, xs, ys);
      return Rcpp::List::create(Rcpp::Named("b") = r[0],
                                  Rcpp::Named("c") = r[1]);
        }

我的代码在哪里创建链接问题？

Answer 1

通过关注@zkoza的答案来解决该错误，该答案指定了一个编译时数组，其中参数数是从数组长度自动推导的。在第81-92行中：

template<auto n>
auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*
{
    auto* result = gsl_vector_alloc(vec.size());
    int i = 0;
    for(const auto e: vec)
    {
        gsl_vector_set(result, i, e);
        i++;
    }
    return result;
}

现在的完整代码看起来像：

#include <RcppGSL.h>
#include <array>
#include <Rcpp.h>
#include <cmath>
#include <iostream>
#include <cstdlib>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <vector>
#include <cassert>
#include <functional>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlinear.h>

using namespace std;
using namespace Rcpp;

// [[Rcpp::depends(RcppGSL)]]

template <typename F, size_t... Is>
auto gen_tuple_impl(F func, std::index_sequence<Is...> )
{
    return std::make_tuple(func(Is)...);
}

template <size_t N, typename F>
auto gen_tuple(F func)
{
    return gen_tuple_impl(func, std::make_index_sequence<N>{} );
}

template <class R, class... ARGS>
struct function_ripper {
    static constexpr size_t n_args = sizeof...(ARGS);
};

template <class R, class... ARGS>
auto constexpr n_params(R (ARGS...) )
{
    return function_ripper<R, ARGS...>();
}

auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params) -> std::vector<double>
{
  // This specifies a trust region method
  const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
  const size_t max_iter = 200;
  const double xtol = 1.0e-8;
  const double gtol = 1.0e-8;
  const double ftol = 1.0e-8;

  auto *work = gsl_multifit_nlinear_alloc(T, params, fdf->n, fdf->p);
  int info;

  // initialize solver
  gsl_multifit_nlinear_init(initial_params, fdf, work);
  //iterate until convergence
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, nullptr, nullptr, &info, work);

  // result will be stored here
  gsl_vector * y    = gsl_multifit_nlinear_position(work);
  auto result = std::vector<double>(initial_params->size);

  for(int i = 0; i < result.size(); i++)
  {
    result[i] = gsl_vector_get(y, i);
  }

  auto niter = gsl_multifit_nlinear_niter(work);
  auto nfev  = fdf->nevalf;
  auto njev  = fdf->nevaldf;
  auto naev  = fdf->nevalfvv;

  gsl_multifit_nlinear_free(work);
  gsl_vector_free(initial_params);
  return result;
}

template<auto n>
auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*
{
    auto* result = gsl_vector_alloc(vec.size());
    int i = 0;
    for(const auto e: vec)
    {
        gsl_vector_set(result, i, e);
        i++;
    }
    return result;
}

template<typename C1>
struct fit_data
{
    const std::vector<double>& t;
    const std::vector<double>& y;
    // the actual function to be fitted
    C1 f;
};


template<typename FitData, int n_params>
int internal_f(const gsl_vector* x, void* params, gsl_vector *f)
{
    auto* d  = static_cast<FitData*>(params);
    // Convert the parameter values from gsl_vector (in x) into std::tuple
    auto init_args = [x](int index)
    {
        return gsl_vector_get(x, index);
    };
    auto parameters = gen_tuple<n_params>(init_args);

    // Calculate the error for each...
    for (size_t i = 0; i < d->t.size(); ++i)
    {
        double ti = d->t[i];
        double yi = d->y[i];
        auto func = [ti, &d](auto ...xs)
        {
            // call the actual function to be fitted
            return d->f(ti, xs...);
        };
        auto y = std::apply(func, parameters);
        gsl_vector_set(f, i, yi - y);
    }
    return GSL_SUCCESS;
}

using func_f_type   = int (*) (const gsl_vector*, void*, gsl_vector*);
using func_df_type  = int (*) (const gsl_vector*, void*, gsl_matrix*);
using func_fvv_type = int (*) (const gsl_vector*, const gsl_vector *, void *, gsl_vector *);

template<auto n>
auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*;

auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params) -> std::vector<double>;

template<typename C1>
auto curve_fit_impl(func_f_type f, func_df_type df, func_fvv_type fvv, gsl_vector* initial_params, fit_data<C1>& fd) -> std::vector<double>
{
    assert(fd.t.size() == fd.y.size());

    auto fdf = gsl_multifit_nlinear_fdf();
    auto fdf_params = gsl_multifit_nlinear_default_parameters();

    fdf.f   = f;
    fdf.df  = df;
    fdf.fvv = fvv;
    fdf.n   = fd.t.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    // "This selects the Levenberg-Marquardt algorithm with geodesic acceleration."
    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    return internal_solve_system(initial_params, &fdf, &fdf_params);
}

template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x, const std::vector<double>& y) -> std::vector<double>
{
    // We can't pass lambdas without convert to std::function.
    //constexpr auto n = 3;//decltype(n_params(f))::n_args - 5;
    //constexpr auto n = 2;
    assert(initial_params.size() == n);

    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x, y, f};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

同一 *.cpp文件中还包括Gaussian函数，Gaussian functor gaussian_fixed_a和>拟合操作函数，如问题所述。

在R中，我将以以下方式测试fittingtest函数：

Sys.setenv("PKG_CXXFLAGS"="-std=c++17")

require(Rcpp)
require(RcppGSL)

sourceCpp('./testingCode.cpp')

x_list = seq(2,14,1)

###START: simulate Data under gaussian function###
res = data.frame(geo=x_list,val=NA)

for(i in 1:length(x_list)){
  res$val[i]= gaussian(x_list[i], a=8, b=0.4, c=5)
}
###END: simulate Data under gaussian function###

fittingTest(xs=res$geo,ys=res$val,a=8)

并且在修复a = 8的同时，我获得的结果是：

$b
[1] 0.4

$c
[1] 5

用于目视检查结果，您可以看到模拟和fitting数据完美重叠（请注意，在此示例中，我没有在模拟数据中添加任何“噪声”。