lipo_solver.hh

#ifndef VTESTBED_CORE_LIPO_SOLVER_HH
#define VTESTBED_CORE_LIPO_SOLVER_HH

#include <limits>
#include <random>
#include <vector>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/core/types.hh>

namespace vtb {

    namespace core {

        template <class T, int N>
        class LIPO_solver {

        public:
            using scalar_type = T;
            using argument_type = Vector<T,N>;
            using value_type = T;

        private:
            value_type _constant{T{}};
            value_type _funceps = 0;
            argument_type _argeps{T{}};
            int _max_evaluations = 1000;
            int _max_iterations = 1000;
            scalar_type _sample_probability{T{}};
            bool _estimate_constant = true;

        public:

            inline value_type constant() const { return this->_constant; }
            inline argument_type argument_precision() const { return this->_argeps; }
            inline value_type function_precision() const { return this->_funceps; }
            inline int max_evaluations() const { return this->_max_evaluations; }
            inline int max_iterations() const { return this->_max_iterations; }
            inline scalar_type sample_probability() const { return this->_sample_probability; }
            inline bool estimate_constant() const { return this->_estimate_constant; }
            inline void constant(value_type rhs) { this->_constant = rhs; }
            inline void argument_precision(argument_type rhs) { this->_argeps = rhs; }
            inline void function_precision(value_type rhs) { this->_funceps = rhs; }
            inline void max_evaluations(int rhs) { this->_max_evaluations = rhs; }
            inline void max_iterations(int rhs) { this->_max_iterations = rhs; }
            inline void sample_probability(scalar_type rhs) { this->_sample_probability = rhs; }
            inline void estimate_constant(bool rhs) { this->_estimate_constant = rhs; }

            template <class Func> inline argument_type
            solve(argument_type lbound, argument_type ubound, Func func) const {
                using distribution_type = std::uniform_real_distribution<scalar_type>;
                using std::min;
                using std::max;
                if (max_evaluations() <= 0) { return argument_type{T{}}; }
                // Choose random point.
                std::random_device dev;
                std::default_random_engine prng(dev());
                argument_type x;
                for (int i=0; i<N; ++i) {
                    distribution_type dist(lbound(i), ubound(i));
                    x(i) = dist(prng);
                }
                if (max_evaluations() <= 1) { return x; }
                std::vector<argument_type> points;
                std::vector<value_type> values;
                points.reserve(max_evaluations());
                points.emplace_back(x);
                values.reserve(max_evaluations());
                values.emplace_back(func(x));
                value_type max_value = values.front();
                int nevaluations = 1;
                std::bernoulli_distribution sample(sample_probability());
                value_type k = estimate_constant() ? value_type{scalar_type{}} : constant();
                for (int i=0; i<max_iterations(); ++i) {
                    // Choose new random point.
                    for (int j=0; j<N; ++j) {
                        distribution_type dist(lbound(j), ubound(j));
                        x(j) = dist(prng);
                    }
                    auto old_nevaluations = nevaluations;
                    if (sample(prng)) {
                        points.emplace_back(x);
                        values.emplace_back(func(x));
                        ++nevaluations;
                    } else {
                        // Find new upper bound.
                        value_type min_ubound{std::numeric_limits<scalar_type>::max()};
                        for (int j=0; j<nevaluations; ++j) {
                            value_type new_ubound = values[j] + k*distance(x,points[j]);
                            if (new_ubound < min_ubound) { min_ubound = new_ubound; }
                        }
                        // Store new point.
                        if (min_ubound >= max_value) {
                            value_type func_x = func(x);
                            if (func_x > max_value) { max_value = func_x; }
                            points.emplace_back(x);
                            values.emplace_back(func_x);
                            ++nevaluations;
                        }
                    }
                    if (old_nevaluations != nevaluations) {
                        if (nevaluations >= max_evaluations()) { break; }
                        // Update the constant.
                        if (estimate_constant()) {
                            argument_type points_n = points[nevaluations-1];
                            value_type values_n = values[nevaluations-1];
                            for (int j=0; j<nevaluations-1; ++j) {
                                scalar_type length = distance(points[j],points_n);
                                value_type new_k = values[j]-values_n;
                                if (length != 0) { new_k /= length; }
                                k = max(k, new_k);
                            }
                        }
                    }
                }
                int max_index = 0;
                for (int i=1; i<nevaluations; ++i) {
                    if (values[i] > values[max_index]) {
                        max_index = i;
                    }
                }
                return points[max_index];
            }

        };

    }

}

#endif // vim:filetype=cpp