Commit f3ccc293 authored by Jason Rhinelander's avatar Jason Rhinelander

single_peak_search: return struct

Returning just the argmax is sort of inconvenient when the max is always
wanted: this returns a struct that has both .max and .arg elements.  It
should retain backwards compatibility in most cases by making the struct
implicitly convertible to a double, yielding the `.arg` value, the same
as before.
parent daf43734
......@@ -60,7 +60,7 @@ 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_;
double single_peak_search(
single_peak_result single_peak_search(
const std::function<double(const double &)> &f,
double left,
double right,
......@@ -70,6 +70,9 @@ double single_peak_search(
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);
......@@ -83,6 +86,7 @@ double single_peak_search(
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);
......@@ -92,6 +96,7 @@ double single_peak_search(
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);
......@@ -108,11 +113,31 @@ double single_peak_search(
}
}
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).
return fl >= fml && fl >= fmr && fl >= fr ? left :
fr >= fmr && fr >= fml ? right :
fml >= fmr ? midleft : midright;
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;
}
}
......@@ -248,6 +248,8 @@ class Stepper final {
};
struct single_peak_result; // forward declaration
/** 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
* work reliably if `f()` has multiple maximum in `[left, right]`; you'll just get a local maximum.
......@@ -279,12 +281,26 @@ class Stepper final {
* \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.
*
* \return a single_peak_result struct with `.arg` set to the peak argument, and `.max` set to the
* value of the function at `.arg`.
*/
double single_peak_search(
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`
};
}
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