allow log hyperparameters

Author: dfm

docs/index.rst

@@ -22,6 +22,7 @@ User Guide
    user/hyper
    user/gp
    user/kernels
+   user/modeling
    user/solvers
 
 

docs/user/kernels.rst

@@ -19,33 +19,50 @@ The standard kernels fall into the following categories:
 sophisticated models and :ref:`new-kernels` explains how you would go about
 incorporating a custom kernel.
 
-**Note:** every kernel takes an optional ``ndim`` keyword that must be set to
-the number of input dimensions for your problem.
+Common parameters
+-----------------
+
+Every kernel accepts the two keyword arguments ``ndim`` and ``axes``. By
+default, kernels are only one dimensional so you must specify the ``ndim``
+argument if you want the kernel to work with higher dimensional inputs.
+By default, higher dimensional kernels are applied to every dimension but you
+can restrict the evaluation to a subspace using the ``axes`` argument.
+For example, if you have a 3 dimensional input space but you want one of the
+kernels to only act in the first dimension, you would do the following:
+
+.. code-block:: python
+
+    from george import kernels
+    kernel = 10.0 * kernels.Matern32Kernel(1.0, ndim=3, axes=0)
+
+Similarly, if you wanted the kernel to act on only the second and third
+dimensions, you could do something like:
+
+.. code-block:: python
+
+    kernel = 10.0 * kernels.ExpSquaredKernel([1.0, 0.5], ndim=3, axes=[1, 2])
+
 
 .. _implementation:
 
-Implementation Details
+Implementation details
 ----------------------
 
 It's worth understanding how these kernels are implemented.
 Most of the hard work is done at a low level (in C++) and the Python is only a
 thin wrapper to this functionality.
 This makes the code fast and consistent across interfaces but it also means
-that it isn't currently possible to implement new kernel functions efficiently
-without recompiling the code.
+that it isn't currently possible to implement new kernel functions without
+recompiling the code.
 Almost every kernel has hyperparameters that you can set to control its
-behavior and these can be accessed via the ``pars`` property.
-The values in this array are in the same order as you specified them when
-initializing the kernel and, in the case of composite kernels (see
-:ref:`combining-kernels`) the order goes from left to right.
-For example,
+behavior and these are controlled using the :ref:`modeling`.
 
 .. code-block:: python
 
     from george import kernels
 
     k = 2.0 * kernels.Matern32Kernel(5.0)
-    print(k.pars)
+    print(k.get_vector())
     # array([ 2.,  5.])
 
 
@@ -156,10 +173,4 @@ addition:
 Implementing New Kernels
 ------------------------
 
-Implementing custom kernels in George is a bit of a pain in the current
-version. For now, the only way to do it is with the :class:`PythonKernel`
-where you provide a Python function that computes the value of the kernel
-function at *a single pair of training points*.
-
-.. autoclass:: george.kernels.PythonKernel
-    :members:
+TO DO.

docs/user/modeling.rst

@@ -1,6 +1,6 @@
 .. _modeling:
 
-Modeling language
+Modeling protocol
 =================
 
 In order to make hyperparameter optimization, george comes with a modeling
@@ -142,9 +142,7 @@ implementation would be something like the following:
             return grad[self.unfrozen]
 
         def freeze_parameter(self, parameter_name):
-            names = self.parameter_names
-            self.unfrozen[names.index(parameter_name)] = False
+            self.unfrozen[self.parameter_names.index(parameter_name)] = False
 
         def thaw_parameter(self, parameter_name):
-            names = self.parameter_names
-            self.unfrozen[names.index(parameter_name)] = True
+            self.unfrozen[self.parameter_names.index(parameter_name)] = True

george/metrics.py

@@ -51,7 +51,7 @@ def __init__(self, metric, ndim=None, axes=None, lower=True):
                     raise ValueError("invalid (negative) metric")
                 for i, v in enumerate(metric):
                     self.parameter_names.append("ln_M_{0}_{0}".format(i))
-                    self.parameters.append(v)
+                    self.parameters.append(np.log(v))
 
             elif len(metric.shape) == 2:
                 self.metric_type = 2
@@ -83,7 +83,7 @@ def __init__(self, metric, ndim=None, axes=None, lower=True):
         else:
             self.metric_type = 0
             self.parameter_names.append("ln_M_0_0")
-            self.parameters.append(metric)
+            self.parameters.append(np.log(metric))
 
         self.parameters = np.array(self.parameters)
         self.unfrozen = np.ones_like(self.parameters, dtype=bool)

george/testing/test_kernels.py

@@ -47,50 +47,16 @@ def do_kernel_t(kernel, N=20, seed=123, eps=1.32e-6):
             "Gradient computation failed in dimension {0}".format(i)
 
 
-#
-# BASIC KERNELS
-#
-
-# def test_custom():
-#     def f(x1, x2, p):
-#         return np.exp(-0.5 * np.dot(x1, x2) / p[0])
-
-#     def g(x1, x2, p):
-#         arg = 0.5 * np.dot(x1, x2) / p[0]
-#         return np.exp(-arg) * arg / p[0]
-
-#     def wrong_g(x1, x2, p):
-#         arg = 0.5 * np.dot(x1, x2) / p[0]
-#         return 10 * np.exp(-arg) * arg / p[0]
-
-#     do_kernel_t(kernels.PythonKernel(f, g, pars=[0.5]))
-#     do_kernel_t(kernels.PythonKernel(f, g, pars=[0.1]))
-
-#     try:
-#         do_kernel_t(kernels.PythonKernel(f, wrong_g, pars=[0.5]))
-#     except AssertionError:
-#         pass
-#     else:
-#         assert False, "This test should fail"
-
-
-# def test_custom_numerical():
-#     def f(x1, x2, p):
-#         return np.exp(-0.5 * np.dot(x1, x2) / p[0])
-#     do_kernel_t(kernels.PythonKernel(f, pars=[0.5]))
-#     do_kernel_t(kernels.PythonKernel(f, pars=[10.0]))
-
-
 def test_constant():
-    do_kernel_t(kernels.ConstantKernel(0.1))
-    do_kernel_t(kernels.ConstantKernel(10.0, 2))
-    do_kernel_t(kernels.ConstantKernel(5.0, 5))
+    do_kernel_t(kernels.ConstantKernel(constant=0.1))
+    do_kernel_t(kernels.ConstantKernel(constant=10.0, ndim=2))
+    do_kernel_t(kernels.ConstantKernel(constant=5.0, ndim=5))
 
 
 def test_dot_prod():
     do_kernel_t(kernels.DotProductKernel())
-    do_kernel_t(kernels.DotProductKernel(2))
-    do_kernel_t(kernels.DotProductKernel(5))
+    do_kernel_t(kernels.DotProductKernel(ndim=2))
+    do_kernel_t(kernels.DotProductKernel(ndim=5, axes=0))
 
 
 #
@@ -145,35 +111,29 @@ def test_matern52():
 
 
 def test_rational_quadratic():
-    do_stationary_t(kernels.RationalQuadraticKernel, ln_alpha=np.log(1.0))
-    do_stationary_t(kernels.RationalQuadraticKernel, ln_alpha=np.log(0.1))
-    do_stationary_t(kernels.RationalQuadraticKernel, ln_alpha=np.log(10.0))
+    do_stationary_t(kernels.RationalQuadraticKernel, alpha=1.0)
+    do_stationary_t(kernels.RationalQuadraticKernel, alpha=0.1)
+    do_stationary_t(kernels.RationalQuadraticKernel, alpha=10.0)
 
 
-#
-# PERIODIC KERNELS
-#
-
 def test_cosine():
-    do_kernel_t(kernels.CosineKernel(1.0))
-    do_kernel_t(kernels.CosineKernel(0.5, ndim=2))
-    do_kernel_t(kernels.CosineKernel(0.5, ndim=2, axes=1))
-    do_kernel_t(kernels.CosineKernel(0.75, ndim=5, axes=[2, 3]))
+    do_kernel_t(kernels.CosineKernel(period=1.0))
+    do_kernel_t(kernels.CosineKernel(period=0.5, ndim=2))
+    do_kernel_t(kernels.CosineKernel(period=0.5, ndim=2, axes=1))
+    do_kernel_t(kernels.CosineKernel(period=0.75, ndim=5, axes=[2, 3]))
 
 
 def test_exp_sine2():
-    do_kernel_t(kernels.ExpSine2Kernel(0.4, 1.0))
-    do_kernel_t(kernels.ExpSine2Kernel(12., 0.5, ndim=2))
-    do_kernel_t(kernels.ExpSine2Kernel(17., 0.5, ndim=2, axes=1))
-    do_kernel_t(kernels.ExpSine2Kernel(13.7, -0.75, ndim=5, axes=[2, 3]))
-    do_kernel_t(kernels.ExpSine2Kernel(-0.7, 0.75, ndim=5, axes=[2, 3]))
-    do_kernel_t(kernels.ExpSine2Kernel(-10, 0.75))
-
+    do_kernel_t(kernels.ExpSine2Kernel(gamma=0.4, period=1.0))
+    do_kernel_t(kernels.ExpSine2Kernel(gamma=12., period=0.5, ndim=2))
+    do_kernel_t(kernels.ExpSine2Kernel(gamma=17., period=0.5, ndim=2, axes=1))
+    do_kernel_t(kernels.ExpSine2Kernel(gamma=13.7, ln_period=-0.75, ndim=5,
+                                       axes=[2, 3]))
+    do_kernel_t(kernels.ExpSine2Kernel(gamma=-0.7, period=0.75, ndim=5,
+                                       axes=[2, 3]))
+    do_kernel_t(kernels.ExpSine2Kernel(gamma=-10, period=0.75))
 
-#
-# COMBINING KERNELS
-#
 
 def test_combine():
-    do_kernel_t(12 * kernels.ExpSine2Kernel(0.4, 1.0, ndim=5) + 0.1)
+    do_kernel_t(12 * kernels.ExpSine2Kernel(gamma=0.4, period=1.0, ndim=5))
     do_kernel_t(12 * kernels.ExpSquaredKernel(0.4, ndim=3) + 0.1)

kernels/Constant.yml

@@ -11,7 +11,7 @@ doc: |
 
     where :math:`c` is a parameter.
 
-    :param lnconstant:
+    :param ln_constant:
         The logarithm of the constant value :math:`c` in the above equation.
 
 value: |

setup.py

@@ -170,6 +170,8 @@ def build_extension(self, ext):
     if len(kernel_specs):
         print("Compiling kernels")
         compile_kernels(kernel_specs)
+        if "kernels" in sys.argv:
+            sys.exit()
 
     # Check for the Cython source (development mode) and compile it if it
     # exists.

templates/kernels.py

@@ -101,15 +101,15 @@ def __len__(self):
 
     def __add__(self, b):
         if not hasattr(b, "is_kernel"):
-            return Sum(ConstantKernel(np.log(float(b)), ndim=self.ndim), self)
+            return Sum(ConstantKernel(constant=float(b), 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"):
-            return Product(ConstantKernel(np.log(float(b)), ndim=self.ndim),
+            return Product(ConstantKernel(constant=float(b), ndim=self.ndim),
                                           self)
         return Product(self, b)
 
@@ -281,6 +281,8 @@ class {{ spec.name }} (Kernel):
 
     def __init__(self,
                  {% for p in spec.params %}{{ p }}=None,
+                 {%- if p.startswith("ln_") %}
+                 {{ p[3:] }}=None,{% endif %}
                  {% endfor -%}
                  {% if spec.stationary -%}
                  metric=None,
@@ -292,11 +294,16 @@ def __init__(self,
         self._unfrozen = np.ones({{ spec.params|length }}, dtype=bool)
 
         {% for p in spec.params -%}
-        if {{ p }} is None:
-            raise ValueError("missing required parameter '{{ p }}'")
+        if {{ p }} is None{% if p.startswith("ln_") %} and {{ p[3:] }} is None{% endif %}:
+            raise ValueError("missing required parameter '{{ p }}'{% if p.startswith("ln_") %} or '{{ p[3:] }}'{% endif %}")
+        {%- if p.startswith("ln_") %}
+        elif {{ p }} is None:
+            if {{ p[3:] }} <= 0.0:
+                raise ValueError("invalid parameter '{{ p[3:] }} <= 0.0'")
+            {{ p }} = np.log({{ p[3:] }})
+        {%- endif %}
         self.{{ p }} = {{ p }}
-        {% endfor -%}
-
+        {% endfor %}
         {% if spec.stationary -%}
         if metric is None:
             raise ValueError("missing required parameter 'metric'")