Basic support vectorised inputs

NEWS.md

@@ -1,3 +1,6 @@
+# boostmath 1.5.0
+  * Functions (excluding distributions and polynomials) now support vectorised inputs
+
 # boostmath 1.4.0
   * Added `logit` and `logistic_sigmoid` functions
   * Added `log_pi` constant

src/boostmath/as_cpp.hpp

@@ -0,0 +1,125 @@
+#ifndef BOOSTMATH_AS_CPP_HPP
+#define BOOSTMATH_AS_CPP_HPP
+
+#include <cpp11.hpp>
+#include <complex>
+#include "./type_traits.hpp"
+
+namespace boostmath {
+  namespace internal {
+    template <typename T, std::size_t... Is, std::enable_if_t<std::negation<is_integral_scalar_type<T>>::value>* = nullptr>
+    T as_cpp_arithmetic_container_impl(SEXP x, std::index_sequence<Is...>) {
+      return {REAL_ELT(x, Is)...};
+    }
+
+    template <typename T, std::size_t... Is, std::enable_if_t<is_integral_scalar_type<T>::value>* = nullptr>
+    T as_cpp_arithmetic_container_impl(SEXP x, std::index_sequence<Is...>) {
+      return {INTEGER_ELT(x, Is)...};
+    }
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::complex<double>>::value>* = nullptr>
+  inline std::complex<double> as_cpp(SEXP x) {
+    Rcomplex r_complex = COMPLEX_ELT(x, 0);
+    return {r_complex.r, r_complex.i};
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::complex<double>>>::value>* = nullptr>
+  inline std::vector<std::complex<double>> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<std::complex<double>> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      Rcomplex r_complex = COMPLEX_ELT(x, i);
+      result[i] = {r_complex.r, r_complex.i};
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<uint64_t>>::value>* = nullptr>
+  inline std::vector<uint64_t> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<uint64_t> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<uint64_t>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int64_t>>::value>* = nullptr>
+  inline std::vector<int64_t> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<int64_t> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<int64_t>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int>>::value>* = nullptr>
+  inline std::vector<int> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<int> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<int>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<double>>::value>* = nullptr>
+  inline std::vector<double> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::vector<double> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result[i] = static_cast<double>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::list<double>>::value>* = nullptr>
+  inline std::list<double> as_cpp(SEXP x) {
+    R_xlen_t n = Rf_xlength(x);
+    std::list<double> result(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      result.push_back(static_cast<double>(Rf_isInteger(x) ? INTEGER_ELT(x, i) : REAL_ELT(x, i)));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::vector<double>>>::value>* = nullptr>
+  inline std::vector<std::vector<double>> as_cpp(SEXP x) {
+    if (Rf_isNewList(x)) {
+      const R_xlen_t n = Rf_xlength(x);
+      std::vector<std::vector<double>> result(n);
+      for (R_xlen_t i = 0; i < n; ++i) {
+        SEXP xi = VECTOR_ELT(x, i);
+        result[i] = as_cpp<std::vector<double>>(xi);
+      }
+      return result;
+    }
+    return {as_cpp<std::vector<double>>(x)};
+  }
+
+  template <typename T, std::enable_if_t<is_arithmetic_container<T>::value>* = nullptr>
+  inline T as_cpp(SEXP x) {
+    return internal::as_cpp_arithmetic_container_impl<T>(x, std::make_index_sequence<std::tuple_size<T>::value>{});
+  }
+
+  template <typename T, std::enable_if_t<is_vector_of_arithmetic_containers<T>::value>* = nullptr>
+  inline T as_cpp(SEXP x) {
+    using array_type = typename T::value_type;
+    R_xlen_t n = Rf_xlength(x);
+    T result;
+    result.reserve(n);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SEXP xi = VECTOR_ELT(x, i);
+      result.push_back(as_cpp<array_type>(xi));
+    }
+    return result;
+  }
+
+  template <typename T, std::enable_if_t<is_cpp11<T>::value>* = nullptr>
+  inline T as_cpp(SEXP x) {
+    return cpp11::as_cpp<T>(x);
+  }
+}
+#endif

src/boostmath/as_sexp.hpp

@@ -0,0 +1,152 @@
+#ifndef BOOSTMATH_AS_SEXP_HPP
+#define BOOSTMATH_AS_SEXP_HPP
+
+#include <cpp11.hpp>
+#include <complex>
+#include "./type_traits.hpp"
+#include "Rinternals.h"
+
+namespace boostmath {
+  namespace internal {
+    template <typename T, std::size_t... Is, std::enable_if_t<std::negation<is_integral_scalar_type<T>>::value>* = nullptr>
+    auto as_sexp_arithmetic_container_impl(const T& x, std::index_sequence<Is...>) {
+      SEXP data = cpp11::safe[Rf_allocVector](REALSXP, sizeof...(Is));
+      double* ptr = REAL(data);
+      ((ptr[Is] = std::get<Is>(x)), ...);
+      return data;
+    }
+
+    template <typename T, std::size_t... Is,  std::enable_if_t<is_integral_scalar_type<T>::value>* = nullptr>
+    auto as_sexp_arithmetic_container_impl(const T& x, std::index_sequence<Is...>) {
+      SEXP data = cpp11::safe[Rf_allocVector](INTSXP, sizeof...(Is));
+      int* ptr = INTEGER(data);
+      ((ptr[Is] = std::get<Is>(x)), ...);
+      return data;
+    }
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::complex<double>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    SEXP data = cpp11::safe[Rf_allocVector](CPLXSXP, 1);
+    Rcomplex* r_complex = COMPLEX(data);
+    r_complex->r = x.real();
+    r_complex->i = x.imag();
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::complex<double>>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](CPLXSXP, n);
+    Rcomplex* r_complex = COMPLEX(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      r_complex[i].r = x[i].real();
+      r_complex[i].i = x[i].imag();
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::list<double>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](REALSXP, n);
+    double* ptr = REAL(data);
+    size_t i = 0;
+    for (const auto& value : x) {
+      ptr[i++] = value;
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::list<double>>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
+    PROTECT(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
+    }
+    UNPROTECT(1);
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<uint64_t>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](INTSXP, n);
+    int* ptr = INTEGER(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<int>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int64_t>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](INTSXP, n);
+    int* ptr = INTEGER(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<int>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<int>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](INTSXP, n);
+    int* ptr = INTEGER(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<int>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<double>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](REALSXP, n);
+    double* ptr = REAL(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      ptr[i] = static_cast<double>(x[i]);
+    }
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<std::is_same<T, std::vector<std::vector<double>>>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
+    PROTECT(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
+    }
+    UNPROTECT(1);
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<is_arithmetic_container<T>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    return internal::as_sexp_arithmetic_container_impl<T>(x, std::make_index_sequence<std::tuple_size<T>::value>{});
+  }
+
+  template <typename T, std::enable_if_t<is_vector_of_arithmetic_containers<T>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    R_xlen_t n = x.size();
+    SEXP data = cpp11::safe[Rf_allocVector](VECSXP, n);
+    PROTECT(data);
+    for (R_xlen_t i = 0; i < n; ++i) {
+      SET_VECTOR_ELT(data, i, as_sexp(x[i]));
+    }
+    UNPROTECT(1);
+    return data;
+  }
+
+  template <typename T, std::enable_if_t<is_cpp11<T>::value>* = nullptr>
+  inline SEXP as_sexp(const T& x) {
+    return cpp11::as_sexp(x);
+  }
+}
+
+#endif // BOOSTMATH_SEXP_HPP

src/boostmath/boost_function.hpp

@@ -0,0 +1,58 @@
+#ifndef BOOSTMATH_UTILS_HPP
+#define BOOSTMATH_UTILS_HPP
+
+#include <tuple>
+#include "as_sexp.hpp"
+#include "scalar_container.hpp"
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace boostmath {
+  namespace internal {
+    template <typename TypesTuple, typename ArgsTupleT, std::size_t... Is>
+    auto make_containers_tuple_impl(ArgsTupleT&& args, std::index_sequence<Is...>) {
+      return std::make_tuple(boostmath::scalar_container<std::tuple_element_t<Is, TypesTuple>>(std::get<Is>(args))...);
+    }
+  }
+
+  template <typename TypesTupleT, typename ArgsTupleT>
+  auto make_containers_tuple(ArgsTupleT&& args) {
+    constexpr auto tuple_index_seq = std::make_index_sequence<std::tuple_size_v<std::decay_t<ArgsTupleT>>>{};
+    return internal::make_containers_tuple_impl<TypesTupleT>(std::forward<ArgsTupleT>(args), tuple_index_seq);
+  }
+
+  template <typename F, typename... ArgsT>
+  struct invoke_rtn {
+    using type = std::invoke_result_t<F, ArgsT...>;
+  };
+
+  template <typename F, typename... ArgsT>
+  struct invoke_rtn<F, std::tuple<ArgsT...>> {
+    using type = std::invoke_result_t<F, ArgsT...>;
+  };
+
+  template <typename F, typename... ArgsT>
+  using invoke_rtn_t = typename invoke_rtn<F, ArgsT...>::type;
+
+  template <typename TypesTuple, typename F, typename ArgsTupleT>
+  SEXP boost_function(F&& f, ArgsTupleT&& args) {
+    const auto containers_tuple = make_containers_tuple<TypesTuple>(std::forward<ArgsTupleT>(args));
+
+    using containers_tuple_t = decltype(containers_tuple);
+    size_t max_n = std::apply([](auto&&... conts) { return std::max({conts.size()...}); },
+                            std::forward<containers_tuple_t>(containers_tuple));
+
+    std::vector<invoke_rtn_t<F, TypesTuple>> results(max_n);
+    for (size_t i = 0; i < max_n; ++i) {
+      results[i] = std::apply([&](auto&&... conts){ return f(conts[i]...); },
+                              std::forward<containers_tuple_t>(containers_tuple));
+    }
+    if (max_n == 1) {
+      return boostmath::as_sexp(results[0]);
+    }
+    return boostmath::as_sexp(results);
+  }
+}
+
+#endif // BOOSTMATH_UTILS_HPP