brent.h

Go to the documentation of this file.

/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2014 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once
#if !defined(__MITSUBA_CORE_BRENT_H_)
#define __MITSUBA_CORE_BRENT_H_

#include <mitsuba/mitsuba.h>
#include <boost/function.hpp>

MTS_NAMESPACE_BEGIN

/**
 * \brief Brent's method nonlinear zero finder
 *
 * The implementation is transcribed from the Apache Commons
 * Java implementation. The supplied function is required to be
 * continuous, but not necessarily smooth.
 *
 * \ingroup libcore
 */
class MTS_EXPORT_CORE BrentSolver {
public:
        /// Return value of \ref BrentSolver::solve()
        struct Result {
                bool success;
                size_t iterations;
                Float x;
                Float y;

                /// Create a new result instance
                inline Result(bool success, size_t iterations, Float x, Float y)
                        : success(success), iterations(iterations), x(x), y(y) { }

                /// Return a string representation of the result
                inline std::string toString() const {
                        std::ostringstream oss;
                        oss << "BrentSolver::Result["
                                << "success=" << success << ", "
                                << "iterations=" << iterations << ", "
                                << "x=" << x << ", "
                                << "y=" << y << "]";
                        return oss.str();
                }
        };

        /**
         * \brief Create a new Brent-style solver with the
         * specified accuracy requirements
         *
         * \param maxIterations
         *      Max. number of successive iterations (default: 100)
         * \param absAccuracy
         *      Absolute accuracy requirement -- the iterations will stop
         *      when |f(x)| < absAccuracy.
         * \param absAccuracyPos
         *      Absolute accuracy requirement of the position -- the
         *      iterations will stop when |minX-maxX| < absAccuracyPos.
         * \param absAccuracyPos
         *      Absolute accuracy requirement of the position -- the
         *      iterations will stop when |minX-maxX|/minX < relAccuracyPos.
         */
        inline BrentSolver(size_t maxIterations = 100,
                        Float absAccuracy = 1e-6f,
                        Float absAccuracyPos = 1e-6f,
                        Float relAccuracyPos = 1e-6f)
                : m_maxIterations(maxIterations),
                  m_absAccuracy(absAccuracy),
                  m_absAccuracyPos(absAccuracyPos),
                  m_relAccuracyPos(relAccuracyPos) { }

        /**
         * \brief Find a zero in the given interval.
         *
         * Requires that the values of the function at the endpoints
         * have opposite signs.
         *
         * \param min the lower bound for the interval.
         * \param max the upper bound for the interval.
         * \return the value where the function is zero
         */
        Result solve(const boost::function<Float (Float)> &func,
                        Float min, Float max) const;

        /**
         * \brief Find a zero in the given interval with an initial guess
         *
         * Requires that the values of the function at the endpoints
         * have opposite signs (note that it is allowed to have endpoints
         * with the same sign if the initial point has opposite sign
         * function-wise).
         *
         * \param min the lower bound for the interval.
         * \param max the upper bound for the interval.
         * \param initial the start value to use (must be set to min
         *    if no initial point is known)
         * \return the value where the function is zero
         */
        Result solve(const boost::function<Float (Float)> &func,
                        Float min, Float max, Float initial) const;

        /**
         * Find a zero starting search according to the three provided points.
         *
         * \param x0 old approximation for the root
         * \param y0 function value at the approximation for the root
         * \param x1 last calculated approximation for the root
         * \param y1 function value at the last calculated approximation
         *    for the root
         * \param x2 bracket point (must be set to x0 if no bracket point is
         *   known, this will force starting with linear interpolation)
         * \param y2 function value at the bracket point.
         * \return the value where the function is zero
         */
        Result solve(const boost::function<Float (Float)> &func,
                                 Float x0, Float y0,
                                 Float x1, Float y1,
                                 Float x2, Float y2) const;
protected:
        size_t m_maxIterations;
        Float m_absAccuracy;
        Float m_absAccuracyPos;
        Float m_relAccuracyPos;
};

MTS_NAMESPACE_END

#endif /* __MITSUBA_CORE_BRENT_H_ */