// -*- C++ -*-
//
// Copyright Sylvain Bougerel 2009 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file COPYING or copy at
// http://www.boost.org/LICENSE_1_0.txt)

/**
 *  \file   spatial_math.hpp
 *  Gather in one file all the mathematical operations, mainly for the
 *  metric types.
 *
 *  Most of the mathematical operations written in this file are not very well
 *  optimized, and for a given compiler or architecture it would be easy to
 *  write more efficient algorithms. Therefore, if you are really looking to
 *  increase the speed of your computation, you might want to write your own
 *  metric optimized for the type you are working with.
 *
 *  \see neighbor
 */

#ifndef SPATIAL_MATH_HPP
#define SPATIAL_MATH_HPP

#include <limits>
#include <cmath>
#include <sstream>
#include "../exception.hpp"
#include "spatial_check_concept.hpp"

namespace spatial
{
  namespace except
  {
    /**
     *  Check that the distance given \c x has a positive value.
     *  \throws negative_distance
     */
    template<typename Tp>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Tp> >::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Tp> >::type
#endif
    check_positive_distance(Tp x)
    {
      if (x < Tp()) // Tp() == 0 by convention
        {
          std::stringstream out;
          out << x << " is negative";
          throw invalid_distance(out.str());
        }
    }

    /**
     *  This arithmetic check is only used when the macro
     *  SPATIAL_SAFER_ARITHMETICS is defined. Check that the absolute of an
     *  element has not led to an error such as an overflow, by forcing the
     *  error itself.
     *
     *  This check is not the best check for arithmetic errors. There are
     *  probably ways to make it faster, but it is intended to be portable and
     *  provide users with the possibility to quickly check the architmectics
     *  during computation with little efforts from their part.
     *
     *  The std::abs() function is working fine for floating point types,
     *  however, for signed integral types (char, wchar_t, short, int, long,
     *  long long) and their signed version, it returns an incorrect value when
     *  trying to compute std::abs(std::numeric_limits<signed int>::min()). To
     *  signal this issue we raise an exception in this case.
     */
    ///@{
    template <typename Tp>
    inline typename enable_if_c
    <std::numeric_limits<Tp>::is_integer
    && std::numeric_limits<Tp>::is_signed, Tp>::type
    check_abs(Tp x)
    {
      if (x == (std::numeric_limits<Tp>::min)())
        {
          std::stringstream out;
          out << "absolute of " << x << " caused overflow";
          throw arithmetic_error(out.str());
        }
      return std::abs(x);
    }
    template <typename Tp>
    inline typename enable_if_c
    <!std::numeric_limits<Tp>::is_integer
    || !std::numeric_limits<Tp>::is_signed, Tp>::type
    check_abs(Tp x)
    { return std::abs(x); }
    ///@}

    /**
     *  This arithmetic check is only used when the macro
     *  SPATIAL_SAFER_ARITHMETICS is defined. Check that the addtion of 2
     *  elements of positive value has not led to an error such as an overflow,
     *  by forcing the error itself.
     *
     *  This check is not the best check for arithmetic errors. There are ways
     *  to make it faster, but it is intended to be portable and provide users
     *  with the possibility to quickly check the architmectics during
     *  computation with little efforts from their part.
     *
     *  In particular, if \c Tp is not a base type, the author of the type must
     *  define the numeric limits \c numeric_limits<Tp>::max() for that type.
     */
    template <typename Tp>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Tp>, Tp>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Tp>, Tp>::type
#endif
    check_positive_add(Tp x, Tp y)
    {
      // The additional bracket is to avoid conflict with MSVC min/max macros
      if (((std::numeric_limits<Tp>::max)() - x) < y)
        {
          std::stringstream out;
          out << x << " + " << y << " caused overflow";
          throw arithmetic_error(out.str());
        }
      return x + y;
    }

    /**
     *  This arithmetic check is only used when the macro
     *  SPATIAL_SAFER_ARITHMETICS is defined. Check that the computation of the
     *  square of an element has not overflown.
     *
     *  This check is not the best check for arithmetic errors. There are ways
     *  to make it faster, but it is intended to be portable and provide users
     *  with the possibility to quickly check the architmectics during
     *  computation with little efforts from their part.
     */
    template <typename Tp>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Tp>, Tp>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Tp>, Tp>::type
#endif
    check_square(Tp x)
    {
      Tp zero = Tp(); // get 0
      if (x != zero)
        {
          Tp abs = check_abs(x);
          if (((std::numeric_limits<Tp>::max)() / abs) < abs)
            {
              std::stringstream out;
              out << "square(" << x << ") caused overflow";
              throw arithmetic_error(out.str());
            }
          return x * x;
        }
      return zero;
    }

    /**
     *  This arithmetic check is only used when the macro
     *  SPATIAL_SAFER_ARITHMETICS is defined. Check that the multiplication of 2
     *  positive elements has not resulted into an arithmetic error such as
     *  overflown.
     *
     *  This check will only work for 2 positive element x and y. This check is
     *  not the best check for arithmetic errors. There are ways to make it
     *  better, but it's hard to make it more portable.
     *
     *  In particular, if \c Tp is not a base type, the author of the type must
     *  define the numeric limits \c std::numeric_limits<Tp>::max() for that
     *  type.
     */
    template <typename Tp>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Tp>, Tp>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Tp>, Tp>::type
#endif
    check_positive_mul(Tp x, Tp y)
    {
      Tp zero = Tp(); // get 0
      if (x != zero)
        {
          if (((std::numeric_limits<Tp>::max)() / x) < y)
            {
              std::stringstream out;
              out << x << " * " << y << " caused overflow";
              throw arithmetic_error(out.str());
            }
          return x * y;
        }
      return zero;
    }
  } // spatial except

  namespace math
  {
    /**
     *  Compute the distance between the \p origin and the closest point to the
     *  plane orthogonal to the axis of dimension \c dim and passing by \c key.
     */
    template <typename Key, typename Difference, typename Unit>
#ifdef __APPLE__
    inline typename enable_if<std::is_floating_point<Unit>, Unit>::type
#else
    inline typename enable_if<std::tr1::is_floating_point<Unit>, Unit>::type
#endif
    euclid_distance_to_plane
    (dimension_type dim, Key origin, Key key, Difference diff)
    {
      return std::abs(diff(dim, origin, key)); // floating types abs is always okay!
    }

    /**
     *  Uses the hypot() algorithm in order to compute the distance: minimize
     *  possibilities of loss of precision due to overflow and underflow.
     *
     *  The trick is to find the maximum value among all the component of the
     *  Key and then divide all other component with this one.
     *
     *  Here is the rationale for the distance:
     *  sqrt( x^2 + y^2 + z^2 + ... ) = abs(x) * sqrt( 1 + (y/x)^2 + (z/x)^2 )
     *
     *  Providing x statisfies x>y, x>z, etc, the second form is less likely to
     *  overflow than the first form.
     */
    template <typename Key, typename Difference, typename Unit>
#ifdef __APPLE__
    inline typename enable_if<std::is_floating_point<Unit>, Unit>::type
#else
    inline typename enable_if<std::tr1::is_floating_point<Unit>, Unit>::type
#endif
    euclid_distance_to_key
    (dimension_type rank, Key origin, Key key, Difference diff)
    {
      // Find a non zero maximum or return 0
      Unit max = euclid_distance_to_plane<Key, Difference, Unit>
        (0, origin, key, diff);
      dimension_type max_dim = 0;
      for (dimension_type i = 1; i < rank; ++i)
        {
          Unit d = euclid_distance_to_plane<Key, Difference, Unit>
            (i, origin, key, diff);
          if (d > max) { max = d; max_dim = i; }
        }
      const Unit zero = Unit();
      if (max == zero) return zero; // they're all zero!
      // Compute the distance
      Unit sum = zero;
      for (dimension_type i = 0; i < rank; ++i)
        {
          if (i == max_dim) continue;
          Unit div = diff(i, origin, key) / max;
          sum += div * div;
        }
      const Unit one = ((Unit) 1.0); // Not sure how to represent it otherwise
#ifdef SPATIAL_SAFER_ARITHMETICS
      return except::check_positive_mul(max, std::sqrt(one + sum));
#else
      return max * std::sqrt(one + sum);
#endif
    }

    /**
     *  Compute the distance between the \p origin and the closest point
     *  to the plane orthogonal to the axis of dimension \c dim and passing by
     *  \c key.
     */
    template <typename Key, typename Difference, typename Unit>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Unit>, Unit>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Unit>, Unit>::type
#endif
    square_euclid_distance_to_plane
    (dimension_type dim, Key origin, Key key, Difference diff)
    {
      Unit d = diff(dim, origin, key);
#ifdef SPATIAL_SAFER_ARITHMETICS
      return except::check_square(d);
#else
      return d * d;
#endif
    }

    /**
     *  Compute the square value of the distance between \p origin and
     *  \p key.
     */
    template <typename Key, typename Difference, typename Unit>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Unit>, Unit>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Unit>, Unit>::type
#endif
    square_euclid_distance_to_key
    (dimension_type rank, const Key& origin, const Key& key,
     const Difference& diff)
    {
      Unit sum = square_euclid_distance_to_plane<Key, Difference, Unit>
        (0, origin, key, diff);
      for (dimension_type i = 1; i < rank; ++i)
        {
#ifdef SPATIAL_SAFER_ARITHMETICS
          sum = except::check_positive_add
            (square_euclid_distance_to_plane<Key, Difference, Unit>
             (i, origin, key, diff), sum);
#else
          sum += square_euclid_distance_to_plane
            <Key, Difference, Unit>(i, origin, key, diff);
#endif
        }
      return sum;
    }

    /**
     *  Compute the distance between the \p origin and the closest point
     *  to the plane orthogonal to the axis of dimension \c dim and passing by
     *  \c key.
     */
    template <typename Key, typename Difference, typename Unit>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Unit>, Unit>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Unit>, Unit>::type
#endif
    manhattan_distance_to_plane
    (dimension_type dim, Key origin, Key key, Difference diff)
    {
#ifdef SPATIAL_SAFER_ARITHMETICS
      return except::check_abs(diff(dim, origin, key));
#else
      return std::abs(diff(dim, origin, key));
#endif
    }

    /**
     *  Compute the manhattan distance between \p origin and \p key.
     */
    template <typename Key, typename Difference, typename Unit>
#ifdef __APPLE__
    inline typename enable_if<std::is_arithmetic<Unit>, Unit>::type
#else
    inline typename enable_if<std::tr1::is_arithmetic<Unit>, Unit>::type
#endif
    manhattan_distance_to_key
    (dimension_type rank, Key origin, Key key, Difference diff)
    {
      Unit sum = manhattan_distance_to_plane<Key, Difference, Unit>
        (0, origin, key, diff);
      for (dimension_type i = 1; i < rank; ++i)
        {
#ifdef SPATIAL_SAFER_ARITHMETICS
          sum = except::check_positive_add
            (manhattan_distance_to_plane<Key, Difference, Unit>
             (i, origin, key, diff), sum);
#else
          sum += manhattan_distance_to_plane
            <Key, Difference, Unit>(i, origin, key, diff);
#endif
        }
      return sum;
    }

    /*
      // For a future implementation where we take earth-like spheroid as an
      // example for non-euclidian spaces, or manifolds.
      math::great_circle_distance_to_key
      math::great_circle_distance_to_plane
      math::vincenty_distance_to_key
      math::vincenty_distance_to_plane
    */
  } // namespace math
}

#endif // SPATIAL_MATH_HPP