algorithms.hpp 29.8 KB
 Jason Rhinelander committed May 31, 2013 1 2 3 4 5 6 #pragma once #include #include #include #include #include  Jason Rhinelander committed Jul 19, 2013 7 #include  Jason Rhinelander committed Aug 26, 2015 8 #include  Jason Rhinelander committed Aug 06, 2017 9 #include  Jason Rhinelander committed May 31, 2013 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39  namespace eris { /** Calculates all combinations of all sizes of the given container elements. For example, given * input {1,2,3} this yields combinations: {}, {1}, {1,2}, {1,2,3}, {1,3}, {2}, {2,3}, {3}. The * specific order of sets is not guaranteed (the order above matches the current implementation, but * could change); within each set, the order will be the same as occurs in the passed-in iterators. * * The specific algorithm used here (which is subject to change) works by keeping a stack of * iterators, starting with just one element at the first iterator value. In each loop iteration: * - the current set of stack elements defines an as-yet unseen combination, so call the function. * - Then we consider the top stack element and: * - if it is at the end, pop it off and increment the (new) end iterator (the stack will always * be in increasing order, so that will never go push the new top element past the end) * - otherwise, add a new iterator to the stack starting at the next position. * * This algorithm is efficient, requiring exactly 1 iteration for each possible non-empty * combination (of which there are 2^n-1 for a given input of size n). * * \param begin a forward iterator pointing to the first element of the container whose elements are to * be recombined. Typically this is begin() or cbegin(), but starting elsewhere (to * recombine only a subset of the container) is permitted. * \param end an iterator pointing at the past-the-end element. Typically this is end() or * cend(), but can be something else for recombining subsets. * \param func a void function that takes a single argument of a const std::vector &, where T is * the same type used by the start and end iterators. This will be called for each (unique) * combination of parameters. The argument is a new unique list of elements of the input * set. Despite being in a vector, these values should really be treated as a set. Values * will be in the same order as encountered in the input iterator. */  Jason Rhinelander committed Nov 03, 2015 40 template  Jason Rhinelander committed May 23, 2017 41 std::enable_if_t::iterator_category>::value>  Jason Rhinelander committed May 31, 2013 42 43 44 all_combinations( const It &begin, const It &end,  Jason Rhinelander committed Nov 03, 2015 45  std::function &)> func  Jason Rhinelander committed May 31, 2013 46 47 48  ) { // Store the current combination being considered  Jason Rhinelander committed Nov 03, 2015 49  std::vector combination; // Will store the current combination being considered  Jason Rhinelander committed May 31, 2013 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81  // The first combination is always the empty set; call with it: func(combination); std::stack it_stack; if (begin != end) it_stack.push(begin); while (!it_stack.empty()) { auto it = it_stack.top(); combination.push_back(*it); // Call the caller-provided function func(combination); auto n = std::next(it); if (n != end) { it_stack.push(n); } else { // We're at the end, so we need to pop ourselves and our value off the stack it_stack.pop(); combination.pop_back(); // Also pop off the previous value and increment it (the next value is going to be // pushed back on in the next iteration). combination.pop_back(); if (!it_stack.empty()) ++it_stack.top(); } } }  Jason Rhinelander committed Jan 12, 2015 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 /** Transforms the range [first, last) into the next strictly-increasing sequence with maximum * value max. To cover all values, the range should initially be in sorted order from min to * min+n, where n is the size of the range. Each permutation will be in sorted, * strictly-increasing order. The last permutation will consist of {max-n, max-n+1, ..., max-1, * max}, where n is the size of the given range. * * Note that this algorithm does not check that values are sorted, so may behave in an undefined * manner if called with initially-unsorted values. * * \returns true if the range was updated to the next permutation, false if no next permutation * exists. * * For example: * * std::vector v({1,2,3}); * do { ... } * while (eris::next_sorted_int_combination(v.begin(), v.end(), 5)); * * will execute the ... code 10 times with v set to: {1,2,3}, {1,2,4}, {1,2,5}, {1,3,4}, * {1,3,5}, {1,4,5}, {2,3,4}, {2,3,5}, {2,4,5}, {3,4,5}. */  Jason Rhinelander committed May 23, 2017 103 104 template ::value &&  Jason Rhinelander committed Jan 12, 2015 105  std::is_base_of::iterator_category>::value  Jason Rhinelander committed May 23, 2017 106 107 , int> = 0> bool next_increasing_integer_permutation(BidirIt first, BidirIt last, typename BidirIt::value_type max) {  Jason Rhinelander committed Jan 12, 2015 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125  auto it = last; --it; while (true) { if (*it < max) { auto last_val = ++*it; for (++it; it != last; ++it) { *it = ++last_val; } return true; } if (it == first) break; --it; --max; } return false; }  Jason Rhinelander committed Jul 19, 2013 126 127 128 129 130 131 /** Generic class for a stepping a value up or down by a certain amount, increasing or decreasing * the step size based on the previous steps. This is often used in optimizers to find an optimal * output/price level. */ class Stepper final { public:  Jason Rhinelander committed Apr 04, 2014 132  /// The default (possibly relative) initial step  Jason Rhinelander committed Jul 25, 2013 133  static constexpr double default_initial_step = 1.0/32.0;  Jason Rhinelander committed Apr 04, 2014 134  /// The default number of same-direction steps before the step size is increased  Jason Rhinelander committed Jul 25, 2013 135  static constexpr int default_increase_count = 4;  Jason Rhinelander committed Apr 04, 2014 136  /// The smallest (possibly relative) step that will be taken  Jason Rhinelander committed Jul 25, 2013 137  static constexpr double default_min_step = std::numeric_limits::epsilon();  Jason Rhinelander committed Apr 04, 2014 138  /// The largest (possibly relative) step that will be taken  Jason Rhinelander committed Jan 31, 2014 139  static constexpr double default_max_step = 0.5;  Jason Rhinelander committed Apr 04, 2014 140  /// The default value for whether steps are relative (true) or absolute (false)  Jason Rhinelander committed Aug 20, 2013 141  static constexpr bool default_relative_steps = true;  Jason Rhinelander committed Jul 25, 2013 142   Jason Rhinelander committed Jul 19, 2013 143 144  /** Constructs a new Stepper object. *  Jason Rhinelander committed Jul 25, 2013 145  * \param initial_step the initial size of a step, relative to the current value. Defaults to 1/32  Jason Rhinelander committed Jul 19, 2013 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162  * (about 3.1%). Over time the step size will change based on the following options. When * increasing, the new value is \f$(1 + step)\f$ times the current value; when decreasing * the value is \f$\frac{1}{1 + step}\f$ times the current value. This ensures that an * increase followed by a decrease will result in the same value (at least within numerical * precision). * * \param increase_count if the value moves consistently in one direction this many times, * the step size will be doubled. Defaults to 4 (i.e. the 4th step will be doubled). Must * be at least 2, though values less than 4 aren't recommended in practice. When * increasing, the previous increase count is halved; thus, with the default value of 4, it * takes only two additional steps in the same direction before increasing the step size * again. With a value of 6, it would take 6 changes for the first increase, then 3 for * subsequent increases, etc. * * \param min_step is the minimum step size that can be taken. The default is equal to * machine epsilon (i.e. the smallest value v such that 1 + v is a value distinct from 1), * which is the smallest value possible for a step.  Jason Rhinelander committed Aug 20, 2013 163  *  Jason Rhinelander committed Jan 31, 2014 164 165 166  * \param max_step is the maximum step size that can be taken. The default is equal to 0.5, * corresponding to an an increase of 50% or a decrease of 33% (when rel_steps is true). *  Jason Rhinelander committed Aug 20, 2013 167 168 169 170 171 172 173  * \param rel_steps specifies whether the steps taken are relative to the current value, or * absolute changes. The default (relative steps) is suitable for things like prices and * quantities, when relative changes are more important than absolute values; specifying * false here is suitable for spatial models, when movement shouldn't be affected by the * current variable value. Note that the default min_step is too small to numerically * affect values much below -1.0 or above +1.0, so choosing a more appropriate min_step is * advised.  Jason Rhinelander committed Jul 19, 2013 174  */  Jason Rhinelander committed Jul 25, 2013 175 176  Stepper(double initial_step = default_initial_step, int increase_count = default_increase_count,  Jason Rhinelander committed Aug 20, 2013 177  double min_step = default_min_step,  Jason Rhinelander committed Jan 31, 2014 178  double max_step = default_max_step,  Jason Rhinelander committed Aug 20, 2013 179  bool rel_steps = default_relative_steps);  Jason Rhinelander committed Jul 19, 2013 180 181  /** Called to signal that a step increase (up=true) or decrease (up=false) should be  Jason Rhinelander committed Aug 20, 2013 182 183 184 185  * taken. * * If relative_steps is true, this returns the relative step value, where 1 indicates no * change, 1.2 indicates a 20% increase, etc. The relative step multiple will always be a  Jason Rhinelander committed Feb 20, 2014 186 187 188  * strictly positive value not equal to 1. If \f$s\f$ is the current step size, the * returned value will be either \f$1+s\f$ for an upward step or \f$\frac{1}{1+s}\f$ for a * downward step.  Jason Rhinelander committed Jul 19, 2013 189  *  Jason Rhinelander committed Aug 20, 2013 190 191 192  * If relative_steps is false, this returns the absolute change in the current value, which * could be any positive or negative value; the returned value is the amount by which the * variable should be changed.  Jason Rhinelander committed Jul 19, 2013 193 194 195 196 197  * * When called, this first checks whether the step size should be changed: if at least * increase steps in the same direction have occured, the step size is doubled (and the * count of previous steps in this direction is halved); if the last step was in the * opposite direction, the step size is halved.  Jason Rhinelander committed Feb 20, 2014 198  *  Jason Rhinelander committed Jul 23, 2015 199  * After this is called, you may optionally consider the oscillating_min value, which will  Jason Rhinelander committed Feb 20, 2014 200 201  * tell you how many of the previous steps have simply oscillated around the minimum step * size (and thus aren't doing anything useful).  Jason Rhinelander committed Jul 19, 2013 202 203 204 205  */ double step(bool up); /// The number of steps in the same direction required to double the step size  Jason Rhinelander committed Feb 20, 2014 206  int increase;  Jason Rhinelander committed Jul 19, 2013 207   Jason Rhinelander committed Feb 20, 2014 208 209  /// The minimum (relative) step size allowed, as specified in the constructor. double min_step;  Jason Rhinelander committed Jul 19, 2013 210   Jason Rhinelander committed Feb 20, 2014 211 212  /// The maximum (relative) step size allowed, as specified in the constructor. double max_step;  Jason Rhinelander committed Jan 31, 2014 213   Jason Rhinelander committed Aug 20, 2013 214 215 216 217 218 219 220 221 222 223 224 225 226 227  /** If true, steps are relative; if false, steps are absolute. * * For example, a relative step size of 1/10 results in a step to either 11/10 or 10/11 for * an upward or downward step, respectively. Thus consecutive downward steps approach but * never reach zero. For variables that have a positive domain (e.g. prices and * quantities), this is exactly what is wanted, because typically relative changes are more * important than absolute changes. * * For some problems, however, steps should be considered as absolute. For example, in * a linear spatial problem, a step of 0.1 should be +0.1 or -0.1; it makes no sense to * scale steps by the current position (since that is mathematically arbitrary). Thus, * when this value is false, the value passed to take_step will be the absolute change, * *not* a multiple of the current value. */  Jason Rhinelander committed Feb 20, 2014 228  bool relative_steps;  Jason Rhinelander committed Aug 20, 2013 229   Jason Rhinelander committed Feb 20, 2014 230  /// The most recent step size. If no step has been taken yet, this is the initial step.  Jason Rhinelander committed Jul 19, 2013 231 232 233 234 235  double step_size; /// The most recent step direction: true if the last step was positive, false if negative. bool prev_up = true;  Jason Rhinelander committed Feb 20, 2014 236 237 238 239 240 241 242  /** Will be set to the number of times the step direction has oscillated back and forth * while at the minimum step size, and thus the normal action of reducing the step size * can't be taken. As soon as two steps occur in the same directory, this will be reset to * 0. */ unsigned int oscillating_min = 0;  Jason Rhinelander committed Jul 19, 2013 243 244 245 246 247  /** The number of steps that have occurred in the same direction. When a step size doubling * occurs, this value is halved (since the previous steps are only half the size of current * steps). */ int same = 0;  Jason Rhinelander committed Aug 20, 2013 248   Jason Rhinelander committed Jul 19, 2013 249 250 };  Jason Rhinelander committed Aug 06, 2017 251 252 /// Struct holding the results of a call to an optimization function such as single_peak_search() or /// constrained_maximum_search()  Jason Rhinelander committed Oct 31, 2017 253 254 255 template struct search_result { domain_t arg; ///< The argument that maximizes the searched function value_t value; ///< The value of the function at .arg  Jason Rhinelander committed Aug 06, 2017 256 257 258 259 260 261  /** 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;  Jason Rhinelander committed Oct 31, 2017 262  operator domain_t() const { return arg; } ///< Implicit conversion to double returns .arg  Jason Rhinelander committed Aug 06, 2017 263   Jason Rhinelander committed Oct 31, 2017 264  search_result(domain_t a, value_t v, bool in, int it = -1) :  Jason Rhinelander committed Aug 21, 2017 265  arg(std::move(a)), value(std::move(v)), inside(in), iterations(it) {}  Jason Rhinelander committed Aug 06, 2017 266 267 };  Jason Rhinelander committed Nov 24, 2017 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 template using non_deduced = std::common_type_t; /** Special class to pass a tolerance into one of the following search functions. This is usually * constructed via a call to either absolute_tolerance, to relative_tolerance, or by implicit * conversion from double (which is equivalent to relative_tolerance). */ template struct search_tolerance { const bool is_relative; const double relative; const AbsTol_t absolute; /** Implicit conversion from a double gives relative tolerance. Negative values are replaced * with 0. */ search_tolerance(double rel) : is_relative{true}, relative{rel > 0. ? rel : 0.}, absolute{0} {} /** Constructs an absolute tolerance; this is usually invoked via the absolute_tolerance() * function. */ search_tolerance(AbsTol_t abs, bool /* unused */) : is_relative{false}, relative{0}, absolute{std::move(abs)} {} search_tolerance(const search_tolerance &) = default; search_tolerance &operator=(const search_tolerance &) = default; search_tolerance(search_tolerance &&) = default; search_tolerance &operator=(search_tolerance &&) = default; /** Implicit conversion to a tolerance with a different domain type; this casts the absolute * value from the foreign to the local type; it is only allowed when the absolute type is * convertible. */ template ::value, int> = 0> search_tolerance(const search_tolerance &tol) : is_relative{tol.is_relative}, relative{tol.relative}, absolute{tol.absolute} {} /** Implicit conversion to a type with a different domain type with non-convertible domain * types; this is only allowed if the copied-from value is a relative tolerance instance (and * throws a domain_error if not). */ template ::value, int> = 0> search_tolerance(const search_tolerance &tol) : is_relative{tol.is_relative}, relative{tol.relative}, absolute{0} { if (!is_relative) throw std::domain_error("Cannot cast absolute search_tolerance types"); } }; /// Constructs a tolerance object that specifies absolute tolerance. template search_tolerance absolute_tolerance(domain_t tol) { return search_tolerance(tol, true); } /// Constructs a tolerance object that specifies relative tolerance. inline search_tolerance relative_tolerance(double tol) { return tol; }  Jason Rhinelander committed Aug 06, 2017 315 316 317 318 319 320 321 322 /// The constant phi. Callers can specialize this template if using custom types with more /// precision than a long double value. template constexpr RealType phi = 1.61803398874989484820458683436563811L; /// The right inner point multiple for a golden section search, phi - 1. template constexpr RealType golden_section_right = phi - RealType(1); /// The left inner point multiple for a golden section search, 1 - right, which also equals 2 - phi. template constexpr RealType golden_section_left = RealType(1) - golden_section_right;  Jason Rhinelander committed Oct 04, 2016 323   Jason Rhinelander committed Apr 29, 2014 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 /** 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). *  Jason Rhinelander committed Aug 06, 2017 342 343  * \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.  Jason Rhinelander committed Apr 29, 2014 344 345  * \param left the left edge of the domain to consider * \param right the right edge of the domain to consider  Jason Rhinelander committed May 09, 2014 346 347  * \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  Jason Rhinelander committed Aug 06, 2017 348 349 350  * 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  Jason Rhinelander committed Oct 31, 2017 351 352 353 354  * maximum precision). * * The domain type may be optionally specified as a template arugment; if not provided it defaults * to double.  Jason Rhinelander committed Oct 04, 2016 355  *  Jason Rhinelander committed Aug 06, 2017 356 357  * \return a eris::search_result struct with .arg set to the peak argument, and .max set to * the value of the function at .arg.  Jason Rhinelander committed Apr 29, 2014 358  */  Jason Rhinelander committed Oct 31, 2017 359 360 361 template ()(std::declval()))> search_result single_peak_search( Func f,  Jason Rhinelander committed Nov 24, 2017 362 363 364  non_deduced left, non_deduced right, search_tolerance> tolerance = 1e-10) {  Jason Rhinelander committed Aug 06, 2017 365   Jason Rhinelander committed Oct 31, 2017 366 367  constexpr domain_t midpoint_right = phi - domain_t(1); constexpr domain_t midpoint_left = domain_t(1) - midpoint_right;  Jason Rhinelander committed Aug 06, 2017 368 369 370 371  // Track whether the peak is strictly inside the initial boundaries: bool inside_left = false, inside_right = false;  Jason Rhinelander committed Oct 31, 2017 372 373 374 375  domain_t span = right - left; domain_t midleft = left + midpoint_left * span; domain_t midright = left + midpoint_right * span; value_t fl = f(left), fml = f(midleft), fmr = f(midright), fr = f(right);  Jason Rhinelander committed Aug 06, 2017 376 377 378 379  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;  Jason Rhinelander committed Nov 24, 2017 380 381  bool done; do {  Jason Rhinelander committed Aug 06, 2017 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413  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); }  Jason Rhinelander committed Nov 24, 2017 414 415 416 417 418  done = tolerance.is_relative ? span <= tolerance.relative * max(abs(left), abs(right)) : span <= tolerance.absolute; } while (!done);  Jason Rhinelander committed Aug 06, 2017 419 420 421 422 423 424 425 426 427 428 429 430  // 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}; }  Jason Rhinelander committed Oct 04, 2016 431   Jason Rhinelander committed Nov 13, 2017 432 433 434 435 436 /// Divides by 2, but with specializations for floating point types that do so by multiplying by /// 0.5 instead. template ::value, int> = 0> T half(T val) { return val * T(0.5); } template ::value, int> = 0> T half(T val) { return val / T(2); }  Jason Rhinelander committed Nov 24, 2017 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 /** Class that specifies a RHS limit when implicitly converted from a domain value, or that * specifies that the RHS limit should be found automatically when default constructed. * * This is not typically invoked directly: instead either specify a value for the right argument, * or specify search_right(). */ template struct search_right_val { /// True if right-hand-side value should be found automatically bool search = true; /// The right value (if search is false). domain_t right; /// Default constructor: specifies a right argument that should be found automatically. search_right_val() = default; /// Implicit conversion from a domain_t value: uses the specified right value directly /// (without searching). If the type supports NaN and the given value is NaN, search mode is /// enabled. search_right_val(domain_t right) : search{false}, right{right} { using std::isnan; if (std::numeric_limits::has_quiet_NaN() && isnan(right)) search = true; } /** Implicit conversion to a search_right_val with a different domain type; this casts the * right value from the foreign to the local type; it is only allowed when the value type is * convertible. */ template ::value, int> = 0> search_right_val(const search_right_val &srv) : search{srv.search}, right{srv.right} {} /** Implicit conversion to a type with a different domain type with non-convertible domain * types; this is only allowed if the copied-from value is a search = true instance (and * throws a domain_error if not). */ template ::value, int> = 0> search_right_val(const search_right_val &srv) : search{srv.search}, right{0} { if (!search) throw std::domain_error("Cannot cast absolute search_right_val types"); } }; /** Constructs a search_right_val that searches for the right-hand side. */ inline search_right_val search_right() { return {}; }  Jason Rhinelander committed Aug 06, 2017 481 482 483 484 485 486 487 488 489 /** Performs a binary search to find the maximum function value that satisfies a constraint given * a pair of values that satisfy and do not satisfy the constraint. * * This algorithm works by evaluating the function at the midpoint between left and right, * eliminating either the left half or right half of the interval at each step depending on whether * the constraint is satisfied or not at the mid-point. * * \param f the function (or function-like) object that can be called with a single double argument * and returns true if the value satisfies the constraint and false otherwise.  Jason Rhinelander committed Nov 13, 2017 490 491 492 493 494  * * \param left the left edge of the domain to consider; the constraint should be satisfied at * f(left) (if not, the algorithms fails immediately). * * \param right the right edge of the domain to consider at which the constraint should not be  Jason Rhinelander committed Nov 24, 2017 495 496 497 498 499 500 501 502 503 504 505 506  * satisfied. Can be specified as a regular value, or as a special search_right() value: if the * latter, the algorithm first starts at x = max{-left, 2*left}, doubling x until a constraint * violation is encountered, at which point this x becomes right. (This attempts to determine * right are not counted in iterations in the returned object). If left is 0, this starts * looking at right at 1.0 (then doubling, etc.). * * \param tolerance the domain value tolerance at which to stop searching: the algoritm proceeds * until the difference between the left and right edge has been narrowed to this value or less. * This is typically specified as a double value for relative tolerance (relative to the larger * absolute value of right or left), or a call to absolute_tolerance(...) for absolute tolerance. * The default, if unspecified, is relative tolerance of 1e-10. It is safe to pass a value of 0 * here: the algorithm will then run to the maximum precision of the domain type.  Jason Rhinelander committed Aug 06, 2017 507 508 509 510 511 512 513 514 515  * * \return a eris::search_result struct. If the initial f(left) is not satisfied, this * immediately returns (with value set to false, inside set to false, and arg set to left). * If the initial f(right) is satisfied, this immediately returns with value set to true, * inside set to false, and arg set to right. Otherwise (i.e. for an interior solution) this * returns with .arg set to the maximum constraint-satisfying argument, .value set to true, and * .inside set to true. (Note that a maximum at exactly left is considered "inside" [left, * right)). */  Jason Rhinelander committed Nov 24, 2017 516 template  Jason Rhinelander committed Oct 31, 2017 517 search_result constrained_maximum_search(  Jason Rhinelander committed Nov 24, 2017 518 519 520 521  Func f, non_deduced left, search_right_val> &&s_right, search_tolerance> tolerance = 1e-10) {  Jason Rhinelander committed Aug 06, 2017 522  static_assert(std::is_same::value,  Jason Rhinelander committed Nov 13, 2017 523  "constrained_maximum_search: given function must return a bool value");  Jason Rhinelander committed Aug 06, 2017 524 525 526 527  bool fl = f(left); if (!fl) return {std::move(left), false, false, 0};  Jason Rhinelander committed Aug 21, 2017 528 529 530 531 532 533  // Don't call std::abs etc. directly (to allow ADL on these funcs) using std::abs; using std::max; using std::isfinite; using std::isnan;  Jason Rhinelander committed Nov 24, 2017 534 535  domain_t right; if (s_right.search) {  Jason Rhinelander committed Nov 13, 2017 536  domain_t x = left < domain_t(0) ? -left : left > domain_t(0) ? domain_t(2)*left : domain_t(1);  Jason Rhinelander committed Nov 24, 2017 537 538 539 540 541  while (isfinite(x) && f(x)) { left = x; fl = true; x *= domain_t(2); }  Jason Rhinelander committed Aug 21, 2017 542 543  right = x; }  Jason Rhinelander committed Nov 24, 2017 544 545  else right = s_right.right;  Jason Rhinelander committed Aug 21, 2017 546   Jason Rhinelander committed Aug 06, 2017 547  bool fr = f(right);  Jason Rhinelander committed Aug 21, 2017 548  if (fr || !isfinite(right)) return {std::move(right), fr, false, 0};  Jason Rhinelander committed Aug 06, 2017 549   Jason Rhinelander committed Oct 31, 2017 550  domain_t span = right - left;  Jason Rhinelander committed Aug 06, 2017 551  int iterations = 0;  Jason Rhinelander committed Nov 24, 2017 552 553  bool done; do {  Jason Rhinelander committed Aug 06, 2017 554  iterations++;  Jason Rhinelander committed Nov 13, 2017 555  domain_t mid = left + half(span);  Jason Rhinelander committed Aug 06, 2017 556 557 558 559 560 561 562 563 564 565 566 567  if (mid == left || mid == right) break; // Numerical precision limit bool fm = f(mid); if (fm) { // The midpoint satisfies the constraint, so throw away the left half left = std::move(mid); fl = fm; } else { right = std::move(mid); fr = fm; } span = right - left;  Jason Rhinelander committed Nov 24, 2017 568 569 570 571 572  done = tolerance.is_relative ? span <= tolerance.relative * max(abs(left), abs(right)) : span <= tolerance.absolute; } while (!done);  Jason Rhinelander committed Aug 06, 2017 573 574 575 576 577 578  return {std::move(left), true, true, iterations}; } /// Helper that converts arguments to use constrained_maximum_search to find the smallest /// satisfying argument rather than the largest.  Jason Rhinelander committed Nov 13, 2017 579 580 581 582 583 /// /// The constraint should be satisfied at f(right) and *not* satisfied at f(left). As with /// constrained_maximum_search, right can be specified as NaN, in which case the same search for /// an initial right value will be done as in constrained_maximum_search, but looking for an /// initial right value that doesn't satisfy the constraint.  Jason Rhinelander committed Nov 24, 2017 584 template  Jason Rhinelander committed Oct 31, 2017 585 search_result constrained_minimum_search(  Jason Rhinelander committed Nov 24, 2017 586 587 588 589  Func f, non_deduced left, search_right_val> &&s_right, search_tolerance> tolerance = 1e-10) {  Jason Rhinelander committed Aug 06, 2017 590 591  static_assert(std::is_same::value, "constrained_minimum_search: given function must return a bool value");  Jason Rhinelander committed Nov 13, 2017 592   Jason Rhinelander committed Nov 24, 2017 593 594  domain_t right; if (s_right.search) {  Jason Rhinelander committed Nov 13, 2017 595  domain_t x = left < domain_t(0) ? -left : left > domain_t(0) ? domain_t(2)*left : domain_t(1);  Jason Rhinelander committed Nov 24, 2017 596 597 598 599  while (isfinite(x) && !f(x)) { left = x; x *= domain_t(2); }  Jason Rhinelander committed Nov 13, 2017 600 601  right = x; }  Jason Rhinelander committed Nov 24, 2017 602 603  else right = s_right.right;  Jason Rhinelander committed Nov 13, 2017 604   Jason Rhinelander committed Aug 06, 2017 605  auto ret = constrained_maximum_search(  Jason Rhinelander committed Oct 31, 2017 606  [f = std::move(f)](domain_t a) { return f(-a); },  Jason Rhinelander committed Nov 24, 2017 607  -right, -left, tolerance);  Jason Rhinelander committed Nov 13, 2017 608  ret.arg = -ret.arg;  Jason Rhinelander committed Aug 06, 2017 609 610 611  return ret; }  Jason Rhinelander committed May 31, 2013 612 }