#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/fisher_f.hpp>
#include <boost/math/distributions/chi_squared.hpp>
#include <boost/math/distributions/exponential.hpp>
#include <boost/math/distributions/gamma.hpp>
#include <boost/math/distributions/lognormal.hpp>
#include <boost/math/distributions/cauchy.hpp>
#include <boost/math/distributions/weibull.hpp>
#include <boost/math/distributions/uniform.hpp>
#include <boost/math/distributions/beta.hpp>
#include <boost/math/distributions/binomial.hpp>
#include <boost/math/distributions/poisson.hpp>
#include <boost/math/distributions/negative_binomial.hpp>

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

#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/iterator/counting_iterator.hpp>

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

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

double _abs(double x) {
  return std::abs(x);
}

void compare(const std::vector<double> & x, const std::vector<double> & y, const std::string & info) {
  std::vector<double> diff;
  std::transform(x.begin(), x.end(), y.begin(), std::back_inserter(diff), boost::lambda::_1 - boost::lambda::_2);
  std::transform(diff.begin(), diff.end(), diff.begin(), _abs);
  std::cout << info << " max diff: " << *std::max_element(diff.begin(), diff.end()) << std::endl;
  std::transform(diff.begin(), diff.end(), x.begin(), diff.begin(), boost::lambda::_1/boost::lambda::_2);
  std::cout << info << " max rel diff: " << *std::max_element(diff.begin(), diff.end()) << std::endl;
}

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;
};


int main()
{
  std::vector<double> xs;
  for(double x = -4.0; x <= 4.0; x += 0.01)
    xs.push_back(x);
  std::vector<double> res1(xs);
  std::vector<double> res2(xs);
//   double (*pf1)(const boost::math::normal&, const double&) = &boost::math::pdf;
//   std::transform(xs.begin(), xs.end(), res1.begin(), boost::lambda::bind(pf1, boost::math::normal(), boost::lambda::_1));
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::normal>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dnorm, boost::lambda::_1, 0.0, 1.0, 0));
  compare(res1, res2, "norm pdf");
  boost::math::students_t t(10);
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::students_t>(t));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dt, boost::lambda::_1, 10, 0));
  compare(res1, res2, "t pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::normal>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pnorm, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "norm cdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::students_t>(t));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pt, boost::lambda::_1, 10, 1, 0));
  compare(res1, res2, "t cdf");

  xs.clear();
  for(double x = 0.01; x <= 100.0; x += 0.01)
    xs.push_back(x);
  res2 = res1 = xs;

  boost::math::fisher_f f(100,200);
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::fisher_f>(f));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(df, boost::lambda::_1, 100, 200, 0));
  compare(res1, res2, "F pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::fisher_f>(f));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pf, boost::lambda::_1, 100., 200., 1, 0));
  compare(res1, res2, "F cdf");

  boost::math::chi_squared chi2(100);
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::chi_squared>(chi2));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dchisq, boost::lambda::_1, 100, 0));
  compare(res1, res2, "chi2 pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::chi_squared>(chi2));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pchisq, boost::lambda::_1, 100, 1, 0));
  compare(res1, res2, "chi2 cdf");

  boost::math::exponential expd(2);
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::exponential>(expd));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dexp, boost::lambda::_1, .5, 0));
  compare(res1, res2, "exp pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::exponential>(expd));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pexp, boost::lambda::_1, .5, 1, 0));
  compare(res1, res2, "exp cdf");

  boost::math::gamma_distribution<> gammad(3,2);
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::gamma_distribution<> >(gammad));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dgamma, boost::lambda::_1, 3, 2, 0));
  compare(res1, res2, "gamma pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::gamma_distribution<> >(gammad));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pgamma, boost::lambda::_1, 3, 2, 1, 0));
  compare(res1, res2, "gamma cdf");

  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::lognormal>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dlnorm, boost::lambda::_1, 0.0, 1.0, 0));
  compare(res1, res2, "lognorm pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::lognormal>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(plnorm, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "lognorm cdf");

  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::cauchy>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dcauchy, boost::lambda::_1, .0, 1., 0));
  compare(res1, res2, "cauchy pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::cauchy>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pcauchy, boost::lambda::_1, .0, 1., 1, 0));
  compare(res1, res2, "cauchy cdf");

  boost::math::weibull weibulld(2);
  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::weibull>(weibulld));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dweibull, boost::lambda::_1, 2., 1., 0));
  compare(res1, res2, "weibull pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::weibull>(weibulld));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(pweibull, boost::lambda::_1, 2., 1., 1, 0));
  compare(res1, res2, "weibull cdf");

  std::vector<double> qs;
  for(double q = 0.001; q <= 0.999; q += 0.001)
    qs.push_back(q);
  res2 = res1 = qs;
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::normal>());
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qnorm, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "norm quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::students_t>(t));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qt, boost::lambda::_1, 10, 1, 0));
  compare(res1, res2, "t quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::fisher_f>(f));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qf, boost::lambda::_1, 100, 200, 1, 0));
  compare(res1, res2, "F quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::fisher_f>(f));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qf, boost::lambda::_1, 100, 200, 1, 0));
  compare(res1, res2, "F quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::chi_squared>(chi2));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qchisq, boost::lambda::_1, 100, 1, 0));
  compare(res1, res2, "chi2 quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::exponential>(expd));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qexp, boost::lambda::_1, .5, 1, 0));
  compare(res1, res2, "exp quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::gamma_distribution<> >(gammad));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qgamma, boost::lambda::_1, 3, 2, 1, 0));
  compare(res1, res2, "gamma quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::lognormal>());
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qlnorm, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "lognorm quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::cauchy>());
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qcauchy, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "cauchy quantile");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::weibull>(weibulld));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qweibull, boost::lambda::_1, 2.0, 1.0, 1, 0));
  compare(res1, res2, "weibull quantile");

  std::transform(xs.begin(), xs.end(), res1.begin(), callpdf<boost::math::uniform>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(dunif, boost::lambda::_1, 0.0, 1.0, 0));
  compare(res1, res2, "uniform pdf");
  std::transform(xs.begin(), xs.end(), res1.begin(), callcdf<boost::math::uniform>());
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(punif, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "uniform cdf");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::uniform>());
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qunif, boost::lambda::_1, 0.0, 1.0, 1, 0));
  compare(res1, res2, "uniform quantile");

  boost::math::beta_distribution<> betad(2,5);
  std::transform(qs.begin(), qs.end(), res1.begin(), callpdf<boost::math::beta_distribution<> >(betad));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(dbeta, boost::lambda::_1, 2, 5, 0));
  compare(res1, res2, "beta pdf");
  std::transform(qs.begin(), qs.end(), res1.begin(), callcdf<boost::math::beta_distribution<> >(betad));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(pbeta, boost::lambda::_1, 2, 5, 1, 0));
  compare(res1, res2, "beta cdf");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::beta_distribution<> >(betad));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qbeta, boost::lambda::_1, 2, 5, 1, 0));
  compare(res1, res2, "beta quantile");

  size_t N = qs.size()-1;
  boost::math::binomial_distribution<> binomd(N, 0.2);
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res1.begin(),
                 callpdf<boost::math::binomial_distribution<> >(binomd));
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res2.begin(),
                 boost::lambda::bind(dbinom, boost::lambda::_1, N, .2, 0));
  compare(res1, res2, "binom pdf");
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res1.begin(),
                 callcdf<boost::math::binomial_distribution<> >(binomd));
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res2.begin(),
                 boost::lambda::bind(pbinom, boost::lambda::_1, N, .2, 1, 0));
  compare(res1, res2, "binom cdf");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::binomial_distribution<> >(binomd));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qbinom, boost::lambda::_1, N, .2, 1, 0));
  compare(res1, res2, "binom quantile");

  boost::math::poisson poissond(0.2);
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res1.begin(),
                 callpdf<boost::math::poisson >(poissond));
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res2.begin(),
                 boost::lambda::bind(dpois, boost::lambda::_1, .2, 0));
  compare(res1, res2, "poisson pdf");
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res1.begin(),
                 callcdf<boost::math::poisson >(poissond));
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res2.begin(),
                 boost::lambda::bind(ppois, boost::lambda::_1, .2, 1, 0));
  compare(res1, res2, "poisson cdf");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::poisson >(poissond));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qpois, boost::lambda::_1, .2, 1, 0));
  compare(res1, res2, "poisson quantile");
  
  boost::math::negative_binomial nbinomd(5, 0.2);
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res1.begin(),
                 callpdf<boost::math::negative_binomial>(nbinomd));
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res2.begin(),
                 boost::lambda::bind(dnbinom, boost::lambda::_1, 5, .2, 0));
  compare(res1, res2, "nbinom pdf");
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res1.begin(),
                 callcdf<boost::math::negative_binomial>(nbinomd));
  std::transform(boost::counting_iterator<int>(0), boost::counting_iterator<int>(N+1), res2.begin(),
                 boost::lambda::bind(pnbinom, boost::lambda::_1, 5, .2, 1, 0));
  compare(res1, res2, "nbinom cdf");
  std::transform(qs.begin(), qs.end(), res1.begin(), callquantile<boost::math::negative_binomial>(nbinomd));
  std::transform(qs.begin(), qs.end(), res2.begin(), boost::lambda::bind(qnbinom, boost::lambda::_1, 5, .2, 1, 0));
  compare(res1, res2, "nbinom quantile");

  res2 = res1 = xs;
  double (*pf_d_d)(double) = &boost::math::lgamma;
  std::transform(xs.begin(), xs.end(), res1.begin(), boost::lambda::bind(pf_d_d, boost::lambda::_1));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(lgammafn, boost::lambda::_1));
  compare(res1, res2, "lgamma");

  pf_d_d = &boost::math::tgamma;
  std::transform(xs.begin(), xs.end(), res1.begin(), boost::lambda::bind(pf_d_d, boost::lambda::_1));
  std::transform(xs.begin(), xs.end(), res2.begin(), boost::lambda::bind(gammafn, boost::lambda::_1));
  compare(res1, res2, "tgamma");


  return 0;
}