Commit 7aa8b705 authored by Jason Rhinelander's avatar Jason Rhinelander

Added single_peak_search code in eris/algorithms

This uses a golden ratio section search to efficiently search a domain
for a function max.  Obviously it only works when the function is
single-peaked and has no perfectly-flat sections.

This also includes a test script.
parent a7c34cfe
......@@ -238,5 +238,35 @@ class Stepper final {
};
/** 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.
* It may also fail if the function has perfectly flat sections in the given domain.
*
* This algorithm works by considering two interior point at proportions \f$2 - \varphi\f$ and
* \f$\varphi - 1\f$ between `left` and `right`, where \f$\varphi = \frac{1+\sqrt{5}}{2}\f$, the
* golden ratio. Whichever interior point yields a larger function value becomes one of the new
* interior points, while the other becomes the new `left` or `right` value as appropriate. This
* process repeats until the tolerance level is reached (see the `tolerance` argument, below), at
* which point whichever of `left`, `right`, and the two midpoints has the greatest `f()` value is
* returned.
*
* Each iteration of this algorithm removes \f$2-\varphi \approx 0.382\f$ of the range and requires
* only a single additional function evaluation due to the use of the golden ratio (because one of
* 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 left the left edge of the domain to consider
* \param right the right edge of the domain to consider
* \param tolerance the relative size of the domain at which the algorithm stops. In particular,
* the algorithm stops once \f$\frac{right - left}{max\{\|left\|, \|right\|\}} \leq tolerance\f$. The default, if the
* argument is omitted, is \f$10^{-12}\f$.
*/
extern double single_peak_search(
const std::function<double(const double &)> &f,
const double &left, const double &right,
const double &tolerance = 1e-12);
}
#include <eris/algorithms.hpp>
#include <cmath>
namespace eris {
Stepper::Stepper(double step, int increase_count, double min_step, double max_step, bool rel_steps)
......@@ -51,4 +53,52 @@ double Stepper::step(bool up) {
else return 1.0 / (1 + step_size);
}
// Relative left midpoint position. Equal to \f$2 - \varphi\f$.
const double sps_left_mid_ = 1.5 - std::sqrt(1.25);
// Relative right midpoint position. Equal to \f$\varphi - 1\f$.
const double sps_right_mid_ = 1 - sps_left_mid_;
double sps_(
const std::function<double(const double &)> &f,
const double &left, const double &midleft, const double &midright, const double &right,
const double &fl, const double &fml, const double &fmr, const double &fr,
const double &tolerance) {
if ((right - left) / std::max(std::fabs(left), std::fabs(right)) <= tolerance) {
// 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 (fml >= fmr) {
// midleft is the higher point, so we can exclude everything right of midright.
// midright -> new right
// midleft -> new midright
// and calculate the new midleft position and evaluate it for the next iteration
const double newmidleft = left + (midright - left) * sps_left_mid_;
return sps_(f, left, newmidleft, midleft, midright, fl, f(newmidleft), fml, fmr, tolerance);
}
else {
// midright is higher, so exclude everything left of midleft.
// midleft -> new left
// midright -> new midleft
// and calculated the new midright position and evaluate for the next iteration
const double newmidright = midleft + (right - midleft) * sps_right_mid_;
return sps_(f, midleft, midright, newmidright, right, fml, fmr, f(newmidright), fr, tolerance);
}
}
double single_peak_search(
const std::function<double(const double &)> &f,
const double &left, const double &right,
const double &tolerance) {
const double size = right - left;
const double midleft = left + sps_left_mid_*size;
const double midright = left + sps_right_mid_*size;
return sps_(f, left, midleft, midright, right, f(left), f(midleft), f(midright), f(right), tolerance);
}
}
......@@ -12,6 +12,7 @@ set(eris_tests
sim-setup-test
lock-test
pos-agent-test
single-peak-search
)
foreach(t ${eris_tests})
......
#include <eris/algorithms.hpp>
#include <gtest/gtest.h>
using eris::single_peak_search;
TEST(Maximize, Quadratic) {
auto f = [] (const double &x) -> double { return -3*x*x + 14*x - 3; };
EXPECT_NEAR(14.0/6.0, single_peak_search(f, -10, 10), 2e-8);
EXPECT_NEAR(14.0/6.0, single_peak_search(f, -100, 3.0), 2e-8);
EXPECT_NEAR(14.0/6.0, single_peak_search(f, 0, 100000), 2e-8);
EXPECT_NEAR(14.0/6.0, single_peak_search(f, 2.3, 2.4), 3e-9);
}
TEST(Maximize, Quartic) {
auto f = [] (const double &x) -> double { const double x2 = x*x; return -21237*x2*x2 + 13*x2 - 1247*x + 3; };
EXPECT_NEAR(-0.24526910870656568, single_peak_search(f, -10, 10), 2e-9);
EXPECT_NEAR(-0.24526910870656568, single_peak_search(f, -1, 0), 4e-9);
EXPECT_NEAR(-0.24526910870656568, single_peak_search(f, -1e10, 1e10), 7e-10);
EXPECT_NEAR(-0.24526910870656568, single_peak_search(f, -0.246, 0.245), 2e-9);
}
TEST(Maximize, LeftEndPoint) {
auto f = [] (const double &x) -> double { return 100 - 12*x; };
// NB: not using DOUBLE_EQ here, these should be exact matches
EXPECT_EQ(-14.675, single_peak_search(f, -14.675, 10000));
EXPECT_EQ(-12, single_peak_search(f, -12, -1));
EXPECT_EQ(2000, single_peak_search(f, 2000, 50000));
}
TEST(Maximize, CubicRightEndPoint) {
auto f = [] (const double &x) -> double { const double x2 = x*x; return x2*x - 2*x2 + 3*x + 17; };
// NB: not using DOUBLE_EQ here, these should be exact matches
EXPECT_EQ(10000, single_peak_search(f, -14.675, 10000));
EXPECT_EQ(-1, single_peak_search(f, -12, -1));
EXPECT_EQ(1e100, single_peak_search(f, -1e100, 1e100));
}
TEST(Maximize, PosQuadraticEnd) {
auto f = [] (const double &x) -> double { return x*x + 14*x + 70; };
// We should always end up at one end-point or the other, but which side depends on the initial
// range. NB: the above has a minimum at -7, but we should never get that back.
EXPECT_EQ(10, single_peak_search(f, -10, 10));
EXPECT_EQ(-11, single_peak_search(f, -11, -6));
EXPECT_EQ(-2.875, single_peak_search(f, -11, -2.875));
}
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