<commit_before>
// Vc. Must be before the ROOT includes for std:: math functions to work...
#include <Vc/Vc>

// ROOT
#include "Math/GenVector/PositionVector3D.h"
#include "Math/GenVector/DisplacementVector3D.h"
#include "Math/GenVector/Plane3D.h"
#include "Math/GenVector/Transform3D.h"
#include "TStopwatch.h"

// STL
#include <random>
#include <vector>
#include <iostream>
#include <string>
#include <typeinfo>
#include <cmath>

// note scale here is > 1 as SIMD and scalar floating point calculations not
// expected to be bit wise identical
int compare(double v1, double v2, const std::string &name = "", double scale = 1000.0)
{
   //  ntest = ntest + 1;

   // numerical double limit for epsilon
   const double eps     = scale * std::numeric_limits<double>::epsilon();
   int          iret    = 0;
   double       delta   = v2 - v1;
   double       d       = 0;
   if (delta < 0) delta = -delta;
   if (v1 == 0 || v2 == 0) {
      if (delta > eps) {
         iret = 1;
      }
   }
   // skip case v1 or v2 is infinity
   else {
      d = v1;

      if (v1 < 0) d = -d;
      // add also case when delta is small by default
      if (delta / d > eps && delta > eps) iret = 1;
   }

   if (iret != 0) {
      int pr = std::cout.precision(18);
      std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << int(delta / d / eps)
                << "   (Allowed discrepancy is " << eps << ")\n";
      std::cout.precision(pr);
   }
   return iret;
}

// randomn generator
static std::default_random_engine gen;
// Distributions for each member
static std::uniform_real_distribution<double> p_x(-800, 800), p_y(-600, 600), p_z(10000, 10500);
static std::uniform_real_distribution<double> d_x(-0.2, 0.2), d_y(-0.1, 0.1), d_z(0.95, 0.99);
static std::uniform_real_distribution<double> c_x(3100, 3200), c_y(10, 15), c_z(3200, 3300);
static std::uniform_real_distribution<double> r_rad(8500, 8600);
static std::uniform_real_distribution<double> p0(-0.002, 0.002), p1(-0.2, 0.2), p2(0.97, 0.99), p3(-1300, 1300);

template <typename POINT, typename VECTOR, typename PLANE, typename FTYPE>
class Data {
public:
   typedef std::vector<Data, Vc::Allocator<Data>> Vector;

public:
   POINT  position;
   VECTOR direction;
   POINT  CoC;
   PLANE  plane;
   FTYPE  radius{0};

public:
   template <typename INDATA>
   Data(const INDATA &ind)
      : position(ind.position.x(), ind.position.y(), ind.position.z()),
        direction(ind.direction.x(), ind.direction.y(), ind.direction.z()), CoC(ind.CoC.x(), ind.CoC.y(), ind.CoC.z()),
        plane(ind.plane.A(), ind.plane.B(), ind.plane.C(), ind.plane.D()), radius(ind.radius)
   {
   }
   Data()
      : position(p_x(gen), p_y(gen), p_z(gen)), direction(d_x(gen), d_y(gen), d_z(gen)),
        CoC(c_x(gen), c_y(gen), c_z(gen)), plane(p0(gen), p1(gen), p2(gen), p3(gen)), radius(r_rad(gen))
   {
   }
};

template <typename INDATA, typename OUTDATA>
void fill(const INDATA &in, OUTDATA &out)
{
   out.clear();
   out.reserve(in.size());
   for (const auto &i : in) {
      out.emplace_back(i);
   }
}

template <typename POINT, typename VECTOR, typename FTYPE>
inline
   typename std::enable_if<std::is_arithmetic<typename POINT::Scalar>::value &&
                              std::is_arithmetic<typename VECTOR::Scalar>::value && std::is_arithmetic<FTYPE>::value,
                           bool>::type
   reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, const FTYPE radius)
{
   constexpr FTYPE zero(0), two(2.0), four(4.0), half(0.5);
   const FTYPE     a     = direction.Mag2();
   const VECTOR    delta = position - CoC;
   const FTYPE     b     = two * direction.Dot(delta);
   const FTYPE     c     = delta.Mag2() - radius * radius;
   const FTYPE     discr = b * b - four * a * c;
   const bool      OK    = discr > zero;
   if (OK) {
      const FTYPE dist = half * (std::sqrt(discr) - b) / a;
      // change position to the intersection point
      position += dist * direction;
      // reflect the vector
      // r = u - 2(u.n)n, r=reflection, u=incident, n=normal
      const VECTOR normal = position - CoC;
      direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal;
   }
   return OK;
}

template <typename POINT, typename VECTOR, typename FTYPE>
inline
   typename std::enable_if<!std::is_arithmetic<typename POINT::Scalar>::value &&
                              !std::is_arithmetic<typename VECTOR::Scalar>::value && !std::is_arithmetic<FTYPE>::value,
                           typename FTYPE::mask_type>::type
   reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, const FTYPE radius)
{
   const FTYPE               two(2.0), four(4.0), half(0.5);
   const FTYPE               a     = direction.Mag2();
   const VECTOR              delta = position - CoC;
   const FTYPE               b     = two * direction.Dot(delta);
   const FTYPE               c     = delta.Mag2() - radius * radius;
   FTYPE                     discr = b * b - four * a * c;
   typename FTYPE::mask_type OK    = discr > FTYPE::Zero();
   if (any_of(OK)) {
      // Zero out the negative values in discr, to prevent sqrt(-ve)
      discr(!OK) = FTYPE::Zero();
      // compute the distance
      const FTYPE dist = half * (sqrt(discr) - b) / a;
      // change position to the intersection point
      position += dist * direction;
      // reflect the vector
      // r = u - 2(u.n)n, r=reflection, u=incident, n=normal
      const VECTOR normal = position - CoC;
      direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal;
   }
   // return the mask indicating which results should be used
   return OK;
}

template <typename POINT, typename VECTOR, typename PLANE>
inline typename std::enable_if<std::is_arithmetic<typename POINT::Scalar>::value &&
                                  std::is_arithmetic<typename VECTOR::Scalar>::value &&
                                  std::is_arithmetic<typename PLANE::Scalar>::value,
                               bool>::type
reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane)
{
   constexpr typename POINT::Scalar two(2.0);
   const bool                       OK = true;
   // Plane normal
   const auto &normal = plane.Normal();
   // compute distance to the plane
   const auto scalar   = direction.Dot(normal);
   const auto distance = -(plane.Distance(position)) / scalar;
   // change position to reflection point and update direction
   position += distance * direction;
   direction -= two * scalar * normal;
   return OK;
}

template <typename POINT, typename VECTOR, typename PLANE, typename FTYPE = typename POINT::Scalar>
inline typename std::enable_if<!std::is_arithmetic<typename POINT::Scalar>::value &&
                                  !std::is_arithmetic<typename VECTOR::Scalar>::value &&
                                  !std::is_arithmetic<typename PLANE::Scalar>::value,
                               typename FTYPE::mask_type>::type
reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane)
{
   const typename POINT::Scalar    two(2.0);
   const typename FTYPE::mask_type OK(true);
   // Plane normal<commit_msg>Update Enable If pattern in GenVector Vc test (#431)<commit_after>
// Vc. Must be before the ROOT includes for std:: math functions to work...
#include <Vc/Vc>

// ROOT
#include "Math/GenVector/PositionVector3D.h"
#include "Math/GenVector/DisplacementVector3D.h"
#include "Math/GenVector/Plane3D.h"
#include "Math/GenVector/Transform3D.h"
#include "TStopwatch.h"

// STL
#include <random>
#include <vector>
#include <iostream>
#include <string>
#include <typeinfo>
#include <cmath>
#include <type_traits>

// note scale here is > 1 as SIMD and scalar floating point calculations not
// expected to be bit wise identical
int compare(double v1, double v2, const std::string &name = "", double scale = 1000.0)
{
   //  ntest = ntest + 1;

   // numerical double limit for epsilon
   const double eps     = scale * std::numeric_limits<double>::epsilon();
   int          iret    = 0;
   double       delta   = v2 - v1;
   double       d       = 0;
   if (delta < 0) delta = -delta;
   if (v1 == 0 || v2 == 0) {
      if (delta > eps) {
         iret = 1;
      }
   }
   // skip case v1 or v2 is infinity
   else {
      d = v1;

      if (v1 < 0) d = -d;
      // add also case when delta is small by default
      if (delta / d > eps && delta > eps) iret = 1;
   }

   if (iret != 0) {
      int pr = std::cout.precision(18);
      std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << int(delta / d / eps)
                << "   (Allowed discrepancy is " << eps << ")\n";
      std::cout.precision(pr);
   }
   return iret;
}

// randomn generator
static std::default_random_engine gen;
// Distributions for each member
static std::uniform_real_distribution<double> p_x(-800, 800), p_y(-600, 600), p_z(10000, 10500);
static std::uniform_real_distribution<double> d_x(-0.2, 0.2), d_y(-0.1, 0.1), d_z(0.95, 0.99);
static std::uniform_real_distribution<double> c_x(3100, 3200), c_y(10, 15), c_z(3200, 3300);
static std::uniform_real_distribution<double> r_rad(8500, 8600);
static std::uniform_real_distribution<double> p0(-0.002, 0.002), p1(-0.2, 0.2), p2(0.97, 0.99), p3(-1300, 1300);

template <typename POINT, typename VECTOR, typename PLANE, typename FTYPE>
class Data {
public:
   typedef std::vector<Data, Vc::Allocator<Data>> Vector;

public:
   POINT  position;
   VECTOR direction;
   POINT  CoC;
   PLANE  plane;
   FTYPE  radius{0};

public:
   template <typename INDATA>
   Data(const INDATA &ind)
      : position(ind.position.x(), ind.position.y(), ind.position.z()),
        direction(ind.direction.x(), ind.direction.y(), ind.direction.z()), CoC(ind.CoC.x(), ind.CoC.y(), ind.CoC.z()),
        plane(ind.plane.A(), ind.plane.B(), ind.plane.C(), ind.plane.D()), radius(ind.radius)
   {
   }
   Data()
      : position(p_x(gen), p_y(gen), p_z(gen)), direction(d_x(gen), d_y(gen), d_z(gen)),
        CoC(c_x(gen), c_y(gen), c_z(gen)), plane(p0(gen), p1(gen), p2(gen), p3(gen)), radius(r_rad(gen))
   {
   }
};

template <typename INDATA, typename OUTDATA>
void fill(const INDATA &in, OUTDATA &out)
{
   out.clear();
   out.reserve(in.size());
   for (const auto &i : in) {
      out.emplace_back(i);
   }
}

template <typename POINT, typename VECTOR, typename FTYPE,
          typename = typename std::enable_if<std::is_arithmetic<typename POINT::Scalar>::value &&
                                             std::is_arithmetic<typename VECTOR::Scalar>::value &&
                                             std::is_arithmetic<FTYPE>::value>::type>
inline bool reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, const FTYPE radius)
{
   constexpr FTYPE zero(0), two(2.0), four(4.0), half(0.5);
   const FTYPE     a     = direction.Mag2();
   const VECTOR    delta = position - CoC;
   const FTYPE     b     = two * direction.Dot(delta);
   const FTYPE     c     = delta.Mag2() - radius * radius;
   const FTYPE     discr = b * b - four * a * c;
   const bool      OK    = discr > zero;
   if (OK) {
      const FTYPE dist = half * (std::sqrt(discr) - b) / a;
      // change position to the intersection point
      position += dist * direction;
      // reflect the vector
      // r = u - 2(u.n)n, r=reflection, u=incident, n=normal
      const VECTOR normal = position - CoC;
      direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal;
   }
   return OK;
}

template <typename POINT, typename VECTOR, typename FTYPE,
          typename = typename std::enable_if<!std::is_arithmetic<typename POINT::Scalar>::value &&
                                             !std::is_arithmetic<typename VECTOR::Scalar>::value &&
                                             !std::is_arithmetic<FTYPE>::value>::type>
inline typename FTYPE::mask_type reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC,
                                                  const FTYPE radius)
{
   const FTYPE               two(2.0), four(4.0), half(0.5);
   const FTYPE               a     = direction.Mag2();
   const VECTOR              delta = position - CoC;
   const FTYPE               b     = two * direction.Dot(delta);
   const FTYPE               c     = delta.Mag2() - radius * radius;
   FTYPE                     discr = b * b - four * a * c;
   typename FTYPE::mask_type OK    = discr > FTYPE::Zero();
   if (any_of(OK)) {
      // Zero out the negative values in discr, to prevent sqrt(-ve)
      discr(!OK) = FTYPE::Zero();
      // compute the distance
      const FTYPE dist = half * (sqrt(discr) - b) / a;
      // change position to the intersection point
      position += dist * direction;
      // reflect the vector
      // r = u - 2(u.n)n, r=reflection, u=incident, n=normal
      const VECTOR normal = position - CoC;
      direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal;
   }
   // return the mask indicating which results should be used
   return OK;
}

template <typename POINT, typename VECTOR, typename PLANE,
          typename = typename std::enable_if<std::is_arithmetic<typename POINT::Scalar>::value &&
                                             std::is_arithmetic<typename VECTOR::Scalar>::value>::type>
inline bool reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane)
{
   constexpr typename POINT::Scalar two(2.0);
   const bool                       OK = true;
   // Plane normal
   const auto &normal = plane.Normal();
   // compute distance to the plane
   const auto scalar   = direction.Dot(normal);
   const auto distance = -(plane.Distance(position)) / scalar;
   // change position to reflection point and update direction
   position += distance * direction;
   direction -= two * scalar * normal;
   return OK;
}

template <typename POINT, typename VECTOR, typename PLANE, typename FTYPE = typename POINT::Scalar,
          typename = typename std::enable_if<!std::is_arithmetic<typename POINT::Scalar>::value &&
                                             !std::is_arithmetic<typename VECTOR::Scalar>::value>::type>
inline typename FTYPE::mask_type reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane)
{
   const typename POINT::Scalar    two(2.0);
   const typename FTYPE::mask_type OK(true);
   // Plane normal
{'fraction_non_alphanumeric': 0.13976995468804462, 'fraction_numerical': 0.022237713489020566, 'mean_word_length': 2.605428499623021, 'pattern_counts': {'":': 0, '<': 126, '<?xml version=': 0, '>': 90, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5663997', 'n_tokens_mistral': 4865, 'n_tokens_neox': 4581, 'n_words': 1265}