bisection.hh

#ifndef VTESTBED_CORE_BISECTION_HH
#define VTESTBED_CORE_BISECTION_HH

#include <vtestbed/base/blitz.hh>

//#define VTB_DEBUG_BISECTION

namespace vtb {

    namespace core {

        template <class T, int N, class F>
        blitz::TinyVector<T,N>
        bisection(
            blitz::TinyVector<T,N> a,
            blitz::TinyVector<T,N> b,
            F func,
            const blitz::TinyVector<T,N>& arg_eps,
            const blitz::TinyVector<T,N>& func_eps,
            int max_iter
        ) {
            typedef blitz::TinyVector<T,N> vec;
            vec fa{func(a)};
            if (all(abs(fa) < func_eps)) {
                #if defined(VTB_DEBUG_BISECTION)
                std::clog
                    << __func__ << ": "
                    << "a=" << a
                    << ",b=" << b
                    << ",fa=" << fa
                    << '\n';
                #endif
                return a;
            }
            vec fb{func(b)};
            if (all(abs(fb) < func_eps)) {
                #if defined(VTB_DEBUG_BISECTION)
                std::clog
                    << __func__ << ": "
                    << "a=" << a
                    << ",b=" << b
                    << ",fb=" << fb
                    << '\n';
                #endif
                return b;
            }
            vec c, fc;
            c = 0, fc = 0;
            for (int i=1; i<=max_iter; ++i) {
                c = T{0.5}*(a + b);
                fc = func(c);
                for (int j=0; j<N; ++j) {
                    if (fa(j)*fc(j) < 0) {
                        b(j) = c(j);
                        fb(j) = fc(j);
                    } else if (fb(j)*fc(j) < 0) {
                        a(j) = c(j);
                        fa(j) = fc(j);
                    } else {
                        if (abs(fa(j)) < abs(fb(j))) {
                            c(j) = a(j);
                        } else {
                            c(j) = b(j);
                        }
                        //return c;
                    }
                }
                #if defined(VTB_DEBUG_BISECTION)
                std::clog
                    << __func__ << ": "
                    << "a=" << a
                    << ",b=" << b
                    << ",c=" << c
                    << ",fa=" << fa
                    << ",fb=" << fb
                    << ",fc=" << fc
                    << '\n';
                #endif
                if (all(abs(b-a) < arg_eps) && all(abs(fc) < func_eps)) {
                    break;
                }
            }
            return c;
        }

        template <class T, class F>
        T
        bisection(T a, T b, F func, T arg_eps, T func_eps, int max_iter) {
            using std::abs;
            T fa = func(a);
            if (abs(fa) < func_eps) {
                return a;
            }
            T fb = func(b);
            if (abs(fb) < func_eps) {
                return b;
            }
            T c{0}, fc{0};
            for (int i=0; i<max_iter; ++i) {
                c = T{0.5}*(a + b);
                fc = func(c);
                if (fa*fc < 0) {
                    b = c;
                    fb = fc;
                } else if (fb*fc < 0) {
                    a = c;
                    fa = fc;
                } else {
                    if (abs(fa) < abs(fb)) {
                        c = a;
                    } else {
                        c = b;
                    }
                    break;
                }
                #if defined(VTB_DEBUG_BISECTON)
                const int w = 12;
                std::clog
                    << std::setw(w) << i
                    << std::setw(w) << a
                    << std::setw(w) << b
                    << std::setw(w) << fa
                    << std::setw(w) << fb
                    << std::setw(w) << fc
                    << std::endl;
                #endif
                if (abs(b-a) < arg_eps && abs(fc) < func_eps) {
                    break;
                }
            }
            return c;
        }

        template <class T, int N>
        class Bisection {

        private:
            using vec = blitz::TinyVector<T,N>;

        private:
            vec _argeps{T{1e-3}};
            vec _funceps{T{1e-3}};
            int _niterations = 100;

        public:

            Bisection() = default;

            inline
            Bisection(
                const vec& argeps,
                const vec& funceps,
                int niterations
            ):
            _argeps(argeps),
            _funceps(funceps),
            _niterations(niterations)
            {}

            template <class Func>
            vec
            operator()(const vec& lbound, const vec& ubound, Func func) const {
                return bisection<T,N,Func>(
                    lbound,
                    ubound,
                    func,
                    this->_argeps,
                    this->_funceps,
                    this->_niterations
                );
            }

            inline void
            argument_precision(const vec& eps) {
                this->_argeps = eps;
            }

            inline void
            function_precision(const vec& eps) {
                this->_funceps = eps;
            }

            inline int max_iterations() const noexcept { return this->_niterations; }
            inline void max_iterations(int n) { this->_niterations = n; }

        };

        template <class T>
        class Bisection<T,1> {

        private:
            T _argeps{1e-3};
            T _funceps{1e-3};
            int _niterations = 100;

        public:

            Bisection() = default;

            inline
            Bisection(const T& argeps, const T& funceps, int niterations):
            _argeps(argeps), _funceps(funceps), _niterations(niterations) {}

            template <class Func>
            T
            operator()(const T& lbound, const T& ubound, Func func) const {
                return bisection<T,Func>(
                    lbound,
                    ubound,
                    func,
                    this->_argeps,
                    this->_funceps,
                    this->_niterations
                );
            }

            inline void
            argument_precision(T eps) {
                this->_argeps = eps;
            }

            inline void
            function_precision(T eps) {
                this->_funceps = eps;
            }

            inline int
            max_iterations() const noexcept {
                return this->_niterations;
            }

            inline void
            max_iterations(int n) {
                this->_niterations = n;
            }

        };

    }

}

#endif // vim:filetype=cpp