#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/fisher_f.hpp>

#include <boost/math/special_functions/gamma.hpp>

#include <boost/timer.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>

#include <string>
#include <vector>
#include <algorithm>
#include <iterator>
#include <functional>
#include <iostream>
#include <cmath>

extern "C" {
#include "Rmath.h"
}



template <class Dist>
struct callpdf : public std::unary_function<double, double> {
  callpdf(Dist dist = Dist()) : _dist(dist) {}
  double operator()(double x) const {
    return boost::math::pdf(_dist, x);
  }

  Dist _dist;
};
template <class Dist>
struct callcdf : public std::unary_function<double, double> {
  callcdf(Dist dist = Dist()) : _dist(dist) {}
  double operator()(double x) const {
    return boost::math::cdf(_dist, x);
  }

  Dist _dist;
};
template <class Dist>
struct callquantile : public std::unary_function<double, double> {
  callquantile(Dist dist = Dist()) : _dist(dist) {}
  double operator()(double x) const {
    return boost::math::quantile(_dist, x);
  }

  Dist _dist;
};

template <typename TF>
double time(TF f, double x, size_t N = 1000000) {
  boost::timer clock;
  clock.restart();
  for(size_t i = 0; i < N; ++i)
    double y = f(x);
  return clock.elapsed();
}

int main()
{
  double x = .0;
  size_t N = 5000000;
  double elapsed = time(boost::lambda::bind(dnorm, boost::lambda::_1, 0.0, 1.0, 0), x, N);
  std::cout << "R norm pdf: " << elapsed << std::endl;
  elapsed = time(callpdf<boost::math::normal>(), x, N);
  std::cout << "boost::math norm pdf: " << elapsed << std::endl;

  elapsed = time(boost::lambda::bind(pnorm, boost::lambda::_1, 0.0, 1.0, 1, 0), x, N);
  std::cout << "R norm cdf: " << elapsed << std::endl;
  elapsed = time(callcdf<boost::math::normal>(), x, N);
  std::cout << "boost::math norm cdf: " << elapsed << std::endl;

  elapsed = time(boost::lambda::bind(qnorm, boost::lambda::_1, 0.0, 1.0, 1, 0), 0.1, N);
  std::cout << "R norm quantile: " << elapsed << std::endl;
  elapsed = time(callquantile<boost::math::normal>(), 0.1, N);
  std::cout << "boost::math norm quantile: " << elapsed << std::endl;

  elapsed = time(boost::lambda::bind(lgammafn, boost::lambda::_1), 1.0, N);
  std::cout << "R lgamma: " << elapsed << std::endl;
  double (*pf_d_d)(double) = &boost::math::lgamma;
  elapsed = time(boost::lambda::bind(pf_d_d, boost::lambda::_1), 1.0, N);
  std::cout << "boost::math lgamma: " << elapsed << std::endl;

  N = 1000000;
  boost::math::students_t t(10);
  elapsed = time(boost::lambda::bind(dt, boost::lambda::_1, 10, 0), x, N);
  std::cout << "R t pdf: " << elapsed << std::endl;
  elapsed = time(callpdf<boost::math::students_t>(t), x, N);
  std::cout << "boost::math t pdf: " << elapsed << std::endl;

  elapsed = time(boost::lambda::bind(pt, boost::lambda::_1, 10, 1, 0), x, N);
  std::cout << "R t cdf: " << elapsed << std::endl;
  elapsed = time(callcdf<boost::math::students_t>(t), x, N);
  std::cout << "boost::math t cdf: " << elapsed << std::endl;

  elapsed = time(boost::lambda::bind(qt, boost::lambda::_1, 10, 1, 0), 0.1, N);
  std::cout << "R t quantile: " << elapsed << std::endl;
  elapsed = time(callquantile<boost::math::students_t>(t), 0.1, N);
  std::cout << "boost::math t quantile: " << elapsed << std::endl;

  boost::math::fisher_f f(100,200);
  elapsed = time(boost::lambda::bind(df, boost::lambda::_1, 100., 200., 0), x, N);
  std::cout << "R F pdf: " << elapsed << std::endl;
  elapsed = time(callpdf<boost::math::fisher_f>(f), x, N);
  std::cout << "boost::math F pdf: " << elapsed << std::endl;
  elapsed = time(boost::lambda::bind(pf, boost::lambda::_1, 100., 200., 1, 0), x, N);
  std::cout << "R F cdf: " << elapsed << std::endl;
  elapsed = time(callcdf<boost::math::fisher_f>(f), x, N);
  std::cout << "boost::math F cdf: " << elapsed << std::endl;
  elapsed = time(boost::lambda::bind(qf, boost::lambda::_1, 100., 200., 1, 0), 0.1, N);
  std::cout << "R F quantile: " << elapsed << std::endl;
  elapsed = time(callquantile<boost::math::fisher_f>(f), 0.1, N);
  std::cout << "boost::math F quantile: " << elapsed << std::endl;

  return 0;
}