CxxWrap accumulator

accumulator.cc

#include "jlcxx/jlcxx.hpp"
#include "xtensor-julia/jltensor.hpp"     // Import the jltensor container definition
#include "xtensor-julia/jlarray.hpp"     // Import the jltensor container definition
#include "xtensor/xmath.hpp"              // xtensor import for the C++ universal functions
#include <inttypes.h>
#include <stdexcept>

#include "accumulator.hpp"
#include <sstream>

#include <boost/accumulators/statistics/min.hpp>
#include <boost/accumulators/statistics/max.hpp>

using boost::accumulators::min;
using boost::accumulators::max;

template<typename accum_t>
inline typename accum_t::sample_type do_mean (accum_t const& accum) {
  return mean (accum);
}

template<typename accum_t>
inline typename accum_t::sample_type do_count (accum_t const& accum) {
  return count (accum);
}

template<typename accum_t>
inline typename decomplexify<typename accum_t::sample_type>::type do_var (accum_t const& accum) {
  return abs_sqr (accum) - norm (mean (accum));
}

template<typename accum_t>
inline typename decomplexify<typename accum_t::sample_type>::type do_rms (accum_t const& accum) {
  return sqrt (do_var (accum));
}

template<typename accum_t>
std::string do_2nd_str (accum_t const& a) {
  std::ostringstream os;
  os << "mean: " << do_mean (a) << " rms: " << do_rms (a);
  return os.str();
}

template<typename accum_t>
std::string do_1st_str (accum_t const& a) {
  std::ostringstream os;
  os << do_mean (a);
  return os.str();
}

template<typename accum_t>
inline typename accum_t::sample_type do_min (accum_t const& accum) {
  return min (accum);
}

template<typename accum_t>
inline typename accum_t::sample_type do_max (accum_t const& accum) {
  return max (accum);
}

template<typename accum_t>
std::string do_extrema_str (accum_t const& a) {
  std::ostringstream os;
  os << "min: " << do_min (a) << " max: " << do_max (a);
  return os.str();
}

template<typename acc_t>
acc_t call (acc_t& a, typename acc_t::sample_type s) {
  a(s);
  return a;
}

template<typename el_t>
void expose_2nd (const char* name, jlcxx::Module& mod) {
  typedef accumulator_set<el_t,stats<tag::mean,tag::abs_sqr> > accum_t;

  mod.add_type<accum_t>(name)
    //.def (py::init<>()) only def constr
    // .method ("call", [](accum_t& a, xt::jlarray<el_t> const& s) {
    //     //   for (auto iter = s.storage_begin(); iter != s.storage_end(); ++iter) a(*iter);
    //     for (auto e : s) a(e);
    //     return a;
    //   })
    //.method ("call", [](accum_t& a, el_t s) { a(s); return a; })
    //.method ("call", &call<accum_t>)
    // .def ("__iadd__", [](accum_t& a, xt::pyarray<el_t> const& s) {
    //     // for (auto iter = s.storage_begin(); iter != s.storage_end(); ++iter) a(*iter);
    //     for (auto e : s) a(e);
    //     return a;
    //   }, py::is_operator())
    // .def ("__iadd__", [](accum_t& a, el_t s) { a(s); return a; }, py::is_operator())
    ;
  mod.method ("call", [](accum_t & a, el_t s) { a(s); return a;});
  mod.method ("mean", [](accum_t const& a) { return mean(a); });
  mod.method ("call", [](accum_t& a, xt::jlarray<el_t> const& s) {
      for (auto e : s) a(e);
      return a;
    });
  mod.method ("var", &do_var<accum_t>);
  // .def_property_readonly ("count", &do_count<accum_t>)
    // .def_property_readonly ("var", &do_var<accum_t>)
    // .def_property_readonly ("rms", &do_rms<accum_t>)
    // .def ("Mean", &do_mean<accum_t>)
    // .def ("Var", &do_var<accum_t>)
    // .def ("RMS", &do_rms<accum_t>)
    // .def ("__repr__", &do_2nd_str<accum_t>)


}

// template<typename el_t>
// static py::object expose_1st (const char* name, py::module m) {
//   typedef accumulator_set<el_t,stats<tag::mean> > accum_t;

//   py::class_<accum_t> pyclass (m, name);

//   pyclass
//     .def (py::init<>())
//     .def ("__call__", [](accum_t& a, xt::pyarray<el_t> const& s) {
//         for (auto e : s) a(e);
//         return a;
//       }, py::is_operator())
//     .def ("__call__", [](accum_t& a, el_t s) { a(s); return a; }, py::is_operator())
//     .def ("__iadd__", [](accum_t& a, xt::pyarray<el_t> const& s) {
//         for (auto e : s) a(e);
//         return a;
//       }, py::is_operator())
//     .def ("__iadd__", [](accum_t& a, el_t s) { a(s); return a; }, py::is_operator())
//     .def_property_readonly ("mean", &do_mean<accum_t>)
//     .def_property_readonly ("count", &do_count<accum_t>)
//     .def ("__repr__", &do_1st_str<accum_t>)
//     ;

//   return pyclass;
// }

// template<typename el_t>
// static void expose_extrema (const char* name, py::module m) {
//   typedef accumulator_set<el_t,stats<tag::min,tag::max> > accum_t;
//   py::class_<accum_t> (m, name)
//     .def (py::init<>())
//     .def ("__call__", [](accum_t& a, xt::pyarray<el_t> const& s) {
//         for (auto e : s) a(e);
//         return a;
//       }, py::is_operator())
//     .def ("__call__", [](accum_t& a, el_t s) { a(s); return a; }, py::is_operator())
//     .def ("__iadd__", [](accum_t& a, xt::pyarray<el_t> const& s) {
//         for (auto e : s) a(e);
//         return a;
//       }, py::is_operator())
//     .def ("__iadd__", [](accum_t& a, el_t s) { a(s); return a; }, py::is_operator())
//     .def_property_readonly ("min", &do_min<accum_t>)
//     .def_property_readonly ("max", &do_max<accum_t>)
//     .def ("__repr__", &do_extrema_str<accum_t>)
//     ;
// }

JULIA_CPP_MODULE_BEGIN(registry)
  jlcxx::Module& acc = registry.create_module("accumulator");

  expose_2nd<double> ("stat2nd_double", acc);
  expose_2nd<complex_t> ("stat2nd_complex", acc);

  // py::object stat1st_double_obj = expose_1st<double> ("stat1st_double", m);
  // py::object stat1st_complex_obj = expose_1st<complex_t> ("stat1st_complex", m);

  // py::object stat2nd_float_obj = expose_2nd<float> ("stat2nd_float", m);
  // py::object stat2nd_complex64_obj = expose_2nd<complex64_t> ("stat2nd_complex64", m);

  // py::object stat1st_float_obj = expose_1st<float> ("stat1st_float", m);
  // py::object stat1st_complex64_obj = expose_1st<complex64_t> ("stat1st_complex64", m);

  // // aliases
  // m.attr ("accumulator_double") = stat2nd_double_obj;
  // m.attr ("accumulator_complex") = stat2nd_complex_obj;

  // expose_extrema<float> ("extrema_float", m);
  // expose_extrema<double> ("extrema_double", m);
  // expose_extrema<int32_t> ("extrema_int32", m);
  // expose_extrema<int64_t> ("extrema_int64", m);

  // return m.ptr();
JULIA_CPP_MODULE_END

accumulator.hpp

#ifndef accumulator_hpp
#define accumulator_hpp

#define BOOST_NUMERIC_FUNCTIONAL_STD_COMPLEX_SUPPORT
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/moment.hpp>
#include <boost/range.hpp>
#include <complex>
#include "Complex.H"		// decomplexify

using namespace boost::accumulators;

typedef std::complex<double> complex_t;

inline double real (double x) { return x; }
inline double imag (double x) { return 0; }
inline double conj (double x) { return x; }
inline double norm (double x) { return x*x; }

template<typename T>
struct value_traits {
  typedef T value_type;
};

template<typename T>
struct value_traits<std::complex<T> > {
  typedef T value_type;
};

namespace boost { namespace accumulators { namespace impl {
      template<typename Sample>
      struct abs_sqr_accumulator
        : accumulator_base {
        //typedef typename Sample::value_type result_type;
	typedef typename value_traits<Sample>::value_type result_type;

        template<typename Args>
        abs_sqr_accumulator (Args const& args) {}

        template<typename Args>
        void operator () (Args const& args) {
          this->sum += real (args[sample]) * real (args[sample]) + imag (args[sample]) * imag (args[sample]);
        }

        template<typename Args>
        result_type result (Args const& args) const {
          return boost::numeric::average (this->sum, count (args));
        }

        result_type sum;
      };
    }}}
namespace boost { namespace accumulators { namespace tag {
      struct abs_sqr
        : depends_on< count >        // depends_on<> to specify dependencies
      {
        // Define a nested typedef called 'impl' that specifies which
        // accumulator implements this feature.
        typedef accumulators::impl::abs_sqr_accumulator< mpl::_1 > impl;
      };
    }}}

namespace boost {
  namespace accumulators {
    namespace extract {
      extractor< tag::abs_sqr> const abs_sqr = {};
    }
    using extract::abs_sqr;
  }
}

#endif

Complex.H

#ifndef Complex_H
#define Complex_H

#include <complex>
typedef std::complex<double> Complex;
typedef std::complex<double> complex_t;
typedef std::complex<float> complex64_t;
typedef std::complex<long double> complex256_t;

// decomplexify ---------------------------------------------------------------
template <typename T>
struct decomplexify_impl
{
  typedef T type;
};

template <typename ELT>
struct decomplexify_impl<std::complex<ELT> >
{
  typedef ELT type;
};

#include <type_traits>

template<typename T>
struct decomplexify
{
  typedef typename std::remove_reference<T>::type TT;
  typedef typename decomplexify_impl<TT>::type type;
};


// complexify -----------------------------------------------------------------
template <typename T>
struct complexify
{
  typedef std::complex<T> type;
};

template <typename ELT>
struct complexify<std::complex<ELT> >
{
  typedef std::complex<ELT> type;
};

#endif

test_acc.jl

using CxxWrap
using Xtensor
wrap_modules(Xtensor._l_tensors)
wrap_modules("./libaccumulator")

using accumulator: stat2nd_double, call, mean, stat2nd_complex, var
s = stat2nd_double()
call(s, 2.0)

#s2 = stat2nd_complex()
#call(s2, complex(2.0))

call(s, ones(10))