golden_section_search.h

Go to the documentation of this file.

/* Copyright (c) 2008-2022 the MRtrix3 contributors.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 * Covered Software is provided under this License on an "as is"
 * basis, without warranty of any kind, either expressed, implied, or
 * statutory, including, without limitation, warranties that the
 * Covered Software is free of defects, merchantable, fit for a
 * particular purpose or non-infringing.
 * See the Mozilla Public License v. 2.0 for more details.
 *
 * For more details, see http://www.mrtrix.org/.
 */
 
#ifndef __math_golden_section_search_h__
#define __math_golden_section_search_h__
 
#include <cmath>
 
#include "progressbar.h"
#include "memory.h"
 
namespace MR
{
  namespace Math
  {
 
    template <class FunctionType, typename ValueType>
      ValueType golden_section_search (FunctionType& function,
                                       const std::string& message,
                                       ValueType min_bound,
                                       ValueType init_estimate,
                                       ValueType max_bound,
                                       ValueType tolerance = 0.01)  {
 
        std::unique_ptr<ProgressBar> progress (message.size() ? new ProgressBar (message) : nullptr);
 
        const ValueType g1 = 0.61803399, g2 = 1 - g1;
        ValueType x0 = min_bound, x1, x2, x3 = max_bound;
 
        if (abs(max_bound - init_estimate) > abs(init_estimate - min_bound)) {
          x1 = init_estimate;
          x2 = init_estimate + g2 * (max_bound - init_estimate);
        } else {
          x2 = init_estimate;
          x1 = init_estimate - g2 * (init_estimate - min_bound);
        }
 
        ValueType f1 = function(x1);
        ValueType f2 = function(x2);
 
        while (tolerance * (abs(x1) + abs(x2)) < abs(x3 - x0)) {
          if (f2 < f1) {
            x0 = x1;
            x1 = x2;
            x2 = g1 * x1 + g2 * x3;
            f1 = f2, f2 = function(x2);
          } else {
            x3 = x2;
            x2 = x1;
            x1 = g1 * x2 + g2 * x0;
            f2 = f1, f1 = function(x1);
          }
          if (progress)
            ++*progress;
        }
        return f1 < f2 ? x1 : x2;
      }
  }
  }
 
#endif