From b853ec83973fc49b61959bd762bd9138c74b3bac Mon Sep 17 00:00:00 2001 From: Dan F-M Date: Mon, 15 Feb 2021 11:37:41 -0500 Subject: [PATCH] generating kernels separately --- .gitignore | 4 - MANIFEST.in | 4 + generate_kernels.py | 49 + pyproject.toml | 2 +- setup.py | 9 - src/george/include/george/kernels.h | 2707 +++++++++++++++++++++++++++ src/george/kernels.py | 1003 ++++++++++ 7 files changed, 3764 insertions(+), 14 deletions(-) create mode 100755 generate_kernels.py create mode 100644 src/george/include/george/kernels.h create mode 100644 src/george/kernels.py diff --git a/.gitignore b/.gitignore index 8222528..3be25da 100644 --- a/.gitignore +++ b/.gitignore @@ -12,8 +12,4 @@ build *.cpp .coverage src/george/george_version.py -src/george/kernels.py -src/george/kerneldefs.pxd -src/george/solvers/kerneldefs.pxd -src/george/include/george/kernels.h kernels/MyLocalGaussian.yml diff --git a/MANIFEST.in b/MANIFEST.in index 4e5e96d..40f3b9f 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,5 +1,9 @@ +graft vendor/eigen/Eigen + exclude .* prune .github prune docs prune paper +prune kernels +prune templates diff --git a/generate_kernels.py b/generate_kernels.py new file mode 100755 index 0000000..32f45ad --- /dev/null +++ b/generate_kernels.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import os +import glob +import yaml +from jinja2 import Template + + +def compile_kernels(fns): + template_dir = "templates" + output_dir = os.path.join("src", "george") + + with open(os.path.join(template_dir, "parser.h")) as f: + PARSER_TEMPLATE = Template(f.read()) + with open(os.path.join(template_dir, "kernels.h")) as f: + CPP_TEMPLATE = Template(f.read()) + with open(os.path.join(template_dir, "kernels.py")) as f: + PYTHON_TEMPLATE = Template(f.read()) + + specs = [] + for i, fn in enumerate(fns): + with open(fn, "r") as f: + spec = yaml.load(f.read(), Loader=yaml.FullLoader) + print("Found kernel '{0}'".format(spec["name"])) + spec["index"] = i + spec["reparams"] = spec.get("reparams", {}) + specs.append(spec) + print("Found {0} kernel specifications".format(len(specs))) + + fn = os.path.join(output_dir, "include", "george", "parser.h") + with open(fn, "w") as f: + print("Saving parser to '{0}'".format(fn)) + f.write(PARSER_TEMPLATE.render(specs=specs)) + fn = os.path.join(output_dir, "include", "george", "kernels.h") + with open(fn, "w") as f: + print("Saving C++ kernels to '{0}'".format(fn)) + f.write(CPP_TEMPLATE.render(specs=specs)) + fn = os.path.join(output_dir, "kernels.py") + with open(fn, "w") as f: + print("Saving Python kernels to '{0}'".format(fn)) + f.write(PYTHON_TEMPLATE.render(specs=specs)) + + +if __name__ == "__main__": + # If the kernel specifications are included (development mode) re-compile + # them first. + kernel_specs = glob.glob(os.path.join("kernels", "*.yml")) + compile_kernels(kernel_specs) diff --git a/pyproject.toml b/pyproject.toml index b00d51f..f22155b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=40.6.0", "wheel", "setuptools_scm[toml]", "PyYAML", "jinja2", "pybind11"] +requires = ["setuptools>=40.6.0", "wheel", "setuptools_scm[toml]", "pybind11"] build-backend = "setuptools.build_meta" [tool.setuptools_scm] diff --git a/setup.py b/setup.py index 2768c97..2d31a0b 100644 --- a/setup.py +++ b/setup.py @@ -52,15 +52,6 @@ def compile_kernels(fns): if not os.path.exists(os.path.join(localincl, "eigen", "Eigen", "Core")): raise RuntimeError("couldn't find Eigen headers") - # If the kernel specifications are included (development mode) re-compile - # them first. - kernel_specs = glob.glob(os.path.join("kernels", "*.yml")) - if len(kernel_specs): - print("Compiling kernels") - compile_kernels(kernel_specs) - if "kernels" in sys.argv: - sys.exit() - include_dirs = [ os.path.join("src", "george", "include"), os.path.join(localincl, "eigen"), diff --git a/src/george/include/george/kernels.h b/src/george/include/george/kernels.h new file mode 100644 index 0000000..add4088 --- /dev/null +++ b/src/george/include/george/kernels.h @@ -0,0 +1,2707 @@ +#ifndef _GEORGE_KERNELS_H_ +#define _GEORGE_KERNELS_H_ + +#include +#include +#include + + + +#include "george/metrics.h" +#include "george/subspace.h" + +#ifndef M_PI +#define M_PI 3.141592653589793238462643383279502884e+00 +#endif + +namespace george { + +namespace kernels { + +class Kernel { +public: + Kernel () {}; + virtual ~Kernel () {}; + virtual double value (const double* x1, const double* x2) { return 0.0; }; + virtual void gradient (const double* x1, const double* x2, + const unsigned* which, double* grad) {}; + virtual void x1_gradient (const double* x1, const double* x2, + double* grad) {}; + virtual void x2_gradient (const double* x1, const double* x2, + double* grad) {}; + + // Parameter vector spec. + virtual size_t size () const { return 0; } + virtual size_t get_ndim () const { return 0; } + virtual void set_parameter (size_t i, double v) {}; + virtual double get_parameter (size_t i) const { return 0.0; }; + virtual void set_metric_parameter (size_t i, double v) {}; + virtual void set_axis (size_t i, size_t v) {}; +}; + + +// +// OPERATORS +// + +class Operator : public Kernel { +public: + Operator (Kernel* k1, Kernel* k2) : kernel1_(k1), kernel2_(k2) {}; + ~Operator () { + delete kernel1_; + delete kernel2_; + }; + Kernel* get_kernel1 () const { return kernel1_; }; + Kernel* get_kernel2 () const { return kernel2_; }; + + // Parameter vector spec. + size_t size () const { return kernel1_->size() + kernel2_->size(); }; + size_t get_ndim () const { return kernel1_->get_ndim(); } + void set_parameter (size_t i, double value) { + size_t n = kernel1_->size(); + if (i < n) kernel1_->set_parameter(i, value); + else kernel2_->set_parameter(i-n, value); + }; + double get_parameter (size_t i) const { + size_t n = kernel1_->size(); + if (i < n) return kernel1_->get_parameter(i); + return kernel2_->get_parameter(i-n); + }; + +protected: + Kernel* kernel1_, * kernel2_; +}; + +class Sum : public Operator { +public: + Sum (Kernel* k1, Kernel* k2) : Operator(k1, k2) {}; + double value (const double* x1, const double* x2) { + return this->kernel1_->value(x1, x2) + this->kernel2_->value(x1, x2); + }; + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + size_t i, n1 = this->kernel1_->size(), n2 = this->size(); + + bool any = false; + for (i = 0; i < n1; ++i) if (which[i]) { any = true; break; } + if (any) this->kernel1_->gradient(x1, x2, which, grad); + + any = false; + for (i = n1; i < n2; ++i) if (which[i]) { any = true; break; } + if (any) this->kernel2_->gradient(x1, x2, &(which[n1]), &(grad[n1])); + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t ndim = this->get_ndim(); + std::vector g1(ndim), g2(ndim); + this->kernel1_->x1_gradient(x1, x2, &(g1[0])); + this->kernel2_->x1_gradient(x1, x2, &(g2[0])); + for (size_t i = 0; i < ndim; ++i) grad[i] = g1[i] + g2[i]; + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t ndim = this->get_ndim(); + std::vector g1(ndim), g2(ndim); + this->kernel1_->x2_gradient(x1, x2, &(g1[0])); + this->kernel2_->x2_gradient(x1, x2, &(g2[0])); + for (size_t i = 0; i < ndim; ++i) grad[i] = g1[i] + g2[i]; + }; +}; + +class Product : public Operator { +public: + Product (Kernel* k1, Kernel* k2) : Operator(k1, k2) {}; + double value (const double* x1, const double* x2) { + return this->kernel1_->value(x1, x2) * this->kernel2_->value(x1, x2); + }; + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + bool any; + size_t i, n1 = this->kernel1_->size(), n2 = this->size(); + double k; + + any = false; + for (i = 0; i < n1; ++i) if (which[i]) { any = true; break; } + if (any) { + k = this->kernel2_->value(x1, x2); + this->kernel1_->gradient(x1, x2, which, grad); + for (i = 0; i < n1; ++i) grad[i] *= k; + } + + any = false; + for (i = n1; i < n2; ++i) if (which[i]) { any = true; break; } + if (any) { + k = this->kernel1_->value(x1, x2); + this->kernel2_->gradient(x1, x2, &(which[n1]), &(grad[n1])); + for (i = n1; i < n2; ++i) grad[i] *= k; + } + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t ndim = this->get_ndim(); + std::vector g1(ndim), g2(ndim); + double k1 = this->kernel1_->value(x1, x2); + double k2 = this->kernel2_->value(x1, x2); + this->kernel1_->x1_gradient(x1, x2, &(g1[0])); + this->kernel2_->x1_gradient(x1, x2, &(g2[0])); + for (size_t i = 0; i < ndim; ++i) { + grad[i] = k2 * g1[i] + k1 * g2[i]; + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t ndim = this->get_ndim(); + std::vector g1(ndim), g2(ndim); + double k1 = this->kernel1_->value(x1, x2); + double k2 = this->kernel2_->value(x1, x2); + this->kernel1_->x2_gradient(x1, x2, &(g1[0])); + this->kernel2_->x2_gradient(x1, x2, &(g2[0])); + for (size_t i = 0; i < ndim; ++i) { + grad[i] = k2 * g1[i] + k1 * g2[i]; + } + }; +}; + + +/* +The linear regression kernel + +.. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = + \frac{(\mathbf{x}_i \cdot \mathbf{x}_j)^P}{\gamma^2} + +:param order: + The power :math:`P`. This parameter is a *constant*; it is not + included in the parameter vector. + +:param log_gamma2: + The scale factor :math:`\gamma^2`. +*/ + +class LinearKernel : public Kernel { +public: + LinearKernel ( + double log_gamma2, + double order, + size_t ndim, + size_t naxes + ) : + size_(1), + subspace_(ndim, naxes) + , param_log_gamma2_(log_gamma2) + , constant_order_(order) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_log_gamma2_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_log_gamma2_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double log_gamma2, + double inv_gamma2, + + double order, + double x1, double x2) { + if (order == 0.0) return inv_gamma2; + return pow(x1 * x2, order) * inv_gamma2; + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_log_gamma2_, + reparam_inv_gamma2_, + + constant_order_, + x1[j], x2[j]); + } + return value; + }; + + double log_gamma2_gradient ( + double log_gamma2, + double inv_gamma2, + + double order, + double x1, double x2) { + if (order == 0.0) return -inv_gamma2; + return -pow(x1 * x2, order) * inv_gamma2; + + }; + double _x1_gradient ( + double log_gamma2, + double inv_gamma2, + + double order, + double x1, double x2) { + if (order == 0.0) return 0.0; + return x2 * order * pow(x1 * x2, order - 1.0) * inv_gamma2; + + }; + + double _x2_gradient ( + double log_gamma2, + double inv_gamma2, + + double order, + double x1, double x2) { + if (order == 0.0) return 0.0; + return x1 * order * pow(x1 * x2, order - 1.0) * inv_gamma2; + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += log_gamma2_gradient( + param_log_gamma2_, + reparam_inv_gamma2_, + + constant_order_, + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_log_gamma2_, + reparam_inv_gamma2_, + + constant_order_, + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_log_gamma2_, + reparam_inv_gamma2_, + + constant_order_, + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_inv_gamma2_ = get_reparam_inv_gamma2 ( + param_log_gamma2_, + constant_order_ + ); + + }; + + double get_reparam_inv_gamma2 ( + double log_gamma2, + double order + ) { + return exp(-log_gamma2); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_log_gamma2_; + + double reparam_inv_gamma2_; + + double constant_order_; +}; + + +/* +You should always document your code.*/ + +class MyLocalGaussianKernel : public Kernel { +public: + MyLocalGaussianKernel ( + double x0, + double log_w, + size_t ndim, + size_t naxes + ) : + size_(2), + subspace_(ndim, naxes) + , param_x0_(x0) + , param_log_w_(log_w) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_x0_; + if (i == 1) return param_log_w_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_x0_ = value; + update_reparams(); + } else + if (i == 1) { + param_log_w_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double x0, + double log_w, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - x0, d2 = x2 - x0; + return exp(-(d1*d1 + d2*d2) * inv_2w); + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_x0_, + param_log_w_, + reparam_inv_2w_, + + x1[j], x2[j]); + } + return value; + }; + + double x0_gradient ( + double x0, + double log_w, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - x0, d2 = x2 - x0; + return 2 * exp(-(d1*d1 + d2*d2) * inv_2w) * inv_2w * (d1 + d2); + + }; + double log_w_gradient ( + double x0, + double log_w, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - x0, d2 = x2 - x0, + arg = (d1*d1 + d2*d2) * inv_2w; + return exp(-arg) * arg; + + }; + double _x1_gradient ( + double x0, + double log_w, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - x0, d2 = x2 - x0; + return -2.0 * exp(-(d1*d1 + d2*d2) * inv_2w) * d1 * inv_2w; + + }; + + double _x2_gradient ( + double x0, + double log_w, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - x0, d2 = x2 - x0; + return -2.0 * exp(-(d1*d1 + d2*d2) * inv_2w) * d2 * inv_2w; + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + grad[1] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += x0_gradient( + param_x0_, + param_log_w_, + reparam_inv_2w_, + + x1[j], x2[j]); + if (which[1]) + grad[1] += log_w_gradient( + param_x0_, + param_log_w_, + reparam_inv_2w_, + + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_x0_, + param_log_w_, + reparam_inv_2w_, + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_x0_, + param_log_w_, + reparam_inv_2w_, + + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_inv_2w_ = get_reparam_inv_2w ( + param_x0_,param_log_w_ + ); + + }; + + double get_reparam_inv_2w ( + double x0,double log_w + ) { + return 0.5 * exp(-log_w); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_x0_; + double param_log_w_; + + double reparam_inv_2w_; + +}; + + +/* +This is equivalent to a "scale mixture" of :class:`ExpSquaredKernel` +kernels with different scale lengths drawn from a gamma distribution. +See R&W for more info but here's the equation: + +.. math:: + k(r^2) = \left[1 - \frac{r^2}{2\,\alpha} \right]^\alpha + +:param log_alpha: + The Gamma distribution parameter. +*/ + +template +class RationalQuadraticKernel : public Kernel { +public: + RationalQuadraticKernel ( + double log_alpha, + int blocked, + const double* min_block, + const double* max_block, + size_t ndim, + size_t naxes + ) : + size_(1), + metric_(ndim, naxes), + blocked_(blocked), + min_block_(naxes), + max_block_(naxes) + , param_log_alpha_(log_alpha) + + { + size_t i; + if (blocked_) { + for (i = 0; i < naxes; ++i) { + min_block_[i] = min_block[i]; + max_block_[i] = max_block[i]; + } + } + update_reparams(); + }; + + size_t get_ndim () const { return metric_.get_ndim(); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_log_alpha_; + return metric_.get_parameter(i - size_); + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_log_alpha_ = value; + update_reparams(); + } else + metric_.set_parameter(i - size_, value); + }; + + double get_metric_parameter (size_t i) const { + return metric_.get_parameter(i); + }; + void set_metric_parameter (size_t i, double value) { + metric_.set_parameter(i, value); + }; + + size_t get_axis (size_t i) const { + return metric_.get_axis(i); + }; + void set_axis (size_t i, size_t value) { + metric_.set_axis(i, value); + }; + + double get_value ( + double log_alpha, + double alpha, + + double r2) { + return pow(1 + 0.5 * r2 / alpha, -alpha); + + }; + + double value (const double* x1, const double* x2) { + if (blocked_) { + size_t i, j; + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) + return 0.0; + } + } + double r2 = metric_.value(x1, x2); + return get_value( + param_log_alpha_, + reparam_alpha_, + + r2); + }; + + double log_alpha_gradient ( + double log_alpha, + double alpha, + + double r2) { + double t1 = 1.0 + 0.5 * r2 / alpha, + t2 = 2.0 * alpha * t1; + return alpha * pow(t1, -alpha) * (r2 / t2 - log(t1)); + + }; + double radial_gradient ( + double log_alpha, + double alpha, + + double r2) { + return -0.5 * pow(1 + 0.5 * r2 / alpha, -alpha-1); + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + bool out = false; + size_t i, j, n = size(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + if (which[0]) + grad[0] = log_alpha_gradient( + param_log_alpha_, + reparam_alpha_, + + r2); + + + bool any = false; + for (i = size_; i < size(); ++i) if (which[i]) { any = true; break; } + if (any) { + double r2grad = radial_gradient( + param_log_alpha_, + + reparam_alpha_, + + r2); + metric_.gradient(x1, x2, &(grad[size_])); + for (i = size_; i < n; ++i) grad[i] *= r2grad; + } + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + param_log_alpha_, + + reparam_alpha_, + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= r2grad; + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + param_log_alpha_, + + reparam_alpha_, + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= -r2grad; + }; + + size_t size () const { return metric_.size() + size_; }; + + void update_reparams () { + reparam_alpha_ = get_reparam_alpha ( + param_log_alpha_ + ); + + }; + + double get_reparam_alpha ( + double log_alpha + ) { + return exp(log_alpha); + } + + +private: + size_t size_; + M metric_; + int blocked_; + std::vector min_block_, max_block_; + double param_log_alpha_; + + double reparam_alpha_; + +}; + + +/* +The exponential kernel is a stationary kernel where the value +at a given radius :math:`r^2` is given by: + +.. math:: + + k(r^2) = \exp \left ( -\sqrt{r^2} \right ) +*/ + +template +class ExpKernel : public Kernel { +public: + ExpKernel ( + int blocked, + const double* min_block, + const double* max_block, + size_t ndim, + size_t naxes + ) : + size_(0), + metric_(ndim, naxes), + blocked_(blocked), + min_block_(naxes), + max_block_(naxes) + + { + size_t i; + if (blocked_) { + for (i = 0; i < naxes; ++i) { + min_block_[i] = min_block[i]; + max_block_[i] = max_block[i]; + } + } + update_reparams(); + }; + + size_t get_ndim () const { return metric_.get_ndim(); }; + + double get_parameter (size_t i) const { + return metric_.get_parameter(i - size_); + }; + void set_parameter (size_t i, double value) { + metric_.set_parameter(i - size_, value); + }; + + double get_metric_parameter (size_t i) const { + return metric_.get_parameter(i); + }; + void set_metric_parameter (size_t i, double value) { + metric_.set_parameter(i, value); + }; + + size_t get_axis (size_t i) const { + return metric_.get_axis(i); + }; + void set_axis (size_t i, size_t value) { + metric_.set_axis(i, value); + }; + + double get_value ( + + double r2) { + return exp(-sqrt(r2)); + }; + + double value (const double* x1, const double* x2) { + if (blocked_) { + size_t i, j; + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) + return 0.0; + } + } + double r2 = metric_.value(x1, x2); + return get_value( + + r2); + }; + + double radial_gradient ( + + double r2) { + if (r2 < DBL_EPSILON) return 0.0; + double r = sqrt(r2); + return -0.5 * exp(-r) / r; + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + bool out = false; + size_t i, j, n = size(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + + + bool any = false; + for (i = size_; i < size(); ++i) if (which[i]) { any = true; break; } + if (any) { + double r2grad = radial_gradient( + + + r2); + metric_.gradient(x1, x2, &(grad[size_])); + for (i = size_; i < n; ++i) grad[i] *= r2grad; + } + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= r2grad; + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= -r2grad; + }; + + size_t size () const { return metric_.size() + size_; }; + + void update_reparams () { + + }; + + + +private: + size_t size_; + M metric_; + int blocked_; + std::vector min_block_, max_block_; + + +}; + + +/* +A local Gaussian kernel. + +.. math:: + k(\mathbf{x}_i,\,\mathbf{x}_j) = \exp\left( + -\frac{(x_i - x_0)^2 + (x_j - x_0)^2}{2\,w} \right)) + +:param location: + The location :math:`x_0` of the Gaussian. + +:param log_width: + The (squared) width :math:`w` of the Gaussian. +*/ + +class LocalGaussianKernel : public Kernel { +public: + LocalGaussianKernel ( + double location, + double log_width, + size_t ndim, + size_t naxes + ) : + size_(2), + subspace_(ndim, naxes) + , param_location_(location) + , param_log_width_(log_width) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_location_; + if (i == 1) return param_log_width_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_location_ = value; + update_reparams(); + } else + if (i == 1) { + param_log_width_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double location, + double log_width, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - location, d2 = x2 - location; + return exp(-(d1*d1 + d2*d2) * inv_2w); + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_location_, + param_log_width_, + reparam_inv_2w_, + + x1[j], x2[j]); + } + return value; + }; + + double location_gradient ( + double location, + double log_width, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - location, d2 = x2 - location; + return 2 * exp(-(d1*d1 + d2*d2) * inv_2w) * inv_2w * (d1 + d2); + + }; + double log_width_gradient ( + double location, + double log_width, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - location, d2 = x2 - location, + arg = (d1*d1 + d2*d2) * inv_2w; + return exp(-arg) * arg; + + }; + double _x1_gradient ( + double location, + double log_width, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - location, d2 = x2 - location; + return -2.0 * exp(-(d1*d1 + d2*d2) * inv_2w) * d1 * inv_2w; + + }; + + double _x2_gradient ( + double location, + double log_width, + double inv_2w, + + double x1, double x2) { + double d1 = x1 - location, d2 = x2 - location; + return -2.0 * exp(-(d1*d1 + d2*d2) * inv_2w) * d2 * inv_2w; + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + grad[1] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += location_gradient( + param_location_, + param_log_width_, + reparam_inv_2w_, + + x1[j], x2[j]); + if (which[1]) + grad[1] += log_width_gradient( + param_location_, + param_log_width_, + reparam_inv_2w_, + + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_location_, + param_log_width_, + reparam_inv_2w_, + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_location_, + param_log_width_, + reparam_inv_2w_, + + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_inv_2w_ = get_reparam_inv_2w ( + param_location_,param_log_width_ + ); + + }; + + double get_reparam_inv_2w ( + double location,double log_width + ) { + return 0.5 * exp(-log_width); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_location_; + double param_log_width_; + + double reparam_inv_2w_; + +}; + + +/* +This kernel is a no-op*/ + +class EmptyKernel : public Kernel { +public: + EmptyKernel ( + size_t ndim, + size_t naxes + ) : + size_(0), + subspace_(ndim, naxes) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + return 0.0; + }; + void set_parameter (size_t i, double value) { + ; + }; + + double get_value ( + + double x1, double x2) { + return 0.0; + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + + x1[j], x2[j]); + } + return value; + }; + + double _x1_gradient ( + + double x1, double x2) { + return 0.0; + }; + + double _x2_gradient ( + + double x1, double x2) { + return 0.0; + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + + + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + + x1[j], x2[j]); + } + }; + + void update_reparams () { + + }; + + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + + +}; + + +/* +The simplest periodic kernel. This + +.. math:: + k(\mathbf{x}_i,\,\mathbf{x}_j) = \cos\left( + \frac{2\,\pi\,|x_i - x_j|}{P} \right) + +where the parameter :math:`P` is the period of the oscillation. This +kernel should probably always be multiplied be a stationary kernel +(e.g. :class:`ExpSquaredKernel`) to allow quasi-periodic variations. + +:param log_period: + The period of the oscillation. +*/ + +class CosineKernel : public Kernel { +public: + CosineKernel ( + double log_period, + size_t ndim, + size_t naxes + ) : + size_(1), + subspace_(ndim, naxes) + , param_log_period_(log_period) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_log_period_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_log_period_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double log_period, + double factor, + + double x1, double x2) { + return cos((x1 - x2) * factor); + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + } + return value; + }; + + double log_period_gradient ( + double log_period, + double factor, + + double x1, double x2) { + double r = factor * (x1 - x2); + return r * sin(r); + + }; + double _x1_gradient ( + double log_period, + double factor, + + double x1, double x2) { + return -factor*sin(factor * (x1-x2)); + + }; + + double _x2_gradient ( + double log_period, + double factor, + + double x1, double x2) { + return factor*sin(factor * (x1-x2)); + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += log_period_gradient( + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_factor_ = get_reparam_factor ( + param_log_period_ + ); + + }; + + double get_reparam_factor ( + double log_period + ) { + return 2 * M_PI * exp(-log_period); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_log_period_; + + double reparam_factor_; + +}; + + +/* +The Matern-5/2 kernel is stationary kernel where the value at a +given radius :math:`r^2` is given by: + +.. math:: + + k(r^2) = \left( 1+\sqrt{5\,r^2}+ \frac{5\,r^2}{3} \right)\, + \exp \left (-\sqrt{5\,r^2} \right ) +*/ + +template +class Matern52Kernel : public Kernel { +public: + Matern52Kernel ( + int blocked, + const double* min_block, + const double* max_block, + size_t ndim, + size_t naxes + ) : + size_(0), + metric_(ndim, naxes), + blocked_(blocked), + min_block_(naxes), + max_block_(naxes) + + { + size_t i; + if (blocked_) { + for (i = 0; i < naxes; ++i) { + min_block_[i] = min_block[i]; + max_block_[i] = max_block[i]; + } + } + update_reparams(); + }; + + size_t get_ndim () const { return metric_.get_ndim(); }; + + double get_parameter (size_t i) const { + return metric_.get_parameter(i - size_); + }; + void set_parameter (size_t i, double value) { + metric_.set_parameter(i - size_, value); + }; + + double get_metric_parameter (size_t i) const { + return metric_.get_parameter(i); + }; + void set_metric_parameter (size_t i, double value) { + metric_.set_parameter(i, value); + }; + + size_t get_axis (size_t i) const { + return metric_.get_axis(i); + }; + void set_axis (size_t i, size_t value) { + metric_.set_axis(i, value); + }; + + double get_value ( + + double r2) { + double r = sqrt(5.0 * r2); + return (1 + r + 5.0 * r2 / 3.0) * exp(-r); + + }; + + double value (const double* x1, const double* x2) { + if (blocked_) { + size_t i, j; + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) + return 0.0; + } + } + double r2 = metric_.value(x1, x2); + return get_value( + + r2); + }; + + double radial_gradient ( + + double r2) { + double r = sqrt(5.0 * r2); + return -5 * (1 + r) * exp(-r) / 6.0; + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + bool out = false; + size_t i, j, n = size(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + + + bool any = false; + for (i = size_; i < size(); ++i) if (which[i]) { any = true; break; } + if (any) { + double r2grad = radial_gradient( + + + r2); + metric_.gradient(x1, x2, &(grad[size_])); + for (i = size_; i < n; ++i) grad[i] *= r2grad; + } + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= r2grad; + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= -r2grad; + }; + + size_t size () const { return metric_.size() + size_; }; + + void update_reparams () { + + }; + + + +private: + size_t size_; + M metric_; + int blocked_; + std::vector min_block_, max_block_; + + +}; + + +/* +The exp-sine-squared kernel is a commonly used periodic kernel. Unlike +the :class:`CosineKernel`, this kernel never has negative covariance +which might be useful for your problem. Here's the equation: + +.. math:: + k(\mathbf{x}_i,\,\mathbf{x}_j) = + \exp \left( -\Gamma\,\sin^2\left[ + \frac{\pi}{P}\,\left|x_i-x_j\right| + \right] \right) + +:param gamma: + The scale :math:`\Gamma` of the correlations. + +:param log_period: + The log of the period :math:`P` of the oscillation (in the same units + as :math:`\mathbf{x}`). +*/ + +class ExpSine2Kernel : public Kernel { +public: + ExpSine2Kernel ( + double gamma, + double log_period, + size_t ndim, + size_t naxes + ) : + size_(2), + subspace_(ndim, naxes) + , param_gamma_(gamma) + , param_log_period_(log_period) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_gamma_; + if (i == 1) return param_log_period_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_gamma_ = value; + update_reparams(); + } else + if (i == 1) { + param_log_period_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double gamma, + double log_period, + double factor, + + double x1, double x2) { + double s = sin((x1 - x2) * factor); + return exp(-gamma * s * s); + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_gamma_, + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + } + return value; + }; + + double gamma_gradient ( + double gamma, + double log_period, + double factor, + + double x1, double x2) { + double s = sin((x1 - x2) * factor), s2 = s * s; + return -s2 * exp(-gamma * s2); + + }; + double log_period_gradient ( + double gamma, + double log_period, + double factor, + + double x1, double x2) { + double arg = (x1 - x2) * factor, + s = sin(arg), c = cos(arg), + A = exp(-gamma * s * s); + return 2 * gamma * arg * c * s * A; + + }; + double _x1_gradient ( + double gamma, + double log_period, + double factor, + + double x1, double x2) { + double d = x1 - x2; + double s = sin(d * factor); + return -exp(-gamma * s * s) * factor * gamma * sin(2.0 * factor * d); + + }; + + double _x2_gradient ( + double gamma, + double log_period, + double factor, + + double x1, double x2) { + double d = x1 - x2; + double s = sin(d * factor); + return exp(-gamma * s * s) * factor * gamma * sin(2.0 * factor * d); + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + grad[1] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += gamma_gradient( + param_gamma_, + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + if (which[1]) + grad[1] += log_period_gradient( + param_gamma_, + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_gamma_, + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_gamma_, + param_log_period_, + reparam_factor_, + + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_factor_ = get_reparam_factor ( + param_gamma_,param_log_period_ + ); + + }; + + double get_reparam_factor ( + double gamma,double log_period + ) { + return M_PI * exp(-log_period); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_gamma_; + double param_log_period_; + + double reparam_factor_; + +}; + + +/* +This kernel returns the constant + +.. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = c + +where :math:`c` is a parameter. + +:param log_constant: + The log of :math:`c` in the above equation. +*/ + +class ConstantKernel : public Kernel { +public: + ConstantKernel ( + double log_constant, + size_t ndim, + size_t naxes + ) : + size_(1), + subspace_(ndim, naxes) + , param_log_constant_(log_constant) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_log_constant_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_log_constant_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double log_constant, + double constant, + + double x1, double x2) { + return constant; + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_log_constant_, + reparam_constant_, + + x1[j], x2[j]); + } + return value; + }; + + double log_constant_gradient ( + double log_constant, + double constant, + + double x1, double x2) { + return constant; + + }; + double _x1_gradient ( + double log_constant, + double constant, + + double x1, double x2) { + return 0.0; + + }; + + double _x2_gradient ( + double log_constant, + double constant, + + double x1, double x2) { + return 0.0; + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += log_constant_gradient( + param_log_constant_, + reparam_constant_, + + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_log_constant_, + reparam_constant_, + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_log_constant_, + reparam_constant_, + + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_constant_ = get_reparam_constant ( + param_log_constant_ + ); + + }; + + double get_reparam_constant ( + double log_constant + ) { + return exp(log_constant); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_log_constant_; + + double reparam_constant_; + +}; + + +/* +The exponential-squared kernel is a stationary kernel where the value +at a given radius :math:`r^2` is given by: + +.. math:: + + k(r^2) = \exp \left ( -\frac{r^2}{2} \right ) +*/ + +template +class ExpSquaredKernel : public Kernel { +public: + ExpSquaredKernel ( + int blocked, + const double* min_block, + const double* max_block, + size_t ndim, + size_t naxes + ) : + size_(0), + metric_(ndim, naxes), + blocked_(blocked), + min_block_(naxes), + max_block_(naxes) + + { + size_t i; + if (blocked_) { + for (i = 0; i < naxes; ++i) { + min_block_[i] = min_block[i]; + max_block_[i] = max_block[i]; + } + } + update_reparams(); + }; + + size_t get_ndim () const { return metric_.get_ndim(); }; + + double get_parameter (size_t i) const { + return metric_.get_parameter(i - size_); + }; + void set_parameter (size_t i, double value) { + metric_.set_parameter(i - size_, value); + }; + + double get_metric_parameter (size_t i) const { + return metric_.get_parameter(i); + }; + void set_metric_parameter (size_t i, double value) { + metric_.set_parameter(i, value); + }; + + size_t get_axis (size_t i) const { + return metric_.get_axis(i); + }; + void set_axis (size_t i, size_t value) { + metric_.set_axis(i, value); + }; + + double get_value ( + + double r2) { + return exp(-0.5 * r2); + }; + + double value (const double* x1, const double* x2) { + if (blocked_) { + size_t i, j; + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) + return 0.0; + } + } + double r2 = metric_.value(x1, x2); + return get_value( + + r2); + }; + + double radial_gradient ( + + double r2) { + return -0.5 * exp(-0.5 * r2); + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + bool out = false; + size_t i, j, n = size(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + + + bool any = false; + for (i = size_; i < size(); ++i) if (which[i]) { any = true; break; } + if (any) { + double r2grad = radial_gradient( + + + r2); + metric_.gradient(x1, x2, &(grad[size_])); + for (i = size_; i < n; ++i) grad[i] *= r2grad; + } + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= r2grad; + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= -r2grad; + }; + + size_t size () const { return metric_.size() + size_; }; + + void update_reparams () { + + }; + + + +private: + size_t size_; + M metric_; + int blocked_; + std::vector min_block_, max_block_; + + +}; + + +/* +The Matern-3/2 kernel is stationary kernel where the value at a +given radius :math:`r^2` is given by: + +.. math:: + + k(r^2) = \left( 1+\sqrt{3\,r^2} \right)\, + \exp \left (-\sqrt{3\,r^2} \right ) +*/ + +template +class Matern32Kernel : public Kernel { +public: + Matern32Kernel ( + int blocked, + const double* min_block, + const double* max_block, + size_t ndim, + size_t naxes + ) : + size_(0), + metric_(ndim, naxes), + blocked_(blocked), + min_block_(naxes), + max_block_(naxes) + + { + size_t i; + if (blocked_) { + for (i = 0; i < naxes; ++i) { + min_block_[i] = min_block[i]; + max_block_[i] = max_block[i]; + } + } + update_reparams(); + }; + + size_t get_ndim () const { return metric_.get_ndim(); }; + + double get_parameter (size_t i) const { + return metric_.get_parameter(i - size_); + }; + void set_parameter (size_t i, double value) { + metric_.set_parameter(i - size_, value); + }; + + double get_metric_parameter (size_t i) const { + return metric_.get_parameter(i); + }; + void set_metric_parameter (size_t i, double value) { + metric_.set_parameter(i, value); + }; + + size_t get_axis (size_t i) const { + return metric_.get_axis(i); + }; + void set_axis (size_t i, size_t value) { + metric_.set_axis(i, value); + }; + + double get_value ( + + double r2) { + double r = sqrt(3.0 * r2); + return (1.0 + r) * exp(-r); + + }; + + double value (const double* x1, const double* x2) { + if (blocked_) { + size_t i, j; + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) + return 0.0; + } + } + double r2 = metric_.value(x1, x2); + return get_value( + + r2); + }; + + double radial_gradient ( + + double r2) { + double r = sqrt(3.0 * r2); + return -3.0 * 0.5 * exp(-r); + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + bool out = false; + size_t i, j, n = size(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + + + bool any = false; + for (i = size_; i < size(); ++i) if (which[i]) { any = true; break; } + if (any) { + double r2grad = radial_gradient( + + + r2); + metric_.gradient(x1, x2, &(grad[size_])); + for (i = size_; i < n; ++i) grad[i] *= r2grad; + } + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= r2grad; + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + bool out = false; + size_t i, j, n = this->get_ndim(); + if (blocked_) { + for (i = 0; i < min_block_.size(); ++i) { + j = metric_.get_axis(i); + if (x1[j] < min_block_[i] || x1[j] > max_block_[i] || + x2[j] < min_block_[i] || x2[j] > max_block_[i]) { + out = true; + break; + } + } + if (out) { + for (i = 0; i < n; ++i) grad[i] = 0.0; + return; + } + } + + double r2 = metric_.value(x1, x2); + double r2grad = 2.0 * radial_gradient( + + + r2); + metric_.x1_gradient(x1, x2, grad); + for (i = 0; i < n; ++i) grad[i] *= -r2grad; + }; + + size_t size () const { return metric_.size() + size_; }; + + void update_reparams () { + + }; + + + +private: + size_t size_; + M metric_; + int blocked_; + std::vector min_block_, max_block_; + + +}; + + +/* +A polynomial kernel + +.. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = + (\mathbf{x}_i \cdot \mathbf{x}_j + \sigma^2)^P + +:param order: + The power :math:`P`. This parameter is a *constant*; it is not + included in the parameter vector. + +:param log_sigma2: + The variance :math:`\sigma^2 > 0`. +*/ + +class PolynomialKernel : public Kernel { +public: + PolynomialKernel ( + double log_sigma2, + double order, + size_t ndim, + size_t naxes + ) : + size_(1), + subspace_(ndim, naxes) + , param_log_sigma2_(log_sigma2) + , constant_order_(order) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + if (i == 0) return param_log_sigma2_; + return 0.0; + }; + void set_parameter (size_t i, double value) { + if (i == 0) { + param_log_sigma2_ = value; + update_reparams(); + } else + ; + }; + + double get_value ( + double log_sigma2, + double sigma2, + + double order, + double x1, double x2) { + if (order == 0.0) return 1.0; + return pow(x1 * x2 + sigma2, order); + + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + param_log_sigma2_, + reparam_sigma2_, + + constant_order_, + x1[j], x2[j]); + } + return value; + }; + + double log_sigma2_gradient ( + double log_sigma2, + double sigma2, + + double order, + double x1, double x2) { + if (order == 0.0) return 0.0; + return sigma2 * pow(x1 * x2 + sigma2, order-1.0) * order; + + }; + double _x1_gradient ( + double log_sigma2, + double sigma2, + + double order, + double x1, double x2) { + if (order == 0.0) return 0.0; + return x2 * order * pow(x1 * x2 + sigma2, order-1.0); + + }; + + double _x2_gradient ( + double log_sigma2, + double sigma2, + + double order, + double x1, double x2) { + if (order == 0.0) return 0.0; + return x1 * order * pow(x1 * x2 + sigma2, order-1.0); + + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + grad[0] = 0.0; + + + size_t i, j, n = subspace_.get_naxes(); + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + if (which[0]) + grad[0] += log_sigma2_gradient( + param_log_sigma2_, + reparam_sigma2_, + + constant_order_, + x1[j], x2[j]); + + } + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + param_log_sigma2_, + reparam_sigma2_, + + constant_order_, + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + param_log_sigma2_, + reparam_sigma2_, + + constant_order_, + x1[j], x2[j]); + } + }; + + void update_reparams () { + reparam_sigma2_ = get_reparam_sigma2 ( + param_log_sigma2_, + constant_order_ + ); + + }; + + double get_reparam_sigma2 ( + double log_sigma2, + double order + ) { + return exp(log_sigma2); + } + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + double param_log_sigma2_; + + double reparam_sigma2_; + + double constant_order_; +}; + + +/* +The dot product kernel + +.. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = \mathbf{x}_i \cdot \mathbf{x}_j + +with no parameters. +*/ + +class DotProductKernel : public Kernel { +public: + DotProductKernel ( + size_t ndim, + size_t naxes + ) : + size_(0), + subspace_(ndim, naxes) + { + update_reparams(); + }; + + size_t get_ndim () const { return subspace_.get_ndim(); }; + size_t get_axis (size_t i) const { return subspace_.get_axis(i); }; + void set_axis (size_t i, size_t value) { subspace_.set_axis(i, value); }; + + double get_parameter (size_t i) const { + return 0.0; + }; + void set_parameter (size_t i, double value) { + ; + }; + + double get_value ( + + double x1, double x2) { + return x1 * x2; + }; + + double value (const double* x1, const double* x2) { + size_t i, j, n = subspace_.get_naxes(); + double value = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + value += get_value( + + x1[j], x2[j]); + } + return value; + }; + + double _x1_gradient ( + + double x1, double x2) { + return x2; + }; + + double _x2_gradient ( + + double x1, double x2) { + return x1; + }; + + void gradient (const double* x1, const double* x2, const unsigned* which, + double* grad) { + + + + }; + + void x1_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x1_gradient( + + x1[j], x2[j]); + } + }; + + void x2_gradient (const double* x1, const double* x2, double* grad) { + size_t i, j, ndim = this->get_ndim(), n = subspace_.get_naxes(); + for (i = 0; i < ndim; ++i) grad[i] = 0.0; + for (i = 0; i < n; ++i) { + j = subspace_.get_axis(i); + grad[j] = _x2_gradient( + + x1[j], x2[j]); + } + }; + + void update_reparams () { + + }; + + + + size_t size () const { return size_; }; + +private: + size_t size_; + george::subspace::Subspace subspace_; + + +}; + +}; // namespace kernels +}; // namespace george + +#endif \ No newline at end of file diff --git a/src/george/kernels.py b/src/george/kernels.py new file mode 100644 index 0000000..c6dbd1c --- /dev/null +++ b/src/george/kernels.py @@ -0,0 +1,1003 @@ +# -*- coding: utf-8 -*- + +from __future__ import division, print_function + +__all__ = [ + "Kernel", "Sum", "Product", + "LinearKernel", + "MyLocalGaussianKernel", + "RationalQuadraticKernel", + "ExpKernel", + "LocalGaussianKernel", + "EmptyKernel", + "CosineKernel", + "Matern52Kernel", + "ExpSine2Kernel", + "ConstantKernel", + "ExpSquaredKernel", + "Matern32Kernel", + "PolynomialKernel", + "DotProductKernel", +] + +import numpy as np + +from .modeling import Model, ModelSet +from .metrics import Metric, Subspace +from .kernel_interface import KernelInterface + + +class Kernel(ModelSet): + """ + The abstract kernel type. Every kernel implemented in George should be + a subclass of this object. + + """ + + is_kernel = True + kernel_type = -1 + + # This function deals with weird behavior when performing arithmetic + # operations with numpy scalars. + def __array_wrap__(self, array, context=None): + if context is None: + raise TypeError("Invalid operation") + ufunc, args, _ = context + if ufunc.__name__ == "multiply": + return float(args[0]) * args[1] + elif ufunc.__name__ == "add": + return float(args[0]) + args[1] + raise TypeError("Invalid operation") + __array_priority__ = np.inf + + def __getstate__(self): + odict = self.__dict__.copy() + odict["_kernel"] = None + return odict + + # We must overload the ModelSet attribute getter to pass the requests + # to the "BaseKernel" + def __getattr__(self, name): + if "models" in self.__dict__: + if name in self.models: + return self.models[name] + if None in self.models: + return getattr(self.models[None], name) + raise AttributeError(name) + + @property + def kernel(self): + return KernelInterface(self) + + def __repr__(self): + kernel = self.models[None] + params = ["{0}={1}".format(k, getattr(kernel, k)) + for k in kernel.parameter_names] + if self.stationary: + params += ["metric={0}".format(repr(self.metric)), + "block={0}".format(repr(self.block))] + else: + params += ["ndim={0}".format(self.ndim), + "axes={0}".format(repr(self.axes))] + return "{0}({1})".format(self.__class__.__name__, ", ".join(params)) + + def __add__(self, b): + if not hasattr(b, "is_kernel"): + return Sum(ConstantKernel(log_constant=np.log(float(b)/self.ndim), + ndim=self.ndim), self) + return Sum(self, b) + + def __radd__(self, b): + return self.__add__(b) + + def __mul__(self, b): + if not hasattr(b, "is_kernel"): + log_constant = np.log(float(b)/self.ndim) + return Product(ConstantKernel(log_constant=log_constant, + ndim=self.ndim), self) + return Product(self, b) + + def __rmul__(self, b): + return self.__mul__(b) + + def get_value(self, x1, x2=None, diag=False): + x1 = np.ascontiguousarray(x1, dtype=np.float64) + if x2 is None: + if diag: + return self.kernel.value_diagonal(x1, x1) + else: + return self.kernel.value_symmetric(x1) + x2 = np.ascontiguousarray(x2, dtype=np.float64) + if diag: + return self.kernel.value_diagonal(x1, x2) + else: + return self.kernel.value_general(x1, x2) + + def get_gradient(self, x1, x2=None, include_frozen=False): + mask = ( + np.ones(self.full_size, dtype=bool) + if include_frozen else self.unfrozen_mask + ) + which = mask.astype(np.uint32) + x1 = np.ascontiguousarray(x1, dtype=np.float64) + if x2 is None: + g = self.kernel.gradient_symmetric(which, x1) + else: + x2 = np.ascontiguousarray(x2, dtype=np.float64) + g = self.kernel.gradient_general(which, x1, x2) + return g[:, :, mask] + + def get_x1_gradient(self, x1, x2=None): + x1 = np.ascontiguousarray(x1, dtype=np.float64) + if x2 is None: + x2 = x1 + else: + x2 = np.ascontiguousarray(x2, dtype=np.float64) + return self.kernel.x1_gradient_general(x1, x2) + + def get_x2_gradient(self, x1, x2=None): + x1 = np.ascontiguousarray(x1, dtype=np.float64) + if x2 is None: + x2 = x1 + else: + x2 = np.ascontiguousarray(x2, dtype=np.float64) + return self.kernel.x2_gradient_general(x1, x2) + + def test_gradient(self, x1, x2=None, eps=1.32e-6, **kwargs): + vector = self.get_parameter_vector() + g0 = self.get_gradient(x1, x2=x2) + + for i, v in enumerate(vector): + vector[i] = v + eps + self.set_parameter_vector(vector) + kp = self.get_value(x1, x2=x2) + + vector[i] = v - eps + self.set_parameter_vector(vector) + km = self.get_value(x1, x2=x2) + + vector[i] = v + self.set_parameter_vector(vector) + + grad = 0.5 * (kp - km) / eps + assert np.allclose(g0[:, :, i], grad, **kwargs), \ + "incorrect gradient for parameter '{0}' ({1})" \ + .format(self.get_parameter_names()[i], i) + + def test_x1_gradient(self, x1, x2=None, eps=1.32e-6, **kwargs): + kwargs["atol"] = kwargs.get("atol", 0.5 * eps) + g0 = self.get_x1_gradient(x1, x2=x2) + if x2 is None: + x2 = np.array(x1) + for i in range(len(x1)): + for k in range(self.ndim): + x1[i, k] += eps + kp = self.get_value(x1, x2=x2) + + x1[i, k] -= 2*eps + km = self.get_value(x1, x2=x2) + + x1[i, k] += eps + + grad = 0.5 * (kp - km) / eps + assert np.allclose(g0[i, :, k], grad[i], **kwargs) + + def test_x2_gradient(self, x1, x2=None, eps=1.32e-6, **kwargs): + kwargs["atol"] = kwargs.get("atol", 0.5 * eps) + g0 = self.get_x2_gradient(x1, x2=x2) + if x2 is None: + x2 = np.array(x1) + for i in range(len(x2)): + for k in range(self.ndim): + x2[i, k] += eps + kp = self.get_value(x1, x2=x2) + + x2[i, k] -= 2*eps + km = self.get_value(x1, x2=x2) + + x2[i, k] += eps + + grad = 0.5 * (kp - km) / eps + assert np.allclose(g0[:, i, k], grad[:, i], **kwargs) + + +class _operator(Kernel): + is_kernel = False + kernel_type = -1 + operator_type = -1 + + def __init__(self, k1, k2): + if k1.ndim != k2.ndim: + raise ValueError("Dimension mismatch") + self.ndim = k1.ndim + self._dirty = True + super(_operator, self).__init__([("k1", k1), ("k2", k2)]) + + @property + def k1(self): + return self.models["k1"] + + @property + def k2(self): + return self.models["k2"] + + @property + def dirty(self): + return self._dirty or self.k1.dirty or self.k2.dirty + + @dirty.setter + def dirty(self, v): + self._dirty = v + self.k1.dirty = False + self.k2.dirty = False + + +class Sum(_operator): + is_kernel = False + operator_type = 0 + + def __repr__(self): + return "{0} + {1}".format(self.k1, self.k2) + + +class Product(_operator): + is_kernel = False + operator_type = 1 + + def __repr__(self): + return "{0} * {1}".format(self.k1, self.k2) + + +class BaseLinearKernel (Model): + parameter_names = ("log_gamma2", ) + + +class LinearKernel (Kernel): + r""" + The linear regression kernel + + .. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = + \frac{(\mathbf{x}_i \cdot \mathbf{x}_j)^P}{\gamma^2} + + :param order: + The power :math:`P`. This parameter is a *constant*; it is not + included in the parameter vector. + + :param log_gamma2: + The scale factor :math:`\gamma^2`. + + + """ + + kernel_type = 0 + stationary = False + + def __init__(self, + log_gamma2=None, + order=None, + bounds=None, + ndim=1, + axes=None): + + if order is None: + raise ValueError("missing required parameter 'order'") + self.order = order + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(log_gamma2=log_gamma2, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseLinearKernel(**kwargs) + super(LinearKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseMyLocalGaussianKernel (Model): + parameter_names = ("x0", "log_w", ) + + +class MyLocalGaussianKernel (Kernel): + r""" + You should always document your code. + + """ + + kernel_type = 1 + stationary = False + + def __init__(self, + x0=None, + log_w=None, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(x0=x0, log_w=log_w, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseMyLocalGaussianKernel(**kwargs) + super(MyLocalGaussianKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseRationalQuadraticKernel (Model): + parameter_names = ("log_alpha", ) + + +class RationalQuadraticKernel (Kernel): + r""" + This is equivalent to a "scale mixture" of :class:`ExpSquaredKernel` + kernels with different scale lengths drawn from a gamma distribution. + See R&W for more info but here's the equation: + + .. math:: + k(r^2) = \left[1 - \frac{r^2}{2\,\alpha} \right]^\alpha + + :param log_alpha: + The Gamma distribution parameter. + + + """ + + kernel_type = 2 + stationary = True + + def __init__(self, + log_alpha=None, + metric=None, + metric_bounds=None, + lower=True, + block=None, + bounds=None, + ndim=1, + axes=None): + + if metric is None: + raise ValueError("missing required parameter 'metric'") + metric = Metric(metric, bounds=metric_bounds, ndim=ndim, + axes=axes, lower=lower) + self.ndim = metric.ndim + self.axes = metric.axes + self.block = block + + kwargs = dict(log_alpha=log_alpha, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseRationalQuadraticKernel(**kwargs) + super(RationalQuadraticKernel, self).__init__([ + (None, base), ("metric", metric) + ]) + + # Common setup. + self.dirty = True + + @property + def block(self): + if not self.blocked: + return None + return list(zip(self.min_block, self.max_block)) + + @block.setter + def block(self, block): + if block is None: + self.blocked = False + self.min_block = -np.inf + np.zeros(len(self.axes)) + self.max_block = np.inf + np.zeros(len(self.axes)) + return + + block = np.atleast_2d(block) + if block.shape != (len(self.axes), 2): + raise ValueError("dimension mismatch in block specification") + self.blocked = True + self.min_block, self.max_block = map(np.array, zip(*block)) + + +class BaseExpKernel (Model): + parameter_names = () + + +class ExpKernel (Kernel): + r""" + The exponential kernel is a stationary kernel where the value + at a given radius :math:`r^2` is given by: + + .. math:: + + k(r^2) = \exp \left ( -\sqrt{r^2} \right ) + + + """ + + kernel_type = 3 + stationary = True + + def __init__(self, + metric=None, + metric_bounds=None, + lower=True, + block=None, + bounds=None, + ndim=1, + axes=None): + + if metric is None: + raise ValueError("missing required parameter 'metric'") + metric = Metric(metric, bounds=metric_bounds, ndim=ndim, + axes=axes, lower=lower) + self.ndim = metric.ndim + self.axes = metric.axes + self.block = block + + kwargs = dict() + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseExpKernel(**kwargs) + super(ExpKernel, self).__init__([ + (None, base), ("metric", metric) + ]) + + # Common setup. + self.dirty = True + + @property + def block(self): + if not self.blocked: + return None + return list(zip(self.min_block, self.max_block)) + + @block.setter + def block(self, block): + if block is None: + self.blocked = False + self.min_block = -np.inf + np.zeros(len(self.axes)) + self.max_block = np.inf + np.zeros(len(self.axes)) + return + + block = np.atleast_2d(block) + if block.shape != (len(self.axes), 2): + raise ValueError("dimension mismatch in block specification") + self.blocked = True + self.min_block, self.max_block = map(np.array, zip(*block)) + + +class BaseLocalGaussianKernel (Model): + parameter_names = ("location", "log_width", ) + + +class LocalGaussianKernel (Kernel): + r""" + A local Gaussian kernel. + + .. math:: + k(\mathbf{x}_i,\,\mathbf{x}_j) = \exp\left( + -\frac{(x_i - x_0)^2 + (x_j - x_0)^2}{2\,w} \right)) + + :param location: + The location :math:`x_0` of the Gaussian. + + :param log_width: + The (squared) width :math:`w` of the Gaussian. + + + """ + + kernel_type = 4 + stationary = False + + def __init__(self, + location=None, + log_width=None, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(location=location, log_width=log_width, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseLocalGaussianKernel(**kwargs) + super(LocalGaussianKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseEmptyKernel (Model): + parameter_names = () + + +class EmptyKernel (Kernel): + r""" + This kernel is a no-op + + """ + + kernel_type = 5 + stationary = False + + def __init__(self, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict() + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseEmptyKernel(**kwargs) + super(EmptyKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseCosineKernel (Model): + parameter_names = ("log_period", ) + + +class CosineKernel (Kernel): + r""" + The simplest periodic kernel. This + + .. math:: + k(\mathbf{x}_i,\,\mathbf{x}_j) = \cos\left( + \frac{2\,\pi\,|x_i - x_j|}{P} \right) + + where the parameter :math:`P` is the period of the oscillation. This + kernel should probably always be multiplied be a stationary kernel + (e.g. :class:`ExpSquaredKernel`) to allow quasi-periodic variations. + + :param log_period: + The period of the oscillation. + + + """ + + kernel_type = 6 + stationary = False + + def __init__(self, + log_period=None, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(log_period=log_period, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseCosineKernel(**kwargs) + super(CosineKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseMatern52Kernel (Model): + parameter_names = () + + +class Matern52Kernel (Kernel): + r""" + The Matern-5/2 kernel is stationary kernel where the value at a + given radius :math:`r^2` is given by: + + .. math:: + + k(r^2) = \left( 1+\sqrt{5\,r^2}+ \frac{5\,r^2}{3} \right)\, + \exp \left (-\sqrt{5\,r^2} \right ) + + + """ + + kernel_type = 7 + stationary = True + + def __init__(self, + metric=None, + metric_bounds=None, + lower=True, + block=None, + bounds=None, + ndim=1, + axes=None): + + if metric is None: + raise ValueError("missing required parameter 'metric'") + metric = Metric(metric, bounds=metric_bounds, ndim=ndim, + axes=axes, lower=lower) + self.ndim = metric.ndim + self.axes = metric.axes + self.block = block + + kwargs = dict() + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseMatern52Kernel(**kwargs) + super(Matern52Kernel, self).__init__([ + (None, base), ("metric", metric) + ]) + + # Common setup. + self.dirty = True + + @property + def block(self): + if not self.blocked: + return None + return list(zip(self.min_block, self.max_block)) + + @block.setter + def block(self, block): + if block is None: + self.blocked = False + self.min_block = -np.inf + np.zeros(len(self.axes)) + self.max_block = np.inf + np.zeros(len(self.axes)) + return + + block = np.atleast_2d(block) + if block.shape != (len(self.axes), 2): + raise ValueError("dimension mismatch in block specification") + self.blocked = True + self.min_block, self.max_block = map(np.array, zip(*block)) + + +class BaseExpSine2Kernel (Model): + parameter_names = ("gamma", "log_period", ) + + +class ExpSine2Kernel (Kernel): + r""" + The exp-sine-squared kernel is a commonly used periodic kernel. Unlike + the :class:`CosineKernel`, this kernel never has negative covariance + which might be useful for your problem. Here's the equation: + + .. math:: + k(\mathbf{x}_i,\,\mathbf{x}_j) = + \exp \left( -\Gamma\,\sin^2\left[ + \frac{\pi}{P}\,\left|x_i-x_j\right| + \right] \right) + + :param gamma: + The scale :math:`\Gamma` of the correlations. + + :param log_period: + The log of the period :math:`P` of the oscillation (in the same units + as :math:`\mathbf{x}`). + + + """ + + kernel_type = 8 + stationary = False + + def __init__(self, + gamma=None, + log_period=None, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(gamma=gamma, log_period=log_period, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseExpSine2Kernel(**kwargs) + super(ExpSine2Kernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseConstantKernel (Model): + parameter_names = ("log_constant", ) + + +class ConstantKernel (Kernel): + r""" + This kernel returns the constant + + .. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = c + + where :math:`c` is a parameter. + + :param log_constant: + The log of :math:`c` in the above equation. + + + """ + + kernel_type = 9 + stationary = False + + def __init__(self, + log_constant=None, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(log_constant=log_constant, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseConstantKernel(**kwargs) + super(ConstantKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseExpSquaredKernel (Model): + parameter_names = () + + +class ExpSquaredKernel (Kernel): + r""" + The exponential-squared kernel is a stationary kernel where the value + at a given radius :math:`r^2` is given by: + + .. math:: + + k(r^2) = \exp \left ( -\frac{r^2}{2} \right ) + + + """ + + kernel_type = 10 + stationary = True + + def __init__(self, + metric=None, + metric_bounds=None, + lower=True, + block=None, + bounds=None, + ndim=1, + axes=None): + + if metric is None: + raise ValueError("missing required parameter 'metric'") + metric = Metric(metric, bounds=metric_bounds, ndim=ndim, + axes=axes, lower=lower) + self.ndim = metric.ndim + self.axes = metric.axes + self.block = block + + kwargs = dict() + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseExpSquaredKernel(**kwargs) + super(ExpSquaredKernel, self).__init__([ + (None, base), ("metric", metric) + ]) + + # Common setup. + self.dirty = True + + @property + def block(self): + if not self.blocked: + return None + return list(zip(self.min_block, self.max_block)) + + @block.setter + def block(self, block): + if block is None: + self.blocked = False + self.min_block = -np.inf + np.zeros(len(self.axes)) + self.max_block = np.inf + np.zeros(len(self.axes)) + return + + block = np.atleast_2d(block) + if block.shape != (len(self.axes), 2): + raise ValueError("dimension mismatch in block specification") + self.blocked = True + self.min_block, self.max_block = map(np.array, zip(*block)) + + +class BaseMatern32Kernel (Model): + parameter_names = () + + +class Matern32Kernel (Kernel): + r""" + The Matern-3/2 kernel is stationary kernel where the value at a + given radius :math:`r^2` is given by: + + .. math:: + + k(r^2) = \left( 1+\sqrt{3\,r^2} \right)\, + \exp \left (-\sqrt{3\,r^2} \right ) + + + """ + + kernel_type = 11 + stationary = True + + def __init__(self, + metric=None, + metric_bounds=None, + lower=True, + block=None, + bounds=None, + ndim=1, + axes=None): + + if metric is None: + raise ValueError("missing required parameter 'metric'") + metric = Metric(metric, bounds=metric_bounds, ndim=ndim, + axes=axes, lower=lower) + self.ndim = metric.ndim + self.axes = metric.axes + self.block = block + + kwargs = dict() + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseMatern32Kernel(**kwargs) + super(Matern32Kernel, self).__init__([ + (None, base), ("metric", metric) + ]) + + # Common setup. + self.dirty = True + + @property + def block(self): + if not self.blocked: + return None + return list(zip(self.min_block, self.max_block)) + + @block.setter + def block(self, block): + if block is None: + self.blocked = False + self.min_block = -np.inf + np.zeros(len(self.axes)) + self.max_block = np.inf + np.zeros(len(self.axes)) + return + + block = np.atleast_2d(block) + if block.shape != (len(self.axes), 2): + raise ValueError("dimension mismatch in block specification") + self.blocked = True + self.min_block, self.max_block = map(np.array, zip(*block)) + + +class BasePolynomialKernel (Model): + parameter_names = ("log_sigma2", ) + + +class PolynomialKernel (Kernel): + r""" + A polynomial kernel + + .. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = + (\mathbf{x}_i \cdot \mathbf{x}_j + \sigma^2)^P + + :param order: + The power :math:`P`. This parameter is a *constant*; it is not + included in the parameter vector. + + :param log_sigma2: + The variance :math:`\sigma^2 > 0`. + + + """ + + kernel_type = 12 + stationary = False + + def __init__(self, + log_sigma2=None, + order=None, + bounds=None, + ndim=1, + axes=None): + + if order is None: + raise ValueError("missing required parameter 'order'") + self.order = order + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict(log_sigma2=log_sigma2, ) + if bounds is not None: + kwargs["bounds"] = bounds + base = BasePolynomialKernel(**kwargs) + super(PolynomialKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True + + +class BaseDotProductKernel (Model): + parameter_names = () + + +class DotProductKernel (Kernel): + r""" + The dot product kernel + + .. math:: + + k(\mathbf{x}_i,\,\mathbf{x}_j) = \mathbf{x}_i \cdot \mathbf{x}_j + + with no parameters. + + + """ + + kernel_type = 13 + stationary = False + + def __init__(self, + bounds=None, + ndim=1, + axes=None): + + self.subspace = Subspace(ndim, axes=axes) + self.ndim = self.subspace.ndim + self.axes = self.subspace.axes + + kwargs = dict() + if bounds is not None: + kwargs["bounds"] = bounds + base = BaseDotProductKernel(**kwargs) + super(DotProductKernel, self).__init__([ + (None, base), + ]) + + # Common setup. + self.dirty = True +