// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "histc.h"
#include <cassert>
#include <iostream>
#include <omp.h>

template <
    typename DerivedX, typename DerivedE, typename DerivedN, typename DerivedB>
inline void igl::histc(
    const Eigen::MatrixBase<DerivedX>& X, const Eigen::MatrixBase<DerivedE>& E,
    Eigen::PlainObjectBase<DerivedN>& N, Eigen::PlainObjectBase<DerivedB>& B) {
  histc(X, E, B);
  const std::size_t n = E.size();
  const std::size_t m = X.size();
  assert(m == B.size());
  N.resize(n, 1);
  N.setConstant(0);
#pragma omp parallel for num_threads(8)
  for (std::size_t j = 0; j < m; ++j) {
    if (B(j) >= 0) {
#pragma omp atomic
      N(B(j))++;
    }
  }
}

template <typename DerivedX, typename DerivedE, typename DerivedB>
inline void igl::histc(
    const Eigen::MatrixBase<DerivedX>& X, const Eigen::MatrixBase<DerivedE>& E,
    Eigen::PlainObjectBase<DerivedB>& B) {
  const std::size_t m = X.size();
  using namespace std;
  assert(
      (E.bottomRightCorner(E.size() - 1, 1) - E.topLeftCorner(E.size() - 1, 1))
              .maxCoeff() >= 0 &&
      "E should be monotonically increasing");
  B.resize(m, 1);
#pragma omp parallel for num_threads(8)
  for (std::size_t j = 0; j < m; ++j) {
    const double x = X(j);
    // Boring one-offs
    if (x < E(0) || x > E(E.size() - 1)) {
      B(j) = -1;
      continue;
    }
    // Find x in E
    std::size_t l = 0;
    std::size_t h = E.size() - 1;
    std::size_t k = l;
    while ((h - l) > 1) {
      assert(x >= E(l));
      assert(x <= E(h));
      k = (h + l) / 2;
      if (x < E(k)) {
        h = k;
      } else {
        l = k;
      }
    }
    if (x == E(h)) {
      k = h;
    } else {
      k = l;
    }
    B(j) = k;
  }
}

template <typename DerivedE>
inline void igl::histc(
    const typename DerivedE::Scalar& x, const Eigen::MatrixBase<DerivedE>& E,
    typename DerivedE::Index& b) {
  Eigen::Matrix<typename DerivedE::Scalar, 1, 1> X;
  X(0) = x;
  Eigen::Matrix<typename DerivedE::Index, 1, 1> B;
  hist(X, E, B);
  b = B(0);
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::histc<
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1> >(
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&,
    Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
template void igl::histc<
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>,
    Eigen::Matrix<int, -1, -1, 0, -1, -1> >(
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
template void igl::histc<
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>,
    Eigen::Matrix<int, -1, 1, 0, -1, 1> >(
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
template void igl::histc<
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<int, -1, 1, 0, -1, 1> >(
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
template void igl::histc<
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<double, -1, 1, 0, -1, 1>,
    Eigen::Matrix<int, -1, -1, 0, -1, -1> >(
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
#else
template void igl::histc<
    Eigen::CwiseNullaryOp<
        Eigen::internal::linspaced_op<int, true>,
        Eigen::Matrix<int, -1, 1, 0, -1, 1> >,
    Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>,
    Eigen::Matrix<int, -1, 1, 0, -1, 1> >(
    Eigen::MatrixBase<Eigen::CwiseNullaryOp<
        Eigen::internal::linspaced_op<int, true>,
        Eigen::Matrix<int, -1, 1, 0, -1, 1> > > const&,
    Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&,
    Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
#endif

#endif
