Нет соответствующей функции для вызова «как» в Rcpp - PullRequest
3 голосов
/ 19 июня 2019

Я пытаюсь написать обертку R для подпрограмм FINUFFT для вычисления БПФ для ряда с неравной выборкой.У меня практически нет опыта работы с C / C ++, поэтому я работаю на примере, который сравнивает традиционное преобразование Фурье с NUFFT.Пример кода следующий.

// this is all you must include for the finufft lib...
#include "finufft.h"
#include <complex>

// also needed for this example...
#include <stdio.h>
#include <stdlib.h>
using namespace std;

int main(int argc, char* argv[])
/* Simple example of calling the FINUFFT library from C++, using plain
   arrays of C++ complex numbers, with a math test. Barnett 3/10/17
   Double-precision version (see example1d1f for single-precision)

   Compile with:
   g++ -fopenmp example1d1.cpp -I ../src ../lib-static/libfinufft.a -o example1d1  -lfftw3 -lfftw3_omp -lm
   or if you have built a single-core version:
   g++ example1d1.cpp -I ../src ../lib-static/libfinufft.a -o example1d1 -lfftw3 -lm

   Usage: ./example1d1
*/
{
  int M = 1e6;            // number of nonuniform points
  int N = 1e6;            // number of modes
  double acc = 1e-9;      // desired accuracy
  nufft_opts opts; finufft_default_opts(&opts);
  complex<double> I = complex<double>(0.0,1.0);  // the imaginary unit

  // generate some random nonuniform points (x) and complex strengths (c):
  double *x = (double *)malloc(sizeof(double)*M);
  complex<double>* c = (complex<double>*)malloc(sizeof(complex<double>)*M);
  for (int j=0; j<M; ++j) {
    x[j] = M_PI*(2*((double)rand()/RAND_MAX)-1);  // uniform random in [-pi,pi)
    c[j] = 2*((double)rand()/RAND_MAX)-1 + I*(2*((double)rand()/RAND_MAX)-1);
  }
  // allocate output array for the Fourier modes:
  complex<double>* F = (complex<double>*)malloc(sizeof(complex<double>)*N);

  // call the NUFFT (with iflag=+1): note N and M are typecast to BIGINT
  int ier = finufft1d1(M,x,c,+1,acc,N,F,opts);

  int n = 142519;   // check the answer just for this mode...
  complex<double> Ftest = complex<double>(0,0);
  for (int j=0; j<M; ++j)
    Ftest += c[j] * exp(I*(double)n*x[j]);
  int nout = n+N/2;        // index in output array for freq mode n
  double Fmax = 0.0;       // compute inf norm of F
  for (int m=0; m<N; ++m) {
    double aF = abs(F[m]);
    if (aF>Fmax) Fmax=aF;
  }
  double err = abs(F[nout] - Ftest)/Fmax;
  printf("1D type-1 NUFFT done. ier=%d, err in F[%d] rel to max(F) is %.3g\n",ier,n,err);

  free(x); free(c); free(F);
  return ier;
}

Многое из этого мне не нужно, например, генерация серии тестов и сравнение с традиционным БПФ.Кроме того, я хочу вернуть значения преобразования, а не просто код ошибки, указывающий на успех.Ниже приведен мой код.

#include "finufft.h"
#include <complex>
#include <Rcpp.h>
#include <stdlib.h>
using namespace Rcpp;
using namespace std;

// [[Rcpp::export]]

ComplexVector finufft(int M, NumericVector x, ComplexVector c, int N) {

  // From example code for finufft, sets precision and default options
  double acc = 1e-9;
  nufft_opts opts; finufft_default_opts(&opts);
  // allocate output array for the finufft routine:
  complex<double>* F = (complex<double>*)malloc(sizeof(complex<double>*)*N);

  // Change vector inputs from R types to C++ types
  double* xd = as< double* >(x);
  complex<double>* cd = as< complex<double>* >(c);      

  // call the NUFFT (with iflag=-1): note N and M are typecast to BIGINT
  int ier = finufft1d1(M,xd,cd,-1,acc,N,F,opts);
  ComplexVector Fd = as<ComplexVector>(*F);

  return Fd;
}

Когда я пытаюсь найти его в Rstudio, я получаю сообщение об ошибке "Нет подходящей функции для вызова" as (std :: complex * &) '",указывая на линию, декларирующую Fd в конце.Я считаю, что ошибка указывает на то, что либо функция «as» не определена (я знаю, что это ложь), либо аргумент «as» не является правильным типом.Примеры здесь включают пример, использующий «as» для преобразования в NumericVector, поэтому, если нет каких-либо сложностей со сложными значениями, я не понимаю, почему это должно быть проблемой здесь.

Я знаюесть потенциальные проблемы с использованием двух пространств имен, но я не думаю, что это проблема здесь.Мое предположение состоит в том, что есть проблема с тем, как я пытаюсь использовать указатели, но у меня нет опыта, чтобы идентифицировать его, и я не могу найти аналогичные примеры в Интернете, чтобы направлять меня.

1 Ответ

5 голосов
/ 19 июня 2019

Rcpp::as<T> преобразуется из типа данных R (SEXP) в тип данных C ++, например, Rcpp::ComplexVector.Это не соответствует вашей ситуации, когда вы пытаетесь преобразовать массив в стиле C в C ++.К счастью, Rcpp::Vector, который является основой для Rcpp::ComplexVector, имеет конструктор для этой задачи: Vector (InputIterator first, InputIterator last).Для другого направления (переход от C ++ к массиву в стиле C) вы можете использовать vector.begin() или &vector[0].

Однако для преобразования между Rcomplex* и std::complex<double>* требуется reinterpret_cast.Это не должно вызывать проблем, поскольку Rcomplex (он же complex double в C) и std::complex<doulbe> совместимы с .

Минимальный пример:

#include <Rcpp.h>
#include <complex>
using namespace Rcpp;

// [[Rcpp::export]]
ComplexVector foo(ComplexVector v) {
    std::complex<double>* F = reinterpret_cast<std::complex<double>*>(v.begin());
    int N = v.length();
    // do something with F
    ComplexVector Fd(reinterpret_cast<Rcomplex*>(F), 
                     reinterpret_cast<Rcomplex*>(F + N));
    return Fd;
}

/*** R
set.seed(42)
foo(runif(4)*(1+1i))
*/

Результат:

> Rcpp::sourceCpp('56675308/code.cpp')

> set.seed(42)

> foo(runif(4)*(1+1i))
[1] 0.9148060+0.9148060i 0.9370754+0.9370754i 0.2861395+0.2861395i 0.8304476+0.8304476i

Кстати, вы можете переместить эти reinterpret_cast с глаз долой, используя std::vector<std::complex<double>> в качестве аргумента и возвращаемого типа для вашей функции.Rcpp сделает все остальное за вас.Это также помогает избавиться от голых malloc:

#include <Rcpp.h>
// dummy function with reduced signature
int finufft1d1(int M, double *xd, std::complex<double> *cd, int N, std::complex<double> *Fd) {
    return 0;
}

// [[Rcpp::export]]
std::vector<std::complex<double>> finufft(int M, 
                                          std::vector<double> x, 
                                          std::vector<std::complex<double>> c, 
                                          int N) {

    // allocate output array for the finufft routine:
    std::vector<std::complex<double>> F(N);
    // Change vector inputs from R types to C++ types
    double* xd = x.data();
    std::complex<double>* cd = c.data();      
    std::complex<double>* Fd = F.data();      

    int ier = finufft1d1(M, xd, cd, N, Fd);

    return F;
}
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...