// (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 | |
#ifdef _MSC_VER | |
#pragma once | |
#endif | |
#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 | |