Commit 2bf6d8d3 authored by Jason Rhinelander's avatar Jason Rhinelander

Make single_peak_search templated; remove abs_tol

This makes single_peak_search templated (so that it doesn't have to be
implemented with double).  It also updates the return type struct with
templated types, a constructor, and an iterator counter.
parent 0df6e3e8
#include <eris/algorithms.hpp>
#include <cmath>
#include <algorithm>
#include <boost/math/constants/constants.hpp>
namespace eris {
......@@ -55,89 +52,4 @@ double Stepper::step(bool up) {
else return 1.0 / (1 + step_size);
}
// Relative right midpoint position. Equal to \f$\varphi - 1\f$.
constexpr double sps_right_mid_ = boost::math::constants::phi<double>() - 1;
// Relative left midpoint position. Equal to 1 minus `sps_right_mid_`, which also equals 2 - phi.
constexpr double sps_left_mid_ = 1 - sps_right_mid_;
single_peak_result single_peak_search(
const std::function<double(const double &)> &f,
double left,
double right,
double tol_rel,
double tol_abs) {
if (tol_rel < 0) tol_rel = 0;
if (tol_abs < 0) tol_abs = 0;
// Track whether the peak is strictly inside the initial boundaries:
bool inside_left = false, inside_right = false;
double midleft = left + sps_left_mid_*(right - left);
double midright = left + sps_right_mid_*(right - left);
double fl = f(left), fml = f(midleft), fmr = f(midright), fr = f(right);
while (
tol_abs < right - left
and
tol_rel < (right - left) / std::max(std::fabs(left), std::fabs(right))
) {
if (fml >= fmr) {
// midleft is the higher point, so we can exclude everything right of midright.
right = midright; fr = fmr;
inside_right = true;
midright = midleft; fmr = fml;
midleft = left + (right - left) * sps_left_mid_;
fml = f(midleft);
// If the midpoint is closer to the endpoint than is numerically distinguishable, finish:
if (midleft == left) break;
}
else {
// midright is higher, so exclude everything left of midleft.
left = midleft; fl = fml;
inside_left = true;
midleft = midright; fml = fmr;
midright = left + (right - left) * sps_right_mid_;
fmr = f(midright);
// If the midpoint is closer to the endpoint than is numerically distinguishable, finish:
if (midright == right) break;
}
// Sometimes, we can run into numerical instability that results in midleft > midright,
// particular when the optimum is close to 0. If we encounter that, swap midleft and
// midright.
if (midleft > midright) {
std::swap(midleft, midright);
std::swap(fml, fmr);
}
}
single_peak_result result;
// Prefer the end-points for ties (the max might legitimate be an end-point), and prefer
// left over right (for no particularly good reason).
if (fl >= fml && fl >= fmr && fl >= fr) {
result.arg = left;
result.max = fl;
result.inside = inside_left;
}
else if (fr >= fmr && fr >= fml) {
result.arg = right;
result.max = fr;
result.inside = inside_right;
}
else if (fml >= fmr) {
result.arg = midleft;
result.max = fml;
result.inside = true;
}
else {
result.arg = midright;
result.max = fmr;
result.inside = true;
}
return result;
}
}
......@@ -6,6 +6,7 @@
#include <vector>
#include <limits>
#include <utility>
#include <cmath>
namespace eris {
......@@ -247,7 +248,30 @@ class Stepper final {
};
struct single_peak_result; // forward declaration
/// Struct holding the results of a call to an optimization function such as single_peak_search() or
/// `constrained_maximum_search()`
template <typename ArgType = double, typename ValueType = double> struct search_result {
ArgType arg; ///< The argument that maximizes the searched function
ValueType value; ///< The value of the function at `.arg`
/** Whether `.arg` is strictly inside the given left/right limits. If false, the found value
* was at one of the end-points, and may not actually be a peak at all.
*/
bool inside;
/// Number of iterations, if applicable (-1 otherwise).
int iterations;
operator ArgType() const { return arg; } ///< Implicit conversion to double returns `.arg`
search_result(ArgType a, ValueType v, bool in, int it = -1) : arg(std::move(a)), value(std::move(v)), inside(in), iterations(it) {}
};
/// The constant phi. Callers can specialize this template if using custom types with more
/// precision than a long double value.
template <typename RealType> constexpr RealType phi = 1.61803398874989484820458683436563811L;
/// The right inner point multiple for a golden section search, phi - 1.
template <typename RealType> constexpr RealType golden_section_right = phi<RealType> - RealType(1);
/// The left inner point multiple for a golden section search, 1 - right, which also equals 2 - phi.
template <typename RealType> constexpr RealType golden_section_left = RealType(1) - golden_section_right<RealType>;
/** Performs a golden section search to find a maximum of a single-peaked function between two
* limits. This function must be called with left and right end points. This function will not
......@@ -267,39 +291,85 @@ struct single_peak_result; // forward declaration
* the prior midpoints will also be the midpoint of the next iteration, and so the function value
* can simply be reused).
*
* \param f a function-like object that can be called with a single double argument and returns the
* function value.
* \param f a function or function-like object that can be called with a single argument value and
* returns the function value at that argument.
* \param left the left edge of the domain to consider
* \param right the right edge of the domain to consider
* \param tol_rel the relative size of the domain at which the algorithm stops. In particular, the
* algorithm stops if \f$\frac{right - left}{max\{\|left\|, \|right\|\}} \leq tol_{rel}\f$. The
* default, if the argument is omitted, is \f$10^{-10}\f$. Note that the algorithm might
* alternatively stop because of the `tol_abs` value. Note also that there is also a numerical
* lower bound on this value: if the midpoint calculation results in a midpoint exactly numerically
* equal to an end point, the algorithm stops.
* \param tol_abs the absolute size of the domain at which the algorithm stops. In particular, the
* algorithm stops if \f$right - left \leq tol_{abs}\f$. Note that the algorithm might
* alternatively stop because of the `tol_rel` value.
* default, if the argument is omitted, is \f$10^{-10}\f$. Note that the algorithm will also stop
* when it reaches the limits of numerical precision, that is, where a calculated midpoint does not
* numerically differ from an end point (and so you can safely specify a tolerance of 0 to get
* maximum double precision).
*
* \return a single_peak_result struct with `.arg` set to the peak argument, and `.max` set to the
* value of the function at `.arg`.
* \return a `eris::search_result` struct with `.arg` set to the peak argument, and `.max` set to
* the value of the function at `.arg`.
*/
single_peak_result single_peak_search(
const std::function<double(const double &)> &f,
double left,
double right,
double tol_rel = 1e-10,
double tol_abs = 1e-20);
/// Struct holding the results of a call to single_peak_search()
struct single_peak_result {
double arg; ///< The argument that maximizes the function given to `single_peak_search`
double max; ///< The value of the function at `.arg`
/** Whether `.arg` is strictly inside the given left/right limits. If false, the peak was at
* one of the end-points, and may not actually be a peak at all.
*/
bool inside;
operator double() const { return arg; } ///< Implicit conversion to double returns `.arg`
};
template <typename ArgT, typename Func, typename ValueT = decltype(std::declval<Func>()(std::declval<ArgT>()))>
search_result<ArgT, ValueT> single_peak_search(
Func f, ArgT left, ArgT right, std::common_type_t<ArgT> tol_rel = 1e-10) {
if (tol_rel < 0) tol_rel = 0;
constexpr ArgT midpoint_right = phi<ArgT> - ArgT(1);
constexpr ArgT midpoint_left = ArgT(1) - midpoint_right;
// Track whether the peak is strictly inside the initial boundaries:
bool inside_left = false, inside_right = false;
ArgT span = right - left;
ArgT midleft = left + midpoint_left * span;
ArgT midright = left + midpoint_right * span;
ValueT fl = f(left), fml = f(midleft), fmr = f(midright), fr = f(right);
int iterations = 1; // Count the above mid calcs as an iteration
using std::abs; // Don't use std::abs directly (to allow ADL on abs)
using std::max;
while (span > tol_rel * max(abs(left), abs(right))) {
iterations++;
if (fml >= fmr) {
// midleft is the higher point, so we can exclude everything right of midright.
right = std::move(midright); fr = std::move(fmr);
inside_right = true;
span = right - left;
midright = std::move(midleft); fmr = std::move(fml);
midleft = left + midpoint_left * span;
fml = f(midleft);
// If the midpoint is closer to the endpoint than is numerically distinguishable, finish:
if (midleft == left) break;
}
else {
// midright is higher, so exclude everything left of midleft.
left = std::move(midleft); fl = std::move(fml);
inside_left = true;
span = right - left;
midleft = std::move(midright); fml = std::move(fmr);
midright = left + midpoint_right * span;
fmr = f(midright);
// If the midpoint is closer to the endpoint than is numerically distinguishable, finish:
if (midright == right) break;
}
// Sometimes we can run into numerical instability that results in midleft > midright,
// particular when the optimum is close to 0. If we encounter that, swap midleft and
// midright.
if (midleft > midright) {
using std::swap;
swap(midleft, midright);
swap(fml, fmr);
}
}
// Prefer the end-points for ties (the max might legitimately be an end-point), and prefer left
// over right (for no particularly good reason).
if (fl >= fml && fl >= fmr && fl >= fr)
return {std::move(left), std::move(fl), inside_left, iterations};
else if (fr >= fmr && fr >= fml)
return {std::move(right), std::move(fr), inside_right, iterations};
else if (fml >= fmr)
return {std::move(midleft), std::move(fml), true, iterations};
else
return {std::move(midright), std::move(fmr), true, iterations};
}
}
#include <eris/algorithms.hpp>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <sstream>
#ifdef SEARCH_F128
extern "C" {
#include <quadmath.h>
}
#endif
template <typename T> std::string tostr(const T &v) {
std::ostringstream s;
s << std::setprecision(20) << std::boolalpha << v;
return s.str();
}
#ifdef SEARCH_F128
template <> std::string tostr(const __float128 &v) {
char buf[128];
quadmath_snprintf(buf, sizeof buf, "%.40Qf", v);
return std::string(buf);
}
#endif
template <typename R> void print(const R &r) {
std::cout << "arg=" << tostr(r.arg) << ", val=" << tostr(r.value) << ", inside=" << r.inside << ", iterations=" << r.iterations << "\n";
}
#ifdef SEARCH_F128
template <> constexpr __float128 eris::phi<__float128> =
1.61803398874989484820458683436563811Q;
#endif
int main() {
#ifndef SEARCH_F128
using FloatType = double;
#else
using FloatType = __float128;
#endif
print(eris::single_peak_search<FloatType>(
[](FloatType x) { return 12 - x*x + 3*x; },
0, 100, 0));
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment