gp_bo.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <random>
#include <stdexcept>
#include <vector>
 
#include "pulsatrix/acquisition_functions.hpp"
#include "pulsatrix/gaussian_process.hpp"
#include "pulsatrix/search_space.hpp"
#include "pulsatrix/trial.hpp"
 
namespace pulsatrix {
 
enum class AcquisitionKind { ExpectedImprovement, ProbabilityOfImprovement, UpperConfidenceBound };
 
inline Configuration UnitCubeToConfiguration(const SearchSpace& space, const std::vector<float>& t) {
    if (t.size() != space.size()) {
        throw std::invalid_argument("UnitCubeToConfiguration: t.size() must equal space.size()");
    }
    Configuration config;
    for (size_t i = 0; i < space.parameters().size(); ++i) {
        const auto& spec = space.parameters()[i];
        float u = t[i];
        if (u < 0.0f || u > 1.0f) {
            throw std::invalid_argument("UnitCubeToConfiguration: every t value must be in [0, 1]");
        }
        switch (spec.kind) {
            case ParameterKind::Continuous:
                config[spec.name] = spec.lower + static_cast<double>(u) * (spec.upper - spec.lower);
                break;
            case ParameterKind::LogUniform: {
                double log_lower = std::log(spec.lower);
                double log_upper = std::log(spec.upper);
                config[spec.name] = std::exp(log_lower + static_cast<double>(u) * (log_upper - log_lower));
                break;
            }
            default:
                throw std::invalid_argument(
                    "UnitCubeToConfiguration: GP-BO supports only Continuous/LogUniform parameters "
                    "(Integer/Categorical are TPE's job)");
        }
    }
    return config;
}
 
template <typename ObjectiveFn, typename RNG>
std::vector<Trial> RunGPBOLoop(const SearchSpace& space, ObjectiveFn objective_fn,
                                size_t num_initial_random, size_t num_iterations,
                                AcquisitionKind acquisition, size_t num_candidates, RNG& rng) {
    if (num_initial_random == 0) {
        throw std::invalid_argument("RunGPBOLoop: num_initial_random must be >= 1");
    }
 
    std::vector<Trial> trials;
    std::vector<std::vector<float>> observed_points;
    std::vector<float> observed_values;
    std::uniform_real_distribution<float> unit_dist(0.0f, 1.0f);
 
    auto evaluate = [&](const std::vector<float>& t) {
        Configuration config = UnitCubeToConfiguration(space, t);
        Trial trial(config);
        double value = objective_fn(config);
        trial.RecordMetric("objective", value, 0);
        trials.push_back(std::move(trial));
        observed_points.push_back(t);
        observed_values.push_back(static_cast<float>(value));
    };
 
    for (size_t i = 0; i < num_initial_random; ++i) {
        std::vector<float> t(space.size());
        for (auto& v : t) {
            v = unit_dist(rng);
        }
        evaluate(t);
    }
 
    for (size_t iter = 0; iter < num_iterations; ++iter) {
        GaussianProcessRegressor gp(/*sigma_f=*/1.0f, /*length_scale=*/0.3f, /*noise_variance=*/1e-6f);
        gp.Fit(observed_points, observed_values);
 
        float best_value = *std::max_element(observed_values.begin(), observed_values.end());
 
        std::vector<float> best_candidate;
        double best_score = -std::numeric_limits<double>::infinity();
        for (size_t c = 0; c < num_candidates; ++c) {
            std::vector<float> candidate(space.size());
            for (auto& v : candidate) {
                v = unit_dist(rng);
            }
            auto posterior = gp.Predict(candidate);
            double score = 0.0;
            switch (acquisition) {
                case AcquisitionKind::ExpectedImprovement:
                    score = ExpectedImprovement(posterior.mean, posterior.variance, best_value);
                    break;
                case AcquisitionKind::ProbabilityOfImprovement:
                    score = ProbabilityOfImprovement(posterior.mean, posterior.variance, best_value);
                    break;
                case AcquisitionKind::UpperConfidenceBound:
                    score = UpperConfidenceBound(posterior.mean, posterior.variance);
                    break;
            }
            if (score > best_score) {
                best_score = score;
                best_candidate = candidate;
            }
        }
        evaluate(best_candidate);
    }
 
    return trials;
}
 
}  // namespace pulsatrix