/*
 * Copyright Nick Thompson, 2024
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0. (See accompanying file
 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 */
#ifndef BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP
#define BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP
#include <algorithm> // for std::sort
#include <cmath>
#include <limits>
#include <sstream>
#include <stdexcept>
#include <random>
#include <type_traits>  // for std::false_type

namespace boost::math::optimization::detail {

template <typename T, typename = void> struct has_resize : std::false_type {};

template <typename T>
struct has_resize<T, std::void_t<decltype(std::declval<T>().resize(size_t{}))>> : std::true_type {};

template <typename T> constexpr bool has_resize_v = has_resize<T>::value;

template <typename ArgumentContainer>
void validate_bounds(ArgumentContainer const &lower_bounds, ArgumentContainer const &upper_bounds) {
  using std::isfinite;
  std::ostringstream oss;
  if (lower_bounds.size() == 0) {
    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
    oss << ": The dimension of the problem cannot be zero.";
    throw std::domain_error(oss.str());
  }
  if (upper_bounds.size() != lower_bounds.size()) {
    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
    oss << ": There must be the same number of lower bounds as upper bounds, but given ";
    oss << upper_bounds.size() << " upper bounds, and " << lower_bounds.size() << " lower bounds.";
    throw std::domain_error(oss.str());
  }
  for (size_t i = 0; i < lower_bounds.size(); ++i) {
    auto lb = lower_bounds[i];
    auto ub = upper_bounds[i];
    if (lb > ub) {
      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
      oss << ": The upper bound must be greater than or equal to the lower bound, but the upper bound is " << ub
          << " and the lower is " << lb << ".";
      throw std::domain_error(oss.str());
    }
    if (!isfinite(lb)) {
      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
      oss << ": The lower bound must be finite, but got " << lb << ".";
      oss << " For infinite bounds, emulate with std::numeric_limits<Real>::lower() or use a standard infinite->finite "
             "transform.";
      throw std::domain_error(oss.str());
    }
    if (!isfinite(ub)) {
      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
      oss << ": The upper bound must be finite, but got " << ub << ".";
      oss << " For infinite bounds, emulate with std::numeric_limits<Real>::max() or use a standard infinite->finite "
             "transform.";
      throw std::domain_error(oss.str());
    }
  }
}

template <typename ArgumentContainer, class URBG>
std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const &lower_bounds,
                                                         ArgumentContainer const &upper_bounds,
                                                         size_t initial_population_size, URBG &&gen) {
  using Real = typename ArgumentContainer::value_type;
  using DimensionlessReal = decltype(Real()/Real());
  constexpr bool has_resize = detail::has_resize_v<ArgumentContainer>;
  std::vector<ArgumentContainer> population(initial_population_size);
  auto const dimension = lower_bounds.size();
  for (size_t i = 0; i < population.size(); ++i) {
    if constexpr (has_resize) {
      population[i].resize(dimension);
    } else {
      // Argument type must be known at compile-time; like std::array:
      if (population[i].size() != dimension) {
        std::ostringstream oss;
        oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
        oss << ": For containers which do not have resize, the default size must be the same as the dimension, ";
        oss << "but the default container size is " << population[i].size() << " and the dimension of the problem is "
            << dimension << ".";
        oss << " The function argument container type is " << typeid(ArgumentContainer).name() << ".\n";
        throw std::runtime_error(oss.str());
      }
    }
  }

  // Why don't we provide an option to initialize with (say) a Gaussian distribution?
  // > If the optimum's location is fairly well known,
  // > a Gaussian distribution may prove somewhat faster, although it
  // > may also increase the probability that the population will converge prematurely.
  // > In general, uniform distributions are preferred, since they best reflect
  // > the lack of knowledge about the optimum's location.
  //  - Differential Evolution: A Practical Approach to Global Optimization
  // That said, scipy uses Latin Hypercube sampling and says self-avoiding sequences are preferable.
  // So this is something that could be investigated and potentially improved.
  using std::uniform_real_distribution;
  uniform_real_distribution<DimensionlessReal> dis(DimensionlessReal(0), DimensionlessReal(1));
  for (size_t i = 0; i < population.size(); ++i) {
    for (size_t j = 0; j < dimension; ++j) {
      auto const &lb = lower_bounds[j];
      auto const &ub = upper_bounds[j];
      population[i][j] = lb + dis(gen) * (ub - lb);
    }
  }

  return population;
}

template <typename ArgumentContainer>
void validate_initial_guess(ArgumentContainer const &initial_guess, ArgumentContainer const &lower_bounds,
                            ArgumentContainer const &upper_bounds) {
  using std::isfinite;
  std::ostringstream oss;
  auto const dimension = lower_bounds.size();
  if (initial_guess.size() != dimension) {
    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
    oss << ": The initial guess must have the same dimensions as the problem,";
    oss << ", but the problem size is " << dimension << " and the initial guess has " << initial_guess.size()
        << " elements.";
    throw std::domain_error(oss.str());
  }
  for (size_t i = 0; i < dimension; ++i) {
    auto lb = lower_bounds[i];
    auto ub = upper_bounds[i];
    if (!isfinite(initial_guess[i])) {
      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
      oss << ": At index " << i << ", the initial guess is " << initial_guess[i]
          << ", make sure all elements of the initial guess are finite.";
      throw std::domain_error(oss.str());
    }
    if (initial_guess[i] < lb || initial_guess[i] > ub) {
      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
      oss << ": At index " << i << " the initial guess " << initial_guess[i] << " is not in the bounds [" << lb << ", "
          << ub << "].";
      throw std::domain_error(oss.str());
    }
  }
}

// Return indices corresponding to the minimum function values.
template <typename Real> std::vector<size_t> best_indices(std::vector<Real> const &function_values) {
  using std::isnan;
  const size_t n = function_values.size();
  std::vector<size_t> indices(n);
  for (size_t i = 0; i < n; ++i) {
    indices[i] = i;
  }

  std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) {
    if (isnan(function_values[a])) {
      return false;
    }
    if (isnan(function_values[b])) {
      return true;
    }
    return function_values[a] < function_values[b];
  });
  return indices;
}

template<typename RandomAccessContainer>
auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessContainer const & weights) {
  using std::isfinite;
  if (values.size() != weights.size()) {
    std::ostringstream oss;
    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
    oss << ": There must be the same number of weights as values, but got " << values.size() << " values and " << weights.size() << " weights.";
    throw std::logic_error(oss.str());
  }
  if (values.size() == 0) {
    std::ostringstream oss;
    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
    oss << ": There must at least one value provided.";
    throw std::logic_error(oss.str());
  }
  using Real = typename RandomAccessContainer::value_type;
  Real numerator = 0;
  Real denominator = 0;
  for (size_t i = 0; i < values.size(); ++i) {
    if (weights[i] < 0 || !isfinite(weights[i])) {
      std::ostringstream oss;
      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
      oss << ": All weights must be positive and finite, but got received weight " << weights[i] << " at index " << i << " of " << weights.size() << ".";
      throw std::domain_error(oss.str());
    }
    Real tmp = weights[i]*values[i];
    numerator += tmp*values[i];
    denominator += tmp;
  }
  return numerator/denominator;
}

} // namespace boost::math::optimization::detail
#endif