#include <scitbx/array_family/boost_python/flex_fwd.h>

#include <scitbx/math/parabolic_cylinder_d.h>
#include <scitbx/math/bessel.h>
#include <scitbx/math/chebyshev.h>
#include <scitbx/math/dihedral.h>
#include <scitbx/math/erf.h>
#include <scitbx/math/euler_angles.h>
#include <scitbx/math/floating_point_epsilon.h>
#include <scitbx/math/gamma.h>
#include <scitbx/math/gcd.h>
#include <scitbx/math/halton.h>
#include <scitbx/math/lambertw.h>
#include <scitbx/math/phase_error.h>
#include <scitbx/math/resample.h>
#include <scitbx/math/superpose.h>
#include <scitbx/math/utils.h>
#include <boost/rational.hpp> // for boost::gcd
#include <scitbx/math/approx_equal.h>
#include <scitbx/math/orthonormal_basis.h>
#include <scitbx/math/gaussian_fit_1d_analytical.h>
#include <scitbx/math/cubic_equation.h>
#include <scitbx/math/distance_difference.h>
#include <boost/python/module.hpp>
#include <boost/python/def.hpp>
#include <boost/python/class.hpp>
#include <boost/python/args.hpp>
#include <scitbx/math/fast_approx_math.h>
#include <scitbx/math/g_function.h>


namespace scitbx { namespace math {
namespace boost_python {

  void wrap_basic_statistics();
  void wrap_gaussian();
  void wrap_golay();
  void wrap_minimum_covering_sphere();
  void wrap_principal_axes_of_inertia();
  void wrap_row_echelon();
  void wrap_tensor_rank_2();
  void wrap_icosahedron();
  void wrap_chebyshev_base();
  void wrap_chebyshev_polynome();
  void wrap_chebyshev_fitter();
  void wrap_chebyshev_lsq();
  void wrap_slatec();
  void wrap_line_search();
  void wrap_r3_rotation();
  void wrap_resample();
  void wrap_quadrature();
  void wrap_unimodular_generator();
  void wrap_halton();
  void wrap_least_squares_plane();
  void wrap_continued_fraction();
  void wrap_numeric_limits();
  void wrap_distributions();
  void wrap_exp_functions();
  void wrap_zernike();
  void wrap_zernike_mom();
  void wrap_2d_zernike_mom();
  void wrap_weighted_covariance();
  void wrap_dmatrix();
  void wrap_correlation();
  void wrap_interpolation();
  void wrap_tetrahedron();
  void wrap_angle_derivative();

namespace {

  int
  time_gcd_int_boost(
    int n)
  {
    int result = 0;
    for(int a=0;a<n;a++) {
      for(int b=0;b<n;b++) {
        int c = boost::gcd(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

  long
  time_gcd_long_boost(
    long n)
  {
    long result = 0;
    for(long a=0;a<n;a++) {
      for(long b=0;b<n;b++) {
        long c = boost::gcd(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

  int
  time_gcd_int_simple(
    int n)
  {
    int result = 0;
    for(int a=0;a<n;a++) {
      for(int b=0;b<n;b++) {
        int c = gcd_int_simple(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

  long
  time_gcd_long_simple(
    long n)
  {
    long result = 0;
    for(long a=0;a<n;a++) {
      for(long b=0;b<n;b++) {
        long c = gcd_long_simple(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

  long
  time_gcd_unsigned_long_binary(
    unsigned long n)
  {
    unsigned long result = 0;
    for(unsigned long a=0;a<n;a++) {
      for(unsigned long b=0;b<n;b++) {
        unsigned long c = gcd_unsigned_long_binary(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

  long
  time_gcd_long_binary(
    long n)
  {
    long result = 0;
    for(long a=0;a<n;a++) {
      for(long b=0;b<n;b++) {
        long c = gcd_long_binary(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

#if defined(SCITBX_MATH_GCD_USING_ASM)
  int
  time_gcd_int32_asm(
    int n)
  {
    int result = 0;
    for(int a=0;a<n;a++) {
      for(int b=0;b<n;b++) {
        int c = gcd_int32_asm(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }

# if defined(__x86_64__)
  long
  time_gcd_int64_asm(
    long n)
  {
    long result = 0;
    for(long a=0;a<n;a++) {
      for(long b=0;b<n;b++) {
        long c = gcd_int64_asm(a, b);
        if (result < c) result = c;
      }
    }
    return result;
  }
# endif
#endif

  template <typename SitesType>
  boost::optional<double>
  dihedral_angle(
    SitesType const& sites,
    bool deg)
  {
    return dihedral(sites).angle(deg);
  }

  mat3<double>
  superpose_kearsley_rotation(
    af::const_ref<vec3<double> > const& reference_sites,
    af::const_ref<vec3<double> > const& other_sites)
  {
    return superpose::superposition<>::kearsley_rotation(
      reference_sites, other_sites);
  }

  mat3< double >
  euler_angles_xyz_matrix(
      const double& ax,
      const double& ay,
      const double& az )
  {
    return euler_angles::xyz_matrix( ax, ay, az );
  }

  vec3< double >
  euler_angles_xyz_angles(
      const mat3< double >& m,
      const double& eps = 1e-12 )
  {
    return euler_angles::xyz_angles( m, eps );
  }

  mat3< double >
  euler_angles_yzx_matrix(
      const double& ay,
      const double& az,
      const double& ax )
  {
    return euler_angles::yzx_matrix( ay, az, ax );
  }

  vec3< double >
  euler_angles_yzx_angles(
      const mat3< double >& m,
      const double& eps = 1e-12 )
  {
    return euler_angles::yzx_angles( m, eps );
  }

  mat3< double >
  euler_angles_zyz_matrix(
      const double& az1,
      const double& ay,
      const double& az3 )
  {
    return euler_angles::zyz_matrix( az1, ay, az3 );
  }

  vec3< double >
  euler_angles_zyz_angles(
      const mat3< double >& m,
      const double& eps = 1e-12 )
  {
    return euler_angles::zyz_angles( m, eps );
  }

  template <typename T>
  struct approx_equal_relatively_wrapper
  {
    typedef math::approx_equal_relatively<T> wt;
    typedef T arg_t;
    typedef typename wt::amplitude_type amplitude_t;
    static bool form_1(arg_t x, arg_t y, amplitude_t relative_error) {
      wt p(relative_error);
      return p(x, y);
    }
    static bool form_2(arg_t x, arg_t y, amplitude_t relative_error,
                       amplitude_t near_zero_threshold)
    {
      wt p(relative_error, near_zero_threshold);
      return p(x, y);
    }

    static void wrap() {
      using namespace boost::python;
      def("approx_equal_relatively", form_1,
          (arg("x"), arg("y"), arg("relative_error")));
      def("approx_equal_relatively", form_2,
          (arg("x"), arg("y"), arg("relative_error"),
           arg("near_zero_threshold")));
    }
  };

  void init_module()
  {
    using namespace boost::python;

    def("time_gcd_int_boost", time_gcd_int_boost);
    def("time_gcd_long_boost", time_gcd_long_boost);
    def("gcd_int_simple", gcd_int_simple, (arg("a"), arg("b")));
    def("time_gcd_int_simple", time_gcd_int_simple);
    def("gcd_long_simple", gcd_long_simple, (arg("a"), arg("b")));
    def("time_gcd_long_simple", time_gcd_long_simple);
    def("time_gcd_unsigned_long_binary", time_gcd_unsigned_long_binary);
    def("gcd_long_binary", gcd_long_binary, (arg("a"), arg("b")));
    def("time_gcd_long_binary", time_gcd_long_binary);
#if defined(SCITBX_MATH_GCD_USING_ASM)
    def("gcd_int32_asm", gcd_int32_asm, (arg("a"), arg("b")));
    def("time_gcd_int32_asm", time_gcd_int32_asm);
# if defined(__x86_64__)
    def("gcd_int64_asm", gcd_int64_asm, (arg("a"), arg("b")));
    def("time_gcd_int64_asm", time_gcd_int64_asm);
# endif
#endif

    def("floating_point_epsilon_float_get",
      &floating_point_epsilon<float>::get);
    def("floating_point_epsilon_double_get",
      &floating_point_epsilon<double>::get);

    def("erf", (double(*)(double const&)) erf);
    def("erf",
      (scitbx::af::shared<double>(*)(
        scitbx::af::const_ref<double> const&)) erf);
    def("erfc", (double(*)(double const&)) erfc);
    def("erfcx", (double(*)(double const&)) erfcx);

    def("parabolic_cylinder_d", (double(*)(double, double))
      parabolic_cylinder_d::dv);

    // G-function
    def("GfuncOfRSsqr_approx", (double(*)(double))
      g_function::GfuncOfRSsqr_approx);

    def("bessel_i1_over_i0", (double(*)(double const&)) bessel::i1_over_i0);
    def("bessel_i1_over_i0",
      (scitbx::af::shared<double>(*)(scitbx::af::const_ref<double> const&))
        bessel::i1_over_i0);
    def("bessel_inverse_i1_over_i0",
      (double(*)(double const&)) bessel::inverse_i1_over_i0);
    def("inverse_bessel_i1_over_i0", (scitbx::af::shared<double>(*)(
         scitbx::af::const_ref<double> const&)) bessel::inverse_i1_over_i0);
    def("bessel_i0", (double(*)(double const&)) bessel::i0);
    def("bessel_i1", (double(*)(double const&)) bessel::i1);
    def("bessel_ln_of_i0", (double(*)(double const&)) bessel::ln_of_i0);
    def("ei1", (double(*)(double const&)) bessel::ei1);
    def("ei0", (double(*)(double const&)) bessel::ei0);

    typedef return_value_policy<return_by_value> rbv;
    namespace smg=scitbx::math::gaussian_fit_1d_analytical;
    class_<smg::compute<> >("gaussian_fit_1d_analytical")
      .def(init<
           af::const_ref<double> const&,
           af::const_ref<double> const&,
           af::const_ref<double> const& >(
             (arg("x"),
              arg("y"),
              arg("z"))))
      .def(init<
           af::const_ref<double> const&,
           af::const_ref<double> const& >(
             (arg("x"),
              arg("y"))))
      .add_property("a", make_getter(&smg::compute<>::a, rbv()))
      .add_property("b", make_getter(&smg::compute<>::b, rbv()))
    ;

    //typedef return_value_policy<return_by_value> rbv;
    namespace cueq=scitbx::math::cubic_equation;
    class_<cueq::real<> >("cubic_equation_real")
      .def(init<
           double const&,
           double const&,
           double const&,
           double const& >(
             (arg("a"),
              arg("b"),
              arg("c"),
              arg("d"))))
      .def("residual", &cueq::real<>::residual)
      .add_property("x", make_getter(&cueq::real<>::x, rbv()))
      .add_property("A", make_getter(&cueq::real<>::A,  rbv()))
      .add_property("B", make_getter(&cueq::real<>::B,  rbv()))
      .add_property("D", make_getter(&cueq::real<>::D,  rbv()))
    ;

#if defined(SCITBX_MATH_BESSEL_HAS_SPHERICAL)
    def("spherical_bessel",
      (double(*)(int const&, double const&))
        bessel::spherical_bessel);
    def("spherical_bessel_array",
      (scitbx::af::shared< double> (*)(
        int const&, scitbx::af::shared<double> const&))
          bessel::spherical_bessel_array);
    def("bessel_J",
      (double(*)(int const&, double const&))
        bessel::bessel_J);
    def("bessel_J_array",
      (scitbx::af::shared< double> (*)(
        int const&, scitbx::af::shared<double> const&))
          bessel::bessel_J_array);
    def ("bessel_J_zeroes",
           (scitbx::af::shared< double> (*)(
        double const&, int const&))
        bessel::bessel_J_zeroes);
    def ("sph_bessel_j_zeroes",
           (scitbx::af::shared< double> (*)(
        double const&, int const&))
        bessel::sph_bessel_j_zeroes);


#endif

    def("gamma_complete", (double(*)(double const&, bool))
      gamma::complete, (
        arg("x"),
        arg("minimax")=true));
    def("gamma_incomplete", (double(*)(double const&,
                                       double const&,
                                       unsigned))
      gamma::incomplete, (arg("a"),
        arg("x"),
        arg("max_iterations")=500));
    def("gamma_incomplete_complement",(double(*)(double const&,
                                                 double const&,
                                                 unsigned))
      gamma::incomplete_complement, (
        arg("a"),
        arg("x"),
        arg("max_iterations")=500));
    def("exponential_integral_e1z", (double(*)(double const&))
         gamma::exponential_integral_e1z );

    def("lambertw", (double(*)(double const&, unsigned))
      lambertw, (
        arg("x"),
        arg("max_iterations")=100));

    wrap_basic_statistics();
    wrap_gaussian();
    wrap_golay();
    wrap_minimum_covering_sphere();
    wrap_principal_axes_of_inertia();
    wrap_row_echelon();
    wrap_tensor_rank_2();
    wrap_icosahedron();
    wrap_chebyshev_base();
    wrap_chebyshev_polynome();
    wrap_chebyshev_fitter();
    wrap_chebyshev_lsq();
    wrap_slatec();
    wrap_line_search();
    wrap_r3_rotation();
    wrap_resample();
    wrap_quadrature();
    wrap_unimodular_generator();
    wrap_halton();
    wrap_least_squares_plane();
    wrap_continued_fraction();
    wrap_numeric_limits();
    wrap_distributions();
    wrap_exp_functions();
    wrap_zernike();
    wrap_zernike_mom();
    wrap_2d_zernike_mom();
    wrap_weighted_covariance();
    wrap_dmatrix();
    wrap_correlation();
    wrap_interpolation();
    wrap_tetrahedron();
    wrap_angle_derivative();

    def("superpose_kearsley_rotation", superpose_kearsley_rotation, (
      arg("reference_sites"), arg("other_sites")));

    def("dihedral_angle",
      dihedral_angle<af::tiny<vec3<double>, 4> >, (
        arg("sites"), arg("deg")=false));
    def("dihedral_angle",
      dihedral_angle<af::const_ref<vec3<double> > >, (
        arg("sites"), arg("deg")=false));

    def( "euler_angles_xyz_matrix", euler_angles_xyz_matrix, (
      arg("ax"), arg("ay"), arg("az")));
    def( "euler_angles_xyz_angles", euler_angles_xyz_angles, (
      arg("m"), arg("eps")=1e-12));

    def( "euler_angles_yzx_matrix", euler_angles_yzx_matrix, (
      arg("ay"), arg("az"), arg("ax")));
    def( "euler_angles_yzx_angles", euler_angles_yzx_angles, (
      arg("m"), arg("eps")=1e-12));

    def( "euler_angles_zyz_matrix", euler_angles_zyz_matrix, (
      arg("az1"), arg("ay"), arg("az3")));
    def( "euler_angles_zyz_angles", euler_angles_zyz_angles, (
      arg("m"), arg("eps")=1e-12));

    def("approx_sqrt", math::approx_sqrt);

    def("cos_table",
      (double(*)(
         af::const_ref<double> const&,
         double,
         double const& ,
         int const& ,
         bool))
           math::cos_table, (
             arg("table"), arg("arg"), arg("step"), arg("n"),
             arg("interpolate")));

    def("sin_table",
      (double(*)(
         af::const_ref<double> const&,
         double,
         double const& ,
         int const& ,
         bool))
           math::sin_table, (
             arg("table"), arg("arg"), arg("step"), arg("n"),
             arg("interpolate")));

    def("signed_phase_error",
      (double(*)(
        double const&, double const&, bool))
          math::signed_phase_error, (
            arg("phi1"), arg("phi2"), arg("deg")=false));
    def("signed_phase_error",
      (af::shared<double>(*)(
        af::const_ref<double> const&, af::const_ref<double> const&, bool))
          math::signed_phase_error, (
            arg("phi1"), arg("phi2"), arg("deg")=false));
    def("phase_error",
      (double(*)(
        double const&, double const&, bool))
          math::phase_error, (
            arg("phi1"), arg("phi2"), arg("deg")=false));
    def("phase_error",
      (af::shared<double>(*)(
        af::const_ref<double> const&, af::const_ref<double> const&, bool))
          math::phase_error, (
            arg("phi1"), arg("phi2"), arg("deg")=false));
    def("nearest_phase",
      (double(*)(
        double const&, double const&, bool))
          math::nearest_phase, (
            arg("reference"), arg("other"), arg("deg")=false));
    def("nearest_phase",
      (af::shared<double>(*)(
        af::const_ref<double> const&, af::const_ref<double> const&, bool))
          math::nearest_phase, (
            arg("reference"), arg("other"), arg("deg")=false));
    def("divmod", math::divmod);
    approx_equal_relatively_wrapper<double>::wrap();
    approx_equal_relatively_wrapper<std::complex<double> >::wrap();
    {
      af::tiny<vec3<double>, 3> (*f1)(vec3<double> const &, vec3<double> const &,
                                      bool) = &orthonormal_basis;
      af::tiny<vec3<double>, 3> (*f2)(vec3<double> const &, int,
                                      vec3<double> const &, int,
                                      bool) = &orthonormal_basis;
      def("orthonormal_basis", f1, (arg("v0"), arg("v1"),
                                    arg("right_handed")=true));
      def("orthonormal_basis", f2, (arg("v0"), arg("axis_index_1"),
                                    arg("v1"), arg("axis_index_2"),
                                    arg("right_handed")=true));
    }
    def("distance_difference_matrix", distance_difference_matrix<double>, (
      arg("sites1"), arg("sites2")));
  }

}}}} // namespace scitbx::math::boost_python::<anonymous>

BOOST_PYTHON_MODULE(scitbx_math_ext)
{
  scitbx::math::boost_python::init_module();
}