|
1
|
|
|
# -*- coding: utf-8 -*- |
|
2
|
|
|
# vim:fileencoding=utf-8 |
|
3
|
|
|
# |
|
4
|
|
|
# Copyright (c) 2017-2018 Stefan Bender |
|
5
|
|
|
# |
|
6
|
|
|
# This module is part of sciapy. |
|
7
|
|
|
# sciapy is free software: you can redistribute it or modify |
|
8
|
|
|
# it under the terms of the GNU General Public License as published |
|
9
|
|
|
# by the Free Software Foundation, version 2. |
|
10
|
|
|
# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html. |
|
11
|
1 |
|
"""SCIAMACHY regression tool Gaussian process kernel options |
|
12
|
|
|
|
|
13
|
|
|
Gaussian process kernel standard parameters for the command line tool. |
|
14
|
|
|
""" |
|
15
|
|
|
|
|
16
|
1 |
|
import numpy as np |
|
17
|
|
|
|
|
18
|
1 |
|
import george |
|
19
|
1 |
|
from george import kernels |
|
20
|
|
|
|
|
21
|
1 |
|
from celerite import terms |
|
22
|
|
|
|
|
23
|
1 |
|
__all__ = ["george_kernels", "george_solvers", |
|
24
|
|
|
"celerite_terms", |
|
25
|
|
|
"setup_george_kernel", "setup_celerite_terms"] |
|
26
|
|
|
|
|
27
|
1 |
|
george_kernels = { |
|
28
|
|
|
"Exp2": kernels.ExpSquaredKernel(10**2), |
|
29
|
|
|
"Exp2ESin2": (kernels.ExpSquaredKernel(10**2) * |
|
30
|
|
|
kernels.ExpSine2Kernel(2 / 1.3**2, 1.0)), |
|
31
|
|
|
"ESin2": kernels.ExpSine2Kernel(2 / 1.3**2, 1.0), |
|
32
|
|
|
"27dESin2": kernels.ExpSine2Kernel(2 / 1.3**2, 27.0 / 365.25), |
|
33
|
|
|
"RatQ": kernels.RationalQuadraticKernel(0.8, 0.1**2), |
|
34
|
|
|
"Mat32": kernels.Matern32Kernel((0.5)**2), |
|
35
|
|
|
"Exp": kernels.ExpKernel((0.5)**2), |
|
36
|
|
|
# "W": kernels.WhiteKernel, # deprecated, delegated to `white_noise` |
|
37
|
|
|
"B": kernels.ConstantKernel, |
|
38
|
|
|
} |
|
39
|
1 |
|
george_solvers = { |
|
40
|
|
|
"basic": george.BasicSolver, |
|
41
|
|
|
"HODLR": george.HODLRSolver, |
|
42
|
|
|
} |
|
43
|
|
|
|
|
44
|
1 |
|
celerite_terms = { |
|
45
|
|
|
"N": terms.Term(), |
|
46
|
|
|
"B": terms.RealTerm(log_a=-6., log_c=-np.inf, |
|
47
|
|
|
bounds={"log_a": [-30, 30], |
|
48
|
|
|
"log_c": [-np.inf, np.inf]}), |
|
49
|
|
|
"W": terms.JitterTerm(log_sigma=-25, |
|
50
|
|
|
bounds={"log_sigma": [-30, 30]}), |
|
51
|
|
|
"Mat32": terms.Matern32Term( |
|
52
|
|
|
log_sigma=1., |
|
53
|
|
|
log_rho=1., |
|
54
|
|
|
bounds={"log_sigma": [-30, 30], |
|
55
|
|
|
# The `celerite` version of the Matern-3/2 |
|
56
|
|
|
# kernel has problems with very large `log_rho` |
|
57
|
|
|
# values. -7.4 is empirical. |
|
58
|
|
|
"log_rho": [-7.4, 15]}), |
|
59
|
|
|
"SHO0": terms.SHOTerm(log_S0=-6, log_Q=1.0 / np.sqrt(2.), log_omega0=0., |
|
60
|
|
|
bounds={"log_S0": [-30, 30], |
|
61
|
|
|
"log_Q": [-30, 30], |
|
62
|
|
|
"log_omega0": [-30, 30]}), |
|
63
|
|
|
"SHO1": terms.SHOTerm(log_S0=-6, log_Q=-2., log_omega0=0., |
|
64
|
|
|
bounds={"log_S0": [-10, 10], |
|
65
|
|
|
"log_omega0": [-10, 10]}), |
|
66
|
|
|
"SHO2": terms.SHOTerm(log_S0=-6, log_Q=0.5, log_omega0=0., |
|
67
|
|
|
bounds={"log_S0": [-10, 10], |
|
68
|
|
|
"log_Q": [-10, 10], |
|
69
|
|
|
"log_omega0": [-10, 10]}), |
|
70
|
|
|
# see Foreman-Mackey et al. 2017, AJ 154, 6, pp. 220 |
|
71
|
|
|
# doi: 10.3847/1538-3881/aa9332 |
|
72
|
|
|
# Eq. (53) |
|
73
|
|
|
"SHO3": terms.SHOTerm(log_S0=-6, log_Q=0., log_omega0=0., |
|
74
|
|
|
bounds={"log_S0": [-15, 5], |
|
75
|
|
|
"log_Q": [-10, 10], |
|
76
|
|
|
"log_omega0": [-10, 10]}) * |
|
77
|
|
|
terms.SHOTerm(log_S0=1, log_Q=1.0 / np.sqrt(2.), log_omega0=0., |
|
78
|
|
|
bounds={"log_omega0": [-5, 5]}), |
|
79
|
|
|
} |
|
80
|
1 |
|
_celerite_terms_freeze_params = { |
|
81
|
|
|
"N": [], |
|
82
|
|
|
"B": ["log_c"], |
|
83
|
|
|
"W": [], |
|
84
|
|
|
"Mat32": [], |
|
85
|
|
|
"SHO0": ["log_Q", "log_omega0"], |
|
86
|
|
|
"SHO1": ["log_Q"], |
|
87
|
|
|
"SHO2": [], |
|
88
|
|
|
"SHO3": ["k2:log_S0", "k2:log_Q"], |
|
89
|
|
|
} |
|
90
|
|
|
|
|
91
|
|
|
|
|
92
|
1 |
|
def setup_george_kernel(kernelnames, kernel_base=1, fit_bias=False): |
|
93
|
|
|
"""Setup the Gaussian Process kernel for george |
|
94
|
|
|
|
|
95
|
|
|
Parameters |
|
96
|
|
|
---------- |
|
97
|
|
|
kernelnames : list of str |
|
98
|
|
|
List of abbreviated names for the kernels, choices: |
|
99
|
|
|
'Exp2' for a `ExpSquaredKernel`, 'ESin2' for a `ExpSine2Kernel`, |
|
100
|
|
|
'Exp2ESin2' for a `ExpSquaredKernel` multiplied by a `ExpSine2Kernel`, |
|
101
|
|
|
'RatQ' for a `RationalQuadraticKernel`, 'Mat32' for a `Matern32Kernel`, |
|
102
|
|
|
'Exp' for `ExpKernel`, and 'B' for a `ConstantKernel` (bias). |
|
103
|
|
|
kernel_base : float, optional |
|
104
|
|
|
The initial "strength" of the kernels. |
|
105
|
|
|
fit_bias : bool, optional |
|
106
|
|
|
Adds a `ConstantKernel` if kernel does not already contain one. |
|
107
|
|
|
|
|
108
|
|
|
Returns |
|
109
|
|
|
------- |
|
110
|
|
|
name : The kernel names concatenated by an underscore and prepended by '_gp' |
|
111
|
|
|
kernel : The covariance kernel for use with george.GP |
|
112
|
|
|
""" |
|
113
|
1 |
|
kernel = None |
|
114
|
1 |
|
kname = "_gp" |
|
115
|
1 |
|
for kn in kernelnames: |
|
116
|
1 |
|
krn = george_kernels.get(kn, None) |
|
117
|
1 |
|
if krn is None: |
|
118
|
|
|
# not found in the list of available kernels |
|
119
|
|
|
continue |
|
120
|
1 |
|
if kn in ["B", "W"]: |
|
121
|
|
|
# don't scale the constant or white kernels |
|
122
|
|
|
krnl = krn(0.25 * kernel_base) |
|
123
|
|
|
else: |
|
124
|
1 |
|
krnl = kernel_base * krn |
|
125
|
1 |
|
kernel = kernel + krnl if hasattr(kernel, "is_kernel") else krnl |
|
126
|
1 |
|
kname += "_" + kn |
|
127
|
|
|
|
|
128
|
1 |
|
if fit_bias and "B" not in kernelnames: |
|
129
|
|
|
krnl = kernels.ConstantKernel(kernel_base) |
|
130
|
|
|
kernel = kernel + krnl if hasattr(kernel, "is_kernel") else krnl |
|
131
|
|
|
kname += "_b" |
|
132
|
|
|
|
|
133
|
1 |
|
return kname, kernel |
|
134
|
|
|
|
|
135
|
|
|
|
|
136
|
1 |
|
def setup_celerite_terms(termnames, fit_bias=False, fit_white=False): |
|
137
|
|
|
"""Setup the Gaussian Process terms for celerite |
|
138
|
|
|
|
|
139
|
|
|
Parameters |
|
140
|
|
|
---------- |
|
141
|
|
|
termnames : list of str |
|
142
|
|
|
List of abbreviated names for the `celerite.terms`, choices: |
|
143
|
|
|
'N' for an empty `Term`, 'B' for a constant term (bias), |
|
144
|
|
|
'W' for a `JitterTerm` (white noise), 'Mat32' for a `Matern32Term`, |
|
145
|
|
|
'SHO0'...'SHO3' for `SHOTerm`s (harmonic oscillators) with different |
|
146
|
|
|
frozen parameters. |
|
147
|
|
|
fit_bias : bool, optional |
|
148
|
|
|
Adds a constant term (`RealTerm` with log_c fixed to -np.inf) if the |
|
149
|
|
|
terms do not already contain one. |
|
150
|
|
|
fit_white : bool, optional |
|
151
|
|
|
Adds a `JitterTerm` if not already included. |
|
152
|
|
|
|
|
153
|
|
|
Returns |
|
154
|
|
|
------- |
|
155
|
|
|
name : The term names concatenated by an underscore and prepended by '_cel' |
|
156
|
|
|
terms : The covariance terms for use with celerite.GP |
|
157
|
|
|
""" |
|
158
|
|
|
# celerite terms |
|
159
|
1 |
|
cel_terms = None |
|
160
|
1 |
|
cel_name = "_cel" |
|
161
|
1 |
|
for tn in termnames: |
|
162
|
1 |
|
trm = celerite_terms.get(tn, None) |
|
163
|
1 |
|
if trm is None: |
|
164
|
|
|
# not found in the list of available terms |
|
165
|
1 |
|
continue |
|
166
|
1 |
|
for freeze_p in _celerite_terms_freeze_params[tn]: |
|
167
|
|
|
trm.freeze_parameter(freeze_p) |
|
168
|
1 |
|
cel_terms = cel_terms + trm if cel_terms is not None else trm |
|
169
|
1 |
|
cel_name += "_" + tn |
|
170
|
|
|
|
|
171
|
|
|
# avoid double initialisation of White and Constant kernels |
|
172
|
|
|
# if already set above |
|
173
|
1 |
|
if fit_white and "W" not in termnames: |
|
174
|
1 |
|
trm = celerite_terms["W"] |
|
175
|
1 |
|
cel_terms = cel_terms + trm if cel_terms is not None else trm |
|
176
|
1 |
|
cel_name += "_w" |
|
177
|
1 |
|
if fit_bias and "B" not in termnames: |
|
178
|
|
|
trm = celerite_terms["B"] |
|
179
|
|
|
trm.freeze_parameter("log_c") |
|
180
|
|
|
cel_terms = cel_terms + trm if cel_terms is not None else trm |
|
181
|
|
|
cel_name += "_b" |
|
182
|
|
|
|
|
183
|
1 |
|
if cel_terms is None: |
|
184
|
1 |
|
cel_terms = terms.Term() |
|
185
|
|
|
|
|
186
|
|
|
return cel_name, cel_terms |
|
187
|
|
|
|
st-bender /
sciapy
