diff --git a/albatross/covariance_functions/covariance_functions.h b/albatross/covariance_functions/covariance_functions.h index e3d7319c..1bc620b9 100644 --- a/albatross/covariance_functions/covariance_functions.h +++ b/albatross/covariance_functions/covariance_functions.h @@ -17,26 +17,27 @@ #include "noise.h" #include "polynomials.h" #include "radial.h" +#include "scaling_function.h" namespace albatross { template struct CovarianceFunction { - Term covariance_term; + Term term; - CovarianceFunction() : covariance_term(){}; - CovarianceFunction(const Term &term) : covariance_term(term){}; + CovarianceFunction() : term(){}; + CovarianceFunction(const Term &term_) : term(term_){}; template inline auto operator+(CovarianceFunction &other) { using Sum = SumOfCovarianceTerms; - auto sum = Sum(covariance_term, other.covariance_term); + auto sum = Sum(term, other.term); return CovarianceFunction{sum}; } template inline auto operator*(CovarianceFunction &other) { using Prod = ProductOfCovarianceTerms; - auto prod = Prod(covariance_term, other.covariance_term); + auto prod = Prod(term, other.term); return CovarianceFunction{prod}; } @@ -50,11 +51,11 @@ template struct CovarianceFunction { typename std::enable_if<(has_call_operator::value), int>::type = 0> inline auto operator()(const X &x, const Y &y) const { - return covariance_term(x, y); + return term(x, y); } /* - * If neither have a valid call method we fail compilation. + * If neither have a valid call method we fail. */ template ::value), @@ -67,23 +68,23 @@ template struct CovarianceFunction { * errors should give you an indication of which types were attempted. */ - inline auto get_name() const { return covariance_term.get_name(); }; - inline auto get_params() const { return covariance_term.get_params(); }; + inline auto get_name() const { return term.get_name(); }; + inline auto get_params() const { return term.get_params(); }; inline auto set_params(const ParameterStore ¶ms) { - return covariance_term.set_params(params); + return term.set_params(params); }; inline auto set_param(const ParameterKey &key, const ParameterValue &value) { - return covariance_term.set_param(key, value); + return term.set_param(key, value); }; - inline auto pretty_string() const { return covariance_term.pretty_string(); }; + inline auto pretty_string() const { return term.pretty_string(); }; inline auto get_params_as_vector() const { - return covariance_term.get_params_as_vector(); + return term.get_params_as_vector(); }; inline auto set_params_from_vector(const std::vector &x) { - return covariance_term.set_params_from_vector(x); + return term.set_params_from_vector(x); }; inline auto unchecked_set_param(const std::string &name, const double value) { - return covariance_term.unchecked_set_param(name, value); + return term.unchecked_set_param(name, value); }; }; diff --git a/albatross/covariance_functions/covariance_term.h b/albatross/covariance_functions/covariance_term.h index a4f66168..75be07f4 100644 --- a/albatross/covariance_functions/covariance_term.h +++ b/albatross/covariance_functions/covariance_term.h @@ -21,8 +21,13 @@ namespace albatross { /* - * An abstract class (though no purely abstract due to templating) - * which holds anything all Covariance terms should have in common. + * An abstract class which holds anything all Covariance + * terms should have in common. + * + * In addition to these abstract methods one or many + * methods with signature, + * operator ()(const X &x, const Y &y) + * should be defined. */ class CovarianceTerm : public ParameterHandlingMixin { public: @@ -135,7 +140,7 @@ class ProductOfCovarianceTerms : public CombinationOfCovarianceTerms { /* * If both LHS and RHS have a valid call method for the types X and Y - * this will return the sum of the two. + * this will return the product of the two. */ template ::value && diff --git a/albatross/covariance_functions/scaling_function.h b/albatross/covariance_functions/scaling_function.h new file mode 100644 index 00000000..9b82c046 --- /dev/null +++ b/albatross/covariance_functions/scaling_function.h @@ -0,0 +1,140 @@ +/* + * Copyright (C) 2018 Swift Navigation Inc. + * Contact: Swift Navigation + * + * This source is subject to the license found in the file 'LICENSE' which must + * be distributed together with this source. All other rights reserved. + * + * THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, + * EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + */ + +#ifndef ALBATROSS_COVARIANCE_FUNCTIONS_SCALING_FUNCTION_H +#define ALBATROSS_COVARIANCE_FUNCTIONS_SCALING_FUNCTION_H + +#include "covariance_term.h" +#include +#include + +namespace albatross { + +class ScalingFunction : public ParameterHandlingMixin { +public: + virtual std::string get_name() const = 0; + + // A scaling function should also implement operators + // for whichever types it is intended to scale using + // the signature: + // double operator(const X &x) const; +}; + +/* + * A scaling term is not actually a covariance function + * in the rigorous sense. It doesn't describe the uncertainty + * between variables, but instead operates deterministically + * on other uncertain variables. For instance, you may have + * some random variable, + * y ~ N(0, S) with S_ij = cov(y_i, y_j) + * And you may then make observations of that random variable + * but through a known transformation, + * z = f(y) * y + * where f is a determinstic function of y that returns a scalar. + * You might then ask what the covariance between two elements in + * z is which is woudl be given by, + * cov(z_i, z_j) = f(y_i) * cov(y_i, y_j) * f(y_j) + * but you might also be interested in the covariance between + * some y_i and an observation z_j, + * cov(y_i, z_j) = cov(y_i, y_j) * f(y_j) + * Here we see that for a typical covariance term, the covariance + * is only defined for two pairs of the same type, in this case + * operator()(Y &y, Y &y) + * but by multiplying by a ScalingTerm we end up with definitions + * for, + * operator()(Y &y, Z &z) + * which provides us with a way of computing the covariance between + * some hidden representation of a variable (y) and the actual + * observations (z) using a single determinstic mapping (f). + * + * This might be better explained by example which can be found + * in the tests (test_scaling_function). + */ +template class ScalingTerm : public CovarianceTerm { +public: + ScalingTerm() : CovarianceTerm(){}; + virtual ~ScalingTerm(){}; + + /* + * The following methods forward any requests dealing with + * the ParameterHandlingMixin to the ScalingFunction. + */ + std::string get_name() const override { return scaling_function_.get_name(); } + + std::string pretty_string() const { + return scaling_function_.pretty_string(); + } + + void set_params(const ParameterStore ¶ms) { + scaling_function_.set_params(params); + } + + virtual ParameterStore get_params() const { + return scaling_function_.get_params(); + } + + template void save(Archive &archive) const { + archive(cereal::make_nvp("base_class", + cereal::base_class(this))); + archive(cereal::make_nvp("scaling_function", scaling_function_)); + } + + template void load(Archive &archive) { + archive(cereal::make_nvp("base_class", + cereal::base_class(this))); + archive(cereal::make_nvp("scaling_function", scaling_function_)); + } + + void unchecked_set_param(const std::string &name, + const double value) override { + scaling_function_.set_param(name, value); + } + + /* + * If both Scaling and Covariance have a valid call method for the types X + * and Y this will return the product of the two. + */ + template < + typename X, typename Y, + typename std::enable_if<(has_call_operator::value && + has_call_operator::value), + int>::type = 0> + double operator()(X &x, Y &y) const { + return this->scaling_function_(x) * this->scaling_function_(y); + } + + /* + * If only one of the types has a scaling function we ignore the other. + */ + template ::value && + has_call_operator::value), + int>::type = 0> + double operator()(X &x, Y &y) const { + return this->scaling_function_(y); + } + + template < + typename X, typename Y, + typename std::enable_if<(has_call_operator::value && + !has_call_operator::value), + int>::type = 0> + double operator()(X &x, Y &y) const { + return this->scaling_function_(x); + } + +private: + ScalingFunction scaling_function_; +}; +} +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7cbf9c2e..2cc10ee6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,6 +14,7 @@ test_parameter_handling_mixin.cc test_models.cc test_core_distribution.cc test_tune.cc +test_scaling_function.cc ) add_dependencies(albatross_unit_tests diff --git a/tests/test_scaling_function.cc b/tests/test_scaling_function.cc new file mode 100644 index 00000000..6d005179 --- /dev/null +++ b/tests/test_scaling_function.cc @@ -0,0 +1,154 @@ +/* + * Copyright (C) 2018 Swift Navigation Inc. + * Contact: Swift Navigation + * + * This source is subject to the license found in the file 'LICENSE' which must + * be distributed together with this source. All other rights reserved. + * + * THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, + * EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + */ + +#include "covariance_functions/covariance_functions.h" +#include "test_utils.h" +#include +#include + +namespace albatross { + +/* + * This test data descibes the following situation. + * + * Assume you have some thin translucent sheet which attenuates + * the intensity of light passing through it but with an unknown + * attenuation. You then make some measurements of the attenuation + * by emitting light and measuring the intesity at different + * locations on the other side of the sheet, + * + * 0 emitter (at x = 1) + * ... + * . . . + * ====.==.==.===== <- thin sheet + * . . . + * . . . + * . . . + * ^ ^ ^ collectors + * + * 0------1------2 + * x axis -> + * + * In this case you have multiple measurements, y, each made at different + * locations, x. Each measurement passes through the sheet at a slightly + * different but deterministically known angle. We can then write + * the amount of observed attenuation, y, as + * + * y = f(x) * z + u + * + * where f(x) descibes how much longer the signal spends in the material + * as a function of the location, x, of the collector. z is the unknown + * attenuation and u is noise. + * + * The amount of material each ray passes through relative to the + * width is given by 1 / (cos(theta)) where theta = atan(x / 1) + * + */ + +double obliquity_function(double x) { return 1. / cos(atan(x - 1.)); } + +class ObliquityScaling : public ScalingFunction { +public: + ObliquityScaling() : ScalingFunction(){}; + + ~ObliquityScaling(){}; + + std::string get_name() const { return "obliquity_scaling"; } + + double operator()(const double &x) const { return obliquity_function(x); } +}; + +auto make_attenuation_data(const double attenuation = 3.14159, + const double sigma_noise = 0.01) { + std::random_device rd{}; + std::mt19937 gen{rd()}; + gen.seed(3); + std::normal_distribution<> d{0., sigma_noise}; + + s32 n = 10; + std::vector features; + Eigen::VectorXd targets(n); + + for (s32 i = 0; i < n; i++) { + // x ranges from 0 to 2 + double x = static_cast(i) * (2. / static_cast(n)); + features.push_back(x); + targets[i] = obliquity_function(x) * attenuation + d(gen); + } + + return RegressionDataset(features, targets); +} + +/* + * This test makes sure that we can make predictions of what + * the attenuation of a signal would be at some unobserved location. + */ +TEST(test_scaling_functions, test_predicts) { + using Feature = double; + using Noise = IndependentNoise; + + CovarianceFunction constant = {Constant(10.)}; + CovarianceFunction noise = {Noise(0.01)}; + using TestScalingTerm = ScalingTerm; + CovarianceFunction scaling = {TestScalingTerm()}; + + auto dataset = make_attenuation_data(); + + // This will create a covariance function that represents some constant + // unknown value, that is scaled according to a known deterministic + // function, then noisy measurements are taken. + auto covariance_function = constant * scaling + noise; + + auto model = gp_from_covariance(covariance_function); + + auto folds = leave_one_out(dataset); + auto cv_scores = + cross_validated_scores(root_mean_square_error, folds, &model); + + EXPECT_LE(cv_scores.mean(), 0.01); +} + +/* + * This test make sure (and illustrates how) we can perform inference + * on the unknown attenuation constant of the material in the test + * case described above. + */ +TEST(test_scaling_functions, test_inference) { + using Feature = double; + using Noise = IndependentNoise; + + double attenuation = 3.14159; + double sigma = 0.01; + + CovarianceFunction constant = {Constant(2 * attenuation)}; + CovarianceFunction noise = {Noise(sigma)}; + using TestScalingTerm = ScalingTerm; + CovarianceFunction scaling = {TestScalingTerm()}; + + auto dataset = make_attenuation_data(attenuation, sigma); + + // This will create a covariance function that represents some constant + // unknown value, that is scaled according to a known deterministic + // function, then noisy measurements are taken. + auto covariance_function = constant * scaling + noise; + + auto model = gp_from_covariance(covariance_function); + + auto state_space = + constant.term.get_state_space_representation(dataset.features); + model.fit(dataset); + auto state_estimate = model.inspect(state_space); + // Make sure our estimate of the attenuation term is close, despite the fact + // that we made scaled observations of it. + EXPECT_LE(fabs(state_estimate.mean[0] - attenuation), 1e-2); +} +}