//  (C) Copyright John Maddock 2006.
//  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_TOOLS_REMEZ_HPP
#define BOOST_MATH_TOOLS_REMEZ_HPP

#include <boost/math/tools/solve.hpp>
#include <boost/math/tools/minima.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/tools/polynomial.hpp>
#include <boost/function/function1.hpp>
#include <boost/scoped_array.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/policy.hpp>

namespace boost{ namespace math{ namespace tools{

namespace detail{

//
// The error function: the difference between F(x) and
// the current approximation.  This is the function
// for which we must find the extema.
//
template <class T>
struct remez_error_function
{
   typedef boost::function1<T, T const &> function_type;
public:
   remez_error_function(
      function_type f_, 
      const polynomial<T>& n, 
      const polynomial<T>& d, 
      bool rel_err)
         : f(f_), numerator(n), denominator(d), rel_error(rel_err) {}

   T operator()(const T& z)const
   {
      T y = f(z);
      T abs = y - (numerator.evaluate(z) / denominator.evaluate(z));
      T err;
      if(rel_error)
      {
         if(y != 0)
            err = abs / fabs(y);
         else if(0 == abs)
         {
            // we must be at a root, or it's not recoverable:
            BOOST_ASSERT(0 == abs);
            err = 0;
         }
         else
         {
            // We have a divide by zero!
            // Lets assume that f(x) is zero as a result of
            // internal cancellation error that occurs as a result
            // of shifting a root at point z to the origin so that
            // the approximation can be "pinned" to pass through
            // the origin: in that case it really
            // won't matter what our approximation calculates here
            // as long as it's a small number, return the absolute error:
            err = abs;
         }
      }
      else
         err = abs;
      return err;
   }
private:
   function_type f;
   polynomial<T> numerator;
   polynomial<T> denominator;
   bool rel_error;
};
//
// This function adapts the error function so that it's minima
// are the extema of the error function.  We can find the minima
// with standard techniques.
//
template <class T>
struct remez_max_error_function
{
   remez_max_error_function(const remez_error_function<T>& f)
      : func(f) {}

   T operator()(const T& x)
   {
      BOOST_MATH_STD_USING
      return -fabs(func(x));
   }
private:
   remez_error_function<T> func;
};

} // detail

template <class T>
class remez_minimax
{
public:
   typedef boost::function1<T, T const &> function_type;
   typedef boost::numeric::ublas::vector<T> vector_type;
   typedef boost::numeric::ublas::matrix<T> matrix_type;

   remez_minimax(function_type f, unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0);
   remez_minimax(function_type f, unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points);

   void reset(unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0);
   void reset(unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points);

   void set_brake(int b)
   {
      BOOST_ASSERT(b < 100);
      BOOST_ASSERT(b >= 0);
      m_brake = b;
   }

   T iterate();

   polynomial<T> denominator()const;
   polynomial<T> numerator()const;

   vector_type const& chebyshev_points()const
   {
      return control_points;
   }

   vector_type const& zero_points()const
   {
      return zeros;
   }

   T error_term()const
   {
      return solution[solution.size() - 1];
   }
   T max_error()const
   {
      return m_max_error;
   }
   T max_change()const
   {
      return m_max_change;
   }
   void rotate()
   {
      --orderN;
      ++orderD;
   }
   void rescale(T a, T b)
   {
      T scale = (b - a) / (max - min);
      for(unsigned i = 0; i < control_points.size(); ++i)
      {
         control_points[i] = (control_points[i] - min) * scale + a;
      }
      min = a;
      max = b;
   }
private:

   void init_chebyshev();

   function_type func;            // The function to approximate.
   vector_type control_points;    // Current control points to be used for the next iteration.
   vector_type solution;          // Solution from the last iteration contains all unknowns including the error term.
   vector_type zeros;             // Location of points of zero error from last iteration, plus the two end points.
   vector_type maxima;            // Location of maxima of the error function, actually contains the control points used for the last iteration.
   T m_max_error;                 // Maximum error found in last approximation.
   T m_max_change;                // Maximum change in location of control points after last iteration.
   unsigned orderN;               // Order of the numerator polynomial.
   unsigned orderD;               // Order of the denominator polynomial.
   T min, max;                    // End points of the range to optimise over.
   bool rel_error;                // If true optimise for relative not absolute error.
   bool pinned;                   // If true the approximation is "pinned" to go through the origin.
   unsigned unknowns;             // Total number of unknowns.
   int m_precision;               // Number of bits precision to which the zeros and maxima are found.
   T m_max_change_history[2];     // Past history of changes to control points.
   int m_brake;                     // amount to break by in percentage points.
   int m_skew;                      // amount to skew starting points by in percentage points: -100-100
};

#ifndef BRAKE
#define BRAKE 0
#endif
#ifndef SKEW
#define SKEW 0
#endif

template <class T>
void remez_minimax<T>::init_chebyshev()
{
   BOOST_MATH_STD_USING
   //
   // Fill in the zeros:
   //
   unsigned terms = pinned ? orderD + orderN : orderD + orderN + 1;

   for(unsigned i = 0; i < terms; ++i)
   {
      T cheb = cos((2 * terms - 1 - 2 * i) * constants::pi<T>() / (2 * terms));
      cheb += 1;
      cheb /= 2;
      if(m_skew != 0)
      {
         T p = static_cast<T>(200 + m_skew) / 200;
         cheb = pow(cheb, p);
      }
      cheb *= (max - min);
      cheb += min;
      zeros[i+1] = cheb;
   }
   zeros[0] = min;
   zeros[unknowns] = max;
   // perform a regular interpolation fit:
   matrix_type A(terms, terms);
   vector_type b(terms);
   // fill in the y values:
   for(unsigned i = 0; i < b.size(); ++i)
   {
      b[i] = func(zeros[i+1]);
   }
   // fill in powers of x evaluated at each of the control points:
   unsigned offsetN = pinned ? 0 : 1;
   unsigned offsetD = offsetN + orderN;
   unsigned maxorder = (std::max)(orderN, orderD);
   for(unsigned i = 0; i < b.size(); ++i)
   {
      T x0 = zeros[i+1];
      T x = x0;
      if(!pinned)
         A(i, 0) = 1;
      for(unsigned j = 0; j < maxorder; ++j)
      {
         if(j < orderN)
            A(i, j + offsetN) = x;
         if(j < orderD)
         {
            A(i, j + offsetD) = -x * b[i];
         }
         x *= x0;
      }
   }
   //
   // Now go ahead and solve the expression to get our solution:
   //
   vector_type l_solution = boost::math::tools::solve(A, b);
   // need to add a "fake" error term:
   l_solution.resize(unknowns);
   l_solution[unknowns-1] = 0;
   solution = l_solution;
   //
   // Now find all the extrema of the error function:
   //
   detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error);
   detail::remez_max_error_function<T> Ex(Err);
   m_max_error = 0;
   int max_err_location = 0;
   for(unsigned i = 0; i < unknowns; ++i)
   {
      std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision);
      maxima[i] = r.first;
      T rel_err = fabs(r.second);
      if(rel_err > m_max_error)
      {
         m_max_error = fabs(r.second);
         max_err_location = i;
      }
   }
   control_points = maxima;
}

template <class T>
void remez_minimax<T>::reset(
         unsigned oN, 
         unsigned oD, 
         T a, 
         T b, 
         bool pin, 
         bool rel_err, 
         int sk,
         int bits)
{
   control_points = vector_type(oN + oD + (pin ? 1 : 2));
   solution = control_points;
   zeros = vector_type(oN + oD + (pin ? 2 : 3));
   maxima = control_points;
   orderN = oN;
   orderD = oD;
   rel_error = rel_err;
   pinned = pin;
   m_skew = sk;
   min = a;
   max = b;
   m_max_error = 0;
   unknowns = orderN + orderD + (pinned ? 1 : 2);
   // guess our initial control points:
   control_points[0] = min;
   control_points[unknowns - 1] = max;
   T interval = (max - min) / (unknowns - 1);
   T spot = min + interval;
   for(unsigned i = 1; i < control_points.size(); ++i)
   {
      control_points[i] = spot;
      spot += interval;
   }
   solution[unknowns - 1] = 0;
   m_max_error = 0;
   if(bits == 0)
   {
      // don't bother about more than float precision:
      m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
   }
   else
   {
      // can't be more accurate than half the bits of T:
      m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
   }
   m_max_change_history[0] = m_max_change_history[1] = 1;
   init_chebyshev();
   // do one iteration whatever:
   //iterate();
}

template <class T>
inline remez_minimax<T>::remez_minimax(
         typename remez_minimax<T>::function_type f, 
         unsigned oN, 
         unsigned oD, 
         T a, 
         T b, 
         bool pin, 
         bool rel_err, 
         int sk,
         int bits)
   : func(f) 
{
   m_brake = 0;
   reset(oN, oD, a, b, pin, rel_err, sk, bits);
}

template <class T>
void remez_minimax<T>::reset(
         unsigned oN, 
         unsigned oD, 
         T a, 
         T b, 
         bool pin, 
         bool rel_err, 
         int sk,
         int bits,
         const vector_type& points)
{
   control_points = vector_type(oN + oD + (pin ? 1 : 2));
   solution = control_points;
   zeros = vector_type(oN + oD + (pin ? 2 : 3));
   maxima = control_points;
   orderN = oN;
   orderD = oD;
   rel_error = rel_err;
   pinned = pin;
   m_skew = sk;
   min = a;
   max = b;
   m_max_error = 0;
   unknowns = orderN + orderD + (pinned ? 1 : 2);
   control_points = points;
   solution[unknowns - 1] = 0;
   m_max_error = 0;
   if(bits == 0)
   {
      // don't bother about more than float precision:
      m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
   }
   else
   {
      // can't be more accurate than half the bits of T:
      m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
   }
   m_max_change_history[0] = m_max_change_history[1] = 1;
   // do one iteration whatever:
   //iterate();
}

template <class T>
inline remez_minimax<T>::remez_minimax(
         typename remez_minimax<T>::function_type f, 
         unsigned oN, 
         unsigned oD, 
         T a, 
         T b, 
         bool pin, 
         bool rel_err, 
         int sk,
         int bits,
         const vector_type& points)
   : func(f)
{
   m_brake = 0;
   reset(oN, oD, a, b, pin, rel_err, sk, bits, points);
}

template <class T>
T remez_minimax<T>::iterate()
{
   BOOST_MATH_STD_USING
   matrix_type A(unknowns, unknowns);
   vector_type b(unknowns);

   // fill in evaluation of f(x) at each of the control points:
   for(unsigned i = 0; i < b.size(); ++i)
   {
      // take care that none of our control points are at the origin:
      if(pinned && (control_points[i] == 0))
      {
         if(i)
            control_points[i] = control_points[i-1] / 3;
         else
            control_points[i] = control_points[i+1] / 3;
      }
      b[i] = func(control_points[i]);
   }

   T err_err;
   unsigned convergence_count = 0;
   do{
      // fill in powers of x evaluated at each of the control points:
      int sign = 1;
      unsigned offsetN = pinned ? 0 : 1;
      unsigned offsetD = offsetN + orderN;
      unsigned maxorder = (std::max)(orderN, orderD);
      T Elast = solution[unknowns - 1];

      for(unsigned i = 0; i < b.size(); ++i)
      {
         T x0 = control_points[i];
         T x = x0;
         if(!pinned)
            A(i, 0) = 1;
         for(unsigned j = 0; j < maxorder; ++j)
         {
            if(j < orderN)
               A(i, j + offsetN) = x;
            if(j < orderD)
            {
               T mult = rel_error ? (b[i] - sign * fabs(b[i]) * Elast): (b[i] - sign * Elast);
               A(i, j + offsetD) = -x * mult;
            }
            x *= x0;
         }
         // The last variable to be solved for is the error term, 
         // sign changes with each control point:
         T E = rel_error ? sign * fabs(b[i]) : sign;
         A(i, unknowns - 1) = E;
         sign = -sign;
      }

   #ifdef BOOST_MATH_INSTRUMENT
      for(unsigned i = 0; i < b.size(); ++i)
         std::cout << b[i] << " ";
      std::cout << "\n\n";
      for(unsigned i = 0; i < b.size(); ++i)
      {
         for(unsigned j = 0; j < b.size(); ++ j)
            std::cout << A(i, j) << " ";
         std::cout << "\n";
      }
      std::cout << std::endl;
   #endif
      //
      // Now go ahead and solve the expression to get our solution:
      //
      solution = boost::math::tools::solve(A, b);

      err_err = (Elast != 0) ? fabs((fabs(solution[unknowns-1]) - fabs(Elast)) / fabs(Elast)) : 1;
   }while(orderD && (convergence_count++ < 80) && (err_err > 0.001));

   //
   // Perform a sanity check to verify that the solution to the equations
   // is not so much in error as to be useless.  The matrix inversion can
   // be very close to singular, so this can be a real problem.
   //
   vector_type sanity = prod(A, solution);
   for(unsigned i = 0; i < b.size(); ++i)
   {
      T err = fabs((b[i] - sanity[i]) / fabs(b[i]));
      if(err > sqrt(epsilon<T>()))
      {
         std::cerr << "Sanity check failed: more than half the digits in the found solution are in error." << std::endl;
      }
   }

   //
   // Next comes another sanity check, we want to verify that all the control
   // points do actually alternate in sign, in practice we may have 
   // additional roots in the error function that cause this to fail.
   // Failure here is always fatal: even though this code attempts to correct
   // the problem it usually only postpones the inevitable.
   //
   polynomial<T> num, denom;
   num = this->numerator();
   denom = this->denominator();
   T e1 = b[0] - num.evaluate(control_points[0]) / denom.evaluate(control_points[0]);
#ifdef BOOST_MATH_INSTRUMENT
   std::cout << e1;
#endif
   for(unsigned i = 1; i < b.size(); ++i)
   {
      T e2 = b[i] - num.evaluate(control_points[i]) / denom.evaluate(control_points[i]);
#ifdef BOOST_MATH_INSTRUMENT
      std::cout << " " << e2;
#endif
      if(e2 * e1 > 0)
      {
         std::cerr << std::flush << "Basic sanity check failed: Error term does not alternate in sign, non-recoverable error may follow..." << std::endl;
         T perturbation = 0.05;
         do{
            T point = control_points[i] * (1 - perturbation) + control_points[i-1] * perturbation;
            e2 = func(point) - num.evaluate(point) / denom.evaluate(point);
            if(e2 * e1 < 0)
            {
               control_points[i] = point;
               break;
            }
            perturbation += 0.05;
         }while(perturbation < 0.8);

         if((e2 * e1 > 0) && (i + 1 < b.size()))
         {
            perturbation = 0.05;
            do{
               T point = control_points[i] * (1 - perturbation) + control_points[i+1] * perturbation;
               e2 = func(point) - num.evaluate(point) / denom.evaluate(point);
               if(e2 * e1 < 0)
               {
                  control_points[i] = point;
                  break;
               }
               perturbation += 0.05;
            }while(perturbation < 0.8);
         }

      }
      e1 = e2;
   }

#ifdef BOOST_MATH_INSTRUMENT
   for(unsigned i = 0; i < solution.size(); ++i)
      std::cout << solution[i] << " ";
   std::cout << std::endl << this->numerator() << std::endl;
   std::cout << this->denominator() << std::endl;
   std::cout << std::endl;
#endif

   //
   // The next step is to find all the intervals in which our maxima
   // lie:
   //
   detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error);
   zeros[0] = min;
   zeros[unknowns] = max;
   for(unsigned i = 1; i < control_points.size(); ++i)
   {
      eps_tolerance<T> tol(m_precision);
      boost::uintmax_t max_iter = 1000;
      std::pair<T, T> p = toms748_solve(
         Err, 
         control_points[i-1], 
         control_points[i], 
         tol, 
         max_iter);
      zeros[i] = (p.first + p.second) / 2;
      //zeros[i] = bisect(Err, control_points[i-1], control_points[i], m_precision);
   }
   //
   // Now find all the extrema of the error function:
   //
   detail::remez_max_error_function<T> Ex(Err);
   m_max_error = 0;
   int max_err_location = 0;
   for(unsigned i = 0; i < unknowns; ++i)
   {
      std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision);
      maxima[i] = r.first;
      T rel_err = fabs(r.second);
      if(rel_err > m_max_error)
      {
         m_max_error = fabs(r.second);
         max_err_location = i;
      }
   }
   //
   // Almost done now! we just need to set our control points
   // to the extrema, and calculate how much each point has changed
   // (this will be our termination condition):
   //
   swap(control_points, maxima);
   m_max_change = 0;
   int max_change_location = 0;
   for(unsigned i = 0; i < unknowns; ++i)
   {
      control_points[i] = (control_points[i] * (100 - m_brake) + maxima[i] * m_brake) / 100;
      T change = fabs((control_points[i] - maxima[i]) / control_points[i]);
#if 0
      if(change > m_max_change_history[1])
      {
         // divergence!!! try capping the change:
         std::cerr << "Possible divergent step, change will be capped!!" << std::endl;
         change = m_max_change_history[1];
         if(control_points[i] < maxima[i])
            control_points[i] = maxima[i] - change * maxima[i];
         else
            control_points[i] = maxima[i] + change * maxima[i];
      }
#endif
      if(change > m_max_change)
      {
         m_max_change = change;
         max_change_location = i;
      }
   }
   //
   // store max change information:
   //
   m_max_change_history[0] = m_max_change_history[1];
   m_max_change_history[1] = fabs(m_max_change);

   return m_max_change;
}

template <class T>
polynomial<T> remez_minimax<T>::numerator()const
{
   boost::scoped_array<T> a(new T[orderN + 1]);
   if(pinned)
      a[0] = 0;
   unsigned terms = pinned ? orderN : orderN + 1;
   for(unsigned i = 0; i < terms; ++i)
      a[pinned ? i+1 : i] = solution[i];
   return boost::math::tools::polynomial<T>(&a[0], orderN);
}

template <class T>
polynomial<T> remez_minimax<T>::denominator()const
{
   unsigned terms = orderD + 1;
   unsigned offsetD = pinned ? orderN : (orderN + 1);
   boost::scoped_array<T> a(new T[terms]);
   a[0] = 1;
   for(unsigned i = 0; i < orderD; ++i)
      a[i+1] = solution[i + offsetD];
   return boost::math::tools::polynomial<T>(&a[0], orderD);
}


}}} // namespaces

#endif // BOOST_MATH_TOOLS_REMEZ_HPP