17445970354

Committed 03 Sep 2025 09:00PM UTC coverage: 87.165%. First build

Build # 17445970354

Build Type

Pull #395

github

Committed by

bjfultn

Commit Message

Remove debug output from GP parameter validation

Pull Request Pull Request #395: testing new ci

Run Details

169 of 236 new or added lines in 6 files covered. (71.61%)

3735 of 4285 relevant lines covered (87.16%)

0.87 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

87.97

/radvel/gp.py

import sys
import radvel
import scipy
from scipy import spatial
import abc
import numpy as np
import warnings

warnings.simplefilter('once')

# implemented kernels & examples of their associated hyperparameters
KERNELS = {
    "SqExp": ['gp_length', 'gp_amp'],
    "Per": ['gp_per', 'gp_length', 'gp_amp'],
    "QuasiPer": ['gp_per', 'gp_perlength', 'gp_explength', 'gp_amp'],
    "Celerite": ['gp_B', 'gp_C', 'gp_L', 'gp_Prot']
}

if sys.version_info[0] < 3:
    ABC = abc.ABCMeta('ABC', (), {})
else:
    ABC = abc.ABC


# celerite is an optional dependency
def _try_celerite():
    try:
        import celerite
        from celerite.solver import CholeskySolver
        return True
    except ImportError:
        warnings.warn("celerite not installed. GP kernals using celerite will not work. \
Try installing celerite using 'pip install celerite'", ImportWarning)
        return False


_has_celerite = _try_celerite()
if _has_celerite:
    import celerite
    from celerite.solver import CholeskySolver


class Kernel(ABC):
    """
    Abstract base class to store kernel info and compute covariance matrix.
    All kernel objects inherit from this class.

    Note:
        To implement your own kernel, create a class that inherits
        from this class. It should have hyperparameters that follow
        the name scheme 'gp_NAME_SUFFIX'.

    """

    @abc.abstractproperty
    def name(self):
        pass

    @abc.abstractmethod
    def compute_distances(self, x1, x2):
        pass

    @abc.abstractmethod
    def compute_covmatrix(self, errors):
        pass


class SqExpKernel(Kernel):
    """
    Class that computes and stores a squared exponential kernel matrix.
    An arbitrary element, :math:`C_{ij}`, of the matrix is:

    .. math::

        C_{ij} = \\eta_1^2 * exp( \\frac{ -|t_i - t_j|^2 }{ \\eta_2^2 } )

    Args:
        hparams (dict of radvel.Parameter): dictionary containing
            radvel.Parameter objects that are GP hyperparameters
            of this kernel. Must contain exactly two objects, 'gp_length*'
            and 'gp_amp*', where * is a suffix identifying
            these hyperparameters with a likelihood object.

    """

    @property
    def name(self):
        return "SqExp"

    def __init__(self, hparams):
        self.covmatrix = None
        self.hparams = {}
        for par in hparams:
            if par.startswith('gp_length'):
                self.hparams['gp_length'] = hparams[par]
            if par.startswith('gp_amp'):
                self.hparams['gp_amp'] = hparams[par]

        assert len(hparams) == 2, \
            "SqExpKernel requires exactly 2 hyperparameters with names" \
            + "'gp_length*' and 'gp_amp*'."

        try:
            self.hparams['gp_length'].value
            self.hparams['gp_amp'].value
        except KeyError:
            raise KeyError("SqExpKernel requires hyperparameters 'gp_length*'" \
                           + " and 'gp_amp*'.")
        except AttributeError:
            raise AttributeError("SqExpKernel requires dictionary of" \
                                 + " radvel.Parameter objects as input.")

    def __repr__(self):
        length = self.hparams['gp_length'].value
        amp = self.hparams['gp_amp'].value
        return "SqExp Kernel with length: {}, amp: {}".format(length, amp)

    def compute_distances(self, x1, x2):
        X1 = np.array([x1]).T
        X2 = np.array([x2]).T
        self.dist = scipy.spatial.distance.cdist(X1, X2, 'sqeuclidean')

    def compute_covmatrix(self, errors):
        """ Compute the covariance matrix, and optionally add errors along
            the diagonal.

            Args:
                errors (float or numpy array): If covariance matrix is non-square,
                    this arg must be set to 0. If covariance matrix is square,
                    this can be a numpy array of observational errors and jitter
                    added in quadrature.
        """
        length = self.hparams['gp_length'].value
        amp = self.hparams['gp_amp'].value

        K = amp**2 * np.exp(-self.dist/(length**2))

        self.covmatrix = K
        # add errors along the diagonal
        try:
            self.covmatrix += (errors**2) * np.identity(K.shape[0])
        except ValueError: # errors can't be added along diagonal to a non-square array
            pass

        return self.covmatrix


class PerKernel(Kernel):
    """
    Class that computes and stores a periodic kernel matrix.
    An arbitrary element, :math:`C_{ij}`, of the matrix is:

    .. math::

        C_{ij} = \\eta_1^2 * exp( \\frac{ -\\sin^2(\\frac{ \\pi|t_i-t_j| }{ \\eta_3^2 } ) }{ 2\\eta_2^2 } )

    Args:
        hparams (dict of radvel.Parameter): dictionary containing
            radvel.Parameter objects that are GP hyperparameters
            of this kernel. Must contain exactly three objects, 'gp_length*',
            'gp_amp*', and 'gp_per*', where * is a suffix identifying
            these hyperparameters with a likelihood object.

    """

    @property
    def name(self):
        return "Per"

    def __init__(self, hparams):
        self.covmatrix = None
        self.hparams = {}
        for par in hparams:
            if par.startswith('gp_length'):
                self.hparams['gp_length'] = hparams[par]
            if par.startswith('gp_amp'):
                self.hparams['gp_amp'] = hparams[par]
            if par.startswith('gp_per'):
                self.hparams['gp_per'] = hparams[par]

        assert len(hparams) == 3, \
            "PerKernel requires exactly 3 hyperparameters with names 'gp_length*'," \
            + " 'gp_amp*', and 'gp_per*'."

        try:
            self.hparams['gp_length'].value
            self.hparams['gp_amp'].value
            self.hparams['gp_per'].value
        except KeyError:
            raise KeyError("PerKernel requires hyperparameters 'gp_length*'," \
                           + " 'gp_amp*', and 'gp_per*'.")
        except AttributeError:
            raise AttributeError("PerKernel requires dictionary of " \
                                 + "radvel.Parameter objects as input.")

    def __repr__(self):
        length = self.hparams['gp_length'].value
        amp = self.hparams['gp_amp'].value
        per = self.hparams['gp_per'].value
        return "Per Kernel with length: {}, amp: {}, per: {}".format(
            length, amp, per
        )

    def compute_distances(self, x1, x2):
        X1 = np.array([x1]).T
        X2 = np.array([x2]).T
        self.dist = scipy.spatial.distance.cdist(X1, X2, 'euclidean')

    def compute_covmatrix(self, errors):
        """ Compute the covariance matrix, and optionally add errors along
            the diagonal.

            Args:
                errors (float or numpy array): If covariance matrix is non-square,
                    this arg must be set to 0. If covariance matrix is square,
                    this can be a numpy array of observational errors and jitter
                    added in quadrature.
        """
        length= self.hparams['gp_length'].value
        amp = self.hparams['gp_amp'].value
        per = self.hparams['gp_per'].value

        K = amp**2 * np.exp(-np.sin(np.pi*self.dist/per)**2. / (2.*length**2))
        self.covmatrix = K
        # add errors along the diagonal
        try:
            self.covmatrix += (errors**2) * np.identity(K.shape[0])
        except ValueError:  # errors can't be added along diagonal to a non-square array
            pass

        return self.covmatrix

class QuasiPerKernel(Kernel):
    """
    Class that computes and stores a quasi periodic kernel matrix.
    An arbitrary element, :math:`C_{ij}`, of the matrix is:

    .. math::

        C_{ij} = \\eta_1^2 * exp( \\frac{ -|t_i - t_j|^2 }{ \\eta_2^2 } -
                 \\frac{ \\sin^2(\\frac{ \\pi|t_i-t_j| }{ \\eta_3 } ) }{ 2\\eta_4^2 } )

    Args:
        hparams (dict of radvel.Parameter): dictionary containing
            radvel.Parameter objects that are GP hyperparameters
            of this kernel. Must contain exactly four objects, 'gp_explength*',
            'gp_amp*', 'gp_per*', and 'gp_perlength*', where * is a suffix
            identifying these hyperparameters with a likelihood object.

    """
    @property
    def name(self):
        return "QuasiPer"

    def __init__(self, hparams):
        self.covmatrix = None
        self.hparams = {}
        for par in hparams:
            if par.startswith('gp_perlength'):
                self.hparams['gp_perlength'] = hparams[par]
            if par.startswith('gp_amp'):
                self.hparams['gp_amp'] = hparams[par]
            if par.startswith('gp_per') and not 'length' in par:
                self.hparams['gp_per'] = hparams[par]
            if par.startswith('gp_explength'):
                self.hparams['gp_explength'] = hparams[par]

        assert len(hparams) == 4, \
            "QuasiPerKernel requires exactly 4 hyperparameters with names" \
            + " 'gp_perlength*', 'gp_amp*', 'gp_per*', and 'gp_explength*'."

        try:
            self.hparams['gp_perlength'].value
            self.hparams['gp_amp'].value
            self.hparams['gp_per'].value
            self.hparams['gp_explength'].value
        except KeyError:
            raise KeyError("QuasiPerKernel requires hyperparameters" \
                           + " 'gp_perlength*', 'gp_amp*', 'gp_per*', " \
                           + "and 'gp_explength*'.")
        except AttributeError:
            raise AttributeError("QuasiPerKernel requires dictionary of" \
                                 + " radvel.Parameter objects as input.")

    def __repr__(self):
        perlength = self.hparams['gp_perlength'].value
        amp = self.hparams['gp_amp'].value
        per = self.hparams['gp_per'].value
        explength = self.hparams['gp_explength'].value

        msg = (
            "QuasiPer Kernel with amp: {}, per length: {}, per: {}, "
            "exp length: {}"
        ).format(amp, perlength, per, explength)
        return msg

    def compute_distances(self, x1, x2):
        X1 = np.array([x1]).T
        X2 = np.array([x2]).T
        self.dist_p = scipy.spatial.distance.cdist(X1, X2, 'euclidean')
        self.dist_se = scipy.spatial.distance.cdist(X1, X2, 'sqeuclidean')

    def compute_covmatrix(self, errors):
        """ Compute the covariance matrix, and optionally add errors along
            the diagonal.

            Args:
                errors (float or numpy array): If covariance matrix is non-square,
                    this arg must be set to 0. If covariance matrix is square,
                    this can be a numpy array of observational errors and jitter
                    added in quadrature.
        """
        perlength = self.hparams['gp_perlength'].value
        amp = self.hparams['gp_amp'].value
        per = self.hparams['gp_per'].value
        explength = self.hparams['gp_explength'].value

        # Check for problematic parameter values
        # Allow 0 values for model comparison (when GP is disabled)
        if not np.isfinite(amp) or amp < 0:
            raise ValueError(f"Invalid gp_amp value: {amp}")
        if not np.isfinite(per) or per < 0:
            raise ValueError(f"Invalid gp_per value: {per}")
        if not np.isfinite(explength) or explength < 0:
            raise ValueError(f"Invalid gp_explength value: {explength}")
        if not np.isfinite(perlength) or perlength < 0:
            raise ValueError(f"Invalid gp_perlength value: {perlength}")
        
        # If any GP parameter is 0, this is likely a "no GP" model comparison
        # In this case, return a zero covariance matrix
        if amp == 0 or per == 0 or explength == 0 or perlength == 0:
            n = len(self.dist_se)
            K = np.zeros((n, n))
            if isinstance(errors, (int, float)) and errors == 0:
                return K
            else:
                K[np.diag_indices_from(K)] += errors**2
                return K

        # Compute covariance matrix step by step
        exp_term1 = -self.dist_se/(explength**2)
        exp_term2 = (-np.sin(np.pi*self.dist_p/per)**2.) / (2.*perlength**2)
        
        # Clamp extreme values to prevent overflow
        exp_term1 = np.clip(exp_term1, -100, 100)
        exp_term2 = np.clip(exp_term2, -100, 100)
        
        # Compute exponential terms
        exp1 = np.exp(exp_term1)
        exp2 = np.exp(exp_term2)
        
        # Compute final covariance matrix
        K = np.array(amp**2 * exp1 * exp2)

        self.covmatrix = K

        # add errors along the diagonal
        try:
            self.covmatrix += (errors**2) * np.identity(K.shape[0])
        except ValueError:  # errors can't be added along diagonal to a non-square array
            pass

        return self.covmatrix


class CeleriteKernel(Kernel):
    """
    Class that computes and stores a matrix approximating the quasi-periodic
    kernel.

    See `radvel/example_planets/k2-131_celerite.py` for an example of a setup
    file that uses this Kernel object.

    See celerite.readthedocs.io and Foreman-Mackey et al. 2017. AJ, 154, 220
    (equation 56) for more details.

    An arbitrary element, :math:`C_{ij}`, of the matrix is:

    .. math::

        C_{ij} = B/(2+C) * exp( -|t_i - t_j| / L) * (\\cos(\\frac{ 2\\pi|t_i-t_j| }{ P_{rot} }) + (1+C) )

    Args:
        hparams (dict of radvel.Parameter): dictionary containing
            radvel.Parameter objects that are GP hyperparameters
            of this kernel. Must contain exactly four objects, 'gp_B*',
            'gp_C*', 'gp_L*', and 'gp_Prot*', where * is a suffix
            identifying these hyperparameters with a likelihood object.
    """

    @property
    def name(self):
        return "Celerite"

    def __init__(self, hparams):

        self.hparams = {}
        for par in hparams:
            if par.startswith('gp_B'):
                self.hparams['gp_B'] = hparams[par]
            if par.startswith('gp_C'):
                self.hparams['gp_C'] = hparams[par]
            if par.startswith('gp_L'):
                self.hparams['gp_L'] = hparams[par]
            if par.startswith('gp_Prot'):
                self.hparams['gp_Prot'] = hparams[par]

        assert len(self.hparams) == 4, """
CeleriteKernel requires exactly 4 hyperparameters with names 'gp_B', 'gp_C', 'gp_L', and 'gp_Prot'.
        """

        try:
            self.hparams['gp_Prot'].value
            self.hparams['gp_C'].value
            self.hparams['gp_B'].value
            self.hparams['gp_L'].value
        except KeyError:
            raise KeyError("""
CeleriteKernel requires hyperparameters 'gp_B*', 'gp_C*', 'gp_L', and 'gp_Prot*'.
                """)
        except AttributeError:
            raise AttributeError("CeleriteKernel requires dictionary of radvel.Parameter objects as input.")

    # get arrays of real and complex parameters
    def compute_real_and_complex_hparams(self):

        self.real = np.zeros((1, 4))
        self.complex = np.zeros((1, 4))

        B = self.hparams['gp_B'].value
        C = self.hparams['gp_C'].value
        L = self.hparams['gp_L'].value
        Prot = self.hparams['gp_Prot'].value

        # Foreman-Mackey et al. (2017) eq 56
        self.real[0,0] = B*(1+C)/(2+C)
        self.real[0,2] = 1/L
        self.complex[0,0] = B/(2+C)
        self.complex[0,1] = 0.
        self.complex[0,2] = 1/L
        self.complex[0,3] = 2*np.pi/Prot

    def __repr__(self):

        B = self.hparams['gp_B'].value
        C = self.hparams['gp_C'].value
        L = self.hparams['gp_L'].value
        Prot = self.hparams['gp_Prot'].value

        msg = (
            "Celerite Kernel with B = {}, C = {}, L = {}, Prot = {}."
        ).format(B, C, L, Prot)
        return msg

    def compute_distances(self, x1, x2):
        """
        The celerite.solver.CholeskySolver object does
        not require distances to be precomputed, so
        this method has been co-opted to define some
        unchanging variables.
        """
        self.x = x1

        # blank matrices (corresponding to Cholesky decomp of kernel) needed for celerite solver
        self.A = np.empty(0)
        self.U = np.empty((0,0))
        self.V = self.U


    def compute_covmatrix(self, errors):
        """ Compute the Cholesky decomposition of a celerite kernel

            Args:
                errors (array of float): observation errors and jitter added
                    in quadrature

            Returns:
                celerite.solver.CholeskySolver: the celerite solver object,
                with Cholesky decomposition computed.
        """
        # initialize celerite solver object
        solver = CholeskySolver()

        self.compute_real_and_complex_hparams()
        solver.compute(
            0., self.real[:,0], self.real[:,2],
            self.complex[:,0], self.complex[:,1],
            self.complex[:,2], self.complex[:,3],
            self.A, self.U, self.V,
            self.x, errors**2
        )

        return solver

1	import sys	1✔
2	import radvel	1✔
3	import scipy	1✔
4	from scipy import spatial	1✔
5	import abc	1✔
6	import numpy as np	1✔
7	import warnings	1✔
8
9	warnings.simplefilter('once')	1✔
10
11	# implemented kernels & examples of their associated hyperparameters
12	KERNELS = {	1✔
13	"SqExp": ['gp_length', 'gp_amp'],
14	"Per": ['gp_per', 'gp_length', 'gp_amp'],
15	"QuasiPer": ['gp_per', 'gp_perlength', 'gp_explength', 'gp_amp'],
16	"Celerite": ['gp_B', 'gp_C', 'gp_L', 'gp_Prot']
17	}
18
19	if sys.version_info[0] < 3:	1✔
20	ABC = abc.ABCMeta('ABC', (), {})	×
21	else:
22	ABC = abc.ABC	1✔
23
24
25	# celerite is an optional dependency
26	def _try_celerite():	1✔
27	try:	1✔
28	import celerite	1✔
29	from celerite.solver import CholeskySolver	1✔
30	return True	1✔
31	except ImportError:	×
32	warnings.warn("celerite not installed. GP kernals using celerite will not work. \	×
33	Try installing celerite using 'pip install celerite'", ImportWarning)
34	return False	×
35
36
37	_has_celerite = _try_celerite()	1✔
38	if _has_celerite:	1✔
39	import celerite	1✔
40	from celerite.solver import CholeskySolver	1✔
41
42
43	class Kernel(ABC):	1✔
44	"""
45	Abstract base class to store kernel info and compute covariance matrix.
46	All kernel objects inherit from this class.
47
48	Note:
49	To implement your own kernel, create a class that inherits
50	from this class. It should have hyperparameters that follow
51	the name scheme 'gp_NAME_SUFFIX'.
52
53	"""
54
55	@abc.abstractproperty	1✔
56	def name(self):	1✔
57	pass	×
58
59	@abc.abstractmethod	1✔
60	def compute_distances(self, x1, x2):	1✔
61	pass	×
62
63	@abc.abstractmethod	1✔
64	def compute_covmatrix(self, errors):	1✔
65	pass	×
66
67
68	class SqExpKernel(Kernel):	1✔
69	"""
70	Class that computes and stores a squared exponential kernel matrix.
71	An arbitrary element, :math:`C_{ij}`, of the matrix is:
72
73	.. math::
74
75	C_{ij} = \\eta_1^2 * exp( \\frac{ -\|t_i - t_j\|^2 }{ \\eta_2^2 } )
76
77	Args:
78	hparams (dict of radvel.Parameter): dictionary containing
79	radvel.Parameter objects that are GP hyperparameters
80	of this kernel. Must contain exactly two objects, 'gp_length*'
81	and 'gp_amp', where is a suffix identifying
82	these hyperparameters with a likelihood object.
83
84	"""
85
86	@property	1✔
87	def name(self):	1✔
88	return "SqExp"	×
89
90	def __init__(self, hparams):	1✔
91	self.covmatrix = None	1✔
92	self.hparams = {}	1✔
93	for par in hparams:	1✔
94	if par.startswith('gp_length'):	1✔
95	self.hparams['gp_length'] = hparams[par]	1✔
96	if par.startswith('gp_amp'):	1✔
97	self.hparams['gp_amp'] = hparams[par]	1✔
98
99	assert len(hparams) == 2, \	1✔
100	"SqExpKernel requires exactly 2 hyperparameters with names" \
101	+ "'gp_length' and 'gp_amp'."
102
103	try:	1✔
104	self.hparams['gp_length'].value	1✔
105	self.hparams['gp_amp'].value	1✔
106	except KeyError:	1✔
107	raise KeyError("SqExpKernel requires hyperparameters 'gp_length*'" \	×
108	+ " and 'gp_amp*'.")
109	except AttributeError:	1✔
110	raise AttributeError("SqExpKernel requires dictionary of" \	1✔
111	+ " radvel.Parameter objects as input.")
112
113	def __repr__(self):	1✔
114	length = self.hparams['gp_length'].value	1✔
115	amp = self.hparams['gp_amp'].value	1✔
116	return "SqExp Kernel with length: {}, amp: {}".format(length, amp)	1✔
117
118	def compute_distances(self, x1, x2):	1✔
119	X1 = np.array([x1]).T	1✔
120	X2 = np.array([x2]).T	1✔
121	self.dist = scipy.spatial.distance.cdist(X1, X2, 'sqeuclidean')	1✔
122
123	def compute_covmatrix(self, errors):	1✔
124	""" Compute the covariance matrix, and optionally add errors along
125	the diagonal.
126
127	Args:
128	errors (float or numpy array): If covariance matrix is non-square,
129	this arg must be set to 0. If covariance matrix is square,
130	this can be a numpy array of observational errors and jitter
131	added in quadrature.
132	"""
133	length = self.hparams['gp_length'].value	1✔
134	amp = self.hparams['gp_amp'].value	1✔
135
136	K = amp*2 np.exp(-self.dist/(length**2))	1✔
137
138	self.covmatrix = K	1✔
139	# add errors along the diagonal
140	try:	1✔
141	self.covmatrix += (errors*2) np.identity(K.shape[0])	1✔
142	except ValueError: # errors can't be added along diagonal to a non-square array	×
143	pass	×
144
145	return self.covmatrix	1✔
146
147
148	class PerKernel(Kernel):	1✔
149	"""
150	Class that computes and stores a periodic kernel matrix.
151	An arbitrary element, :math:`C_{ij}`, of the matrix is:
152
153	.. math::
154
155	C_{ij} = \\eta_1^2 * exp( \\frac{ -\\sin^2(\\frac{ \\pi\|t_i-t_j\| }{ \\eta_3^2 } ) }{ 2\\eta_2^2 } )
156
157	Args:
158	hparams (dict of radvel.Parameter): dictionary containing
159	radvel.Parameter objects that are GP hyperparameters
160	of this kernel. Must contain exactly three objects, 'gp_length*',
161	'gp_amp', and 'gp_per', where * is a suffix identifying
162	these hyperparameters with a likelihood object.
163
164	"""
165
166	@property	1✔
167	def name(self):	1✔
168	return "Per"	×
169
170	def __init__(self, hparams):	1✔
171	self.covmatrix = None	1✔
172	self.hparams = {}	1✔
173	for par in hparams:	1✔
174	if par.startswith('gp_length'):	1✔
175	self.hparams['gp_length'] = hparams[par]	1✔
176	if par.startswith('gp_amp'):	1✔
177	self.hparams['gp_amp'] = hparams[par]	1✔
178	if par.startswith('gp_per'):	1✔
179	self.hparams['gp_per'] = hparams[par]	1✔
180
181	assert len(hparams) == 3, \	1✔
182	"PerKernel requires exactly 3 hyperparameters with names 'gp_length*'," \
183	+ " 'gp_amp', and 'gp_per'."
184
185	try:	1✔
186	self.hparams['gp_length'].value	1✔
187	self.hparams['gp_amp'].value	1✔
188	self.hparams['gp_per'].value	1✔
189	except KeyError:	1✔
190	raise KeyError("PerKernel requires hyperparameters 'gp_length*'," \	×
191	+ " 'gp_amp', and 'gp_per'.")
192	except AttributeError:	1✔
193	raise AttributeError("PerKernel requires dictionary of " \	1✔
194	+ "radvel.Parameter objects as input.")
195
196	def __repr__(self):	1✔
197	length = self.hparams['gp_length'].value	1✔
198	amp = self.hparams['gp_amp'].value	1✔
199	per = self.hparams['gp_per'].value	1✔
200	return "Per Kernel with length: {}, amp: {}, per: {}".format(	1✔
201	length, amp, per
202	)
203
204	def compute_distances(self, x1, x2):	1✔
205	X1 = np.array([x1]).T	1✔
206	X2 = np.array([x2]).T	1✔
207	self.dist = scipy.spatial.distance.cdist(X1, X2, 'euclidean')	1✔
208
209	def compute_covmatrix(self, errors):	1✔
210	""" Compute the covariance matrix, and optionally add errors along
211	the diagonal.
212
213	Args:
214	errors (float or numpy array): If covariance matrix is non-square,
215	this arg must be set to 0. If covariance matrix is square,
216	this can be a numpy array of observational errors and jitter
217	added in quadrature.
218	"""
219	length= self.hparams['gp_length'].value	1✔
220	amp = self.hparams['gp_amp'].value	1✔
221	per = self.hparams['gp_per'].value	1✔
222
223	K = amp*2 np.exp(-np.sin(np.piself.dist/per)2. / (2.length**2))	1✔
224	self.covmatrix = K	1✔
225	# add errors along the diagonal
226	try:	1✔
227	self.covmatrix += (errors*2) np.identity(K.shape[0])	1✔
228	except ValueError: # errors can't be added along diagonal to a non-square array	×
229	pass	×
230
231	return self.covmatrix	1✔
232
233	class QuasiPerKernel(Kernel):	1✔
234	"""
235	Class that computes and stores a quasi periodic kernel matrix.
236	An arbitrary element, :math:`C_{ij}`, of the matrix is:
237
238	.. math::
239
240	C_{ij} = \\eta_1^2 * exp( \\frac{ -\|t_i - t_j\|^2 }{ \\eta_2^2 } -
241	\\frac{ \\sin^2(\\frac{ \\pi\|t_i-t_j\| }{ \\eta_3 } ) }{ 2\\eta_4^2 } )
242
243	Args:
244	hparams (dict of radvel.Parameter): dictionary containing
245	radvel.Parameter objects that are GP hyperparameters
246	of this kernel. Must contain exactly four objects, 'gp_explength*',
247	'gp_amp', 'gp_per', and 'gp_perlength', where is a suffix
248	identifying these hyperparameters with a likelihood object.
249
250	"""
251	@property	1✔
252	def name(self):	1✔
253	return "QuasiPer"	×
254
255	def __init__(self, hparams):	1✔
256	self.covmatrix = None	1✔
257	self.hparams = {}	1✔
258	for par in hparams:	1✔
259	if par.startswith('gp_perlength'):	1✔
260	self.hparams['gp_perlength'] = hparams[par]	1✔
261	if par.startswith('gp_amp'):	1✔
262	self.hparams['gp_amp'] = hparams[par]	1✔
263	if par.startswith('gp_per') and not 'length' in par:	1✔
264	self.hparams['gp_per'] = hparams[par]	1✔
265	if par.startswith('gp_explength'):	1✔
266	self.hparams['gp_explength'] = hparams[par]	1✔
267
268	assert len(hparams) == 4, \	1✔
269	"QuasiPerKernel requires exactly 4 hyperparameters with names" \
270	+ " 'gp_perlength', 'gp_amp', 'gp_per', and 'gp_explength'."
271
272	try:	1✔
273	self.hparams['gp_perlength'].value	1✔
274	self.hparams['gp_amp'].value	1✔
275	self.hparams['gp_per'].value	1✔
276	self.hparams['gp_explength'].value	1✔
277	except KeyError:	1✔
278	raise KeyError("QuasiPerKernel requires hyperparameters" \	×
279	+ " 'gp_perlength', 'gp_amp', 'gp_per*', " \
280	+ "and 'gp_explength*'.")
281	except AttributeError:	1✔
282	raise AttributeError("QuasiPerKernel requires dictionary of" \	1✔
283	+ " radvel.Parameter objects as input.")
284
285	def __repr__(self):	1✔
286	perlength = self.hparams['gp_perlength'].value	1✔
287	amp = self.hparams['gp_amp'].value	1✔
288	per = self.hparams['gp_per'].value	1✔
289	explength = self.hparams['gp_explength'].value	1✔
290
291	msg = (	1✔
292	"QuasiPer Kernel with amp: {}, per length: {}, per: {}, "
293	"exp length: {}"
294	).format(amp, perlength, per, explength)
295	return msg	1✔
296
297	def compute_distances(self, x1, x2):	1✔
298	X1 = np.array([x1]).T	1✔
299	X2 = np.array([x2]).T	1✔
300	self.dist_p = scipy.spatial.distance.cdist(X1, X2, 'euclidean')	1✔
301	self.dist_se = scipy.spatial.distance.cdist(X1, X2, 'sqeuclidean')	1✔
302
303	def compute_covmatrix(self, errors):	1✔
304	""" Compute the covariance matrix, and optionally add errors along
305	the diagonal.
306
307	Args:
308	errors (float or numpy array): If covariance matrix is non-square,
309	this arg must be set to 0. If covariance matrix is square,
310	this can be a numpy array of observational errors and jitter
311	added in quadrature.
312	"""
313	perlength = self.hparams['gp_perlength'].value	1✔
314	amp = self.hparams['gp_amp'].value	1✔
315	per = self.hparams['gp_per'].value	1✔
316	explength = self.hparams['gp_explength'].value	1✔
317
318	# Check for problematic parameter values
319	# Allow 0 values for model comparison (when GP is disabled)
320	if not np.isfinite(amp) or amp < 0:	1✔
NEW 321	raise ValueError(f"Invalid gp_amp value: {amp}")	×
322	if not np.isfinite(per) or per < 0:	1✔
NEW 323	raise ValueError(f"Invalid gp_per value: {per}")	×
324	if not np.isfinite(explength) or explength < 0:	1✔
NEW 325	raise ValueError(f"Invalid gp_explength value: {explength}")	×
326	if not np.isfinite(perlength) or perlength < 0:	1✔
NEW 327	raise ValueError(f"Invalid gp_perlength value: {perlength}")	×
328
329	# If any GP parameter is 0, this is likely a "no GP" model comparison
330	# In this case, return a zero covariance matrix
331	if amp == 0 or per == 0 or explength == 0 or perlength == 0:	1✔
NEW 332	n = len(self.dist_se)	×
NEW 333	K = np.zeros((n, n))	×
NEW 334	if isinstance(errors, (int, float)) and errors == 0:	×
NEW 335	return K	×
336	else:
NEW 337	K[np.diag_indices_from(K)] += errors**2	×
NEW 338	return K	×
339
340	# Compute covariance matrix step by step
341	exp_term1 = -self.dist_se/(explength**2)	1✔
342	exp_term2 = (-np.sin(np.piself.dist_p/per)2.) / (2.perlength**2)	1✔
343
344	# Clamp extreme values to prevent overflow
345	exp_term1 = np.clip(exp_term1, -100, 100)	1✔
346	exp_term2 = np.clip(exp_term2, -100, 100)	1✔
347
348	# Compute exponential terms
349	exp1 = np.exp(exp_term1)	1✔
350	exp2 = np.exp(exp_term2)	1✔
351
352	# Compute final covariance matrix
353	K = np.array(amp*2 exp1 * exp2)	1✔
354
355	self.covmatrix = K	1✔
356
357	# add errors along the diagonal
358	try:	1✔
359	self.covmatrix += (errors*2) np.identity(K.shape[0])	1✔
360	except ValueError: # errors can't be added along diagonal to a non-square array	1✔
361	pass	1✔
362
363	return self.covmatrix	1✔
364
365
366	class CeleriteKernel(Kernel):	1✔
367	"""
368	Class that computes and stores a matrix approximating the quasi-periodic
369	kernel.
370
371	See `radvel/example_planets/k2-131_celerite.py` for an example of a setup
372	file that uses this Kernel object.
373
374	See celerite.readthedocs.io and Foreman-Mackey et al. 2017. AJ, 154, 220
375	(equation 56) for more details.
376
377	An arbitrary element, :math:`C_{ij}`, of the matrix is:
378
379	.. math::
380
381	C_{ij} = B/(2+C) * exp( -\|t_i - t_j\| / L) * (\\cos(\\frac{ 2\\pi\|t_i-t_j\| }{ P_{rot} }) + (1+C) )
382
383	Args:
384	hparams (dict of radvel.Parameter): dictionary containing
385	radvel.Parameter objects that are GP hyperparameters
386	of this kernel. Must contain exactly four objects, 'gp_B*',
387	'gp_C', 'gp_L', and 'gp_Prot', where is a suffix
388	identifying these hyperparameters with a likelihood object.
389	"""
390
391	@property	1✔
392	def name(self):	1✔
393	return "Celerite"	×
394
395	def __init__(self, hparams):	1✔
396
397	self.hparams = {}	1✔
398	for par in hparams:	1✔
399	if par.startswith('gp_B'):	1✔
400	self.hparams['gp_B'] = hparams[par]	1✔
401	if par.startswith('gp_C'):	1✔
402	self.hparams['gp_C'] = hparams[par]	1✔
403	if par.startswith('gp_L'):	1✔
404	self.hparams['gp_L'] = hparams[par]	1✔
405	if par.startswith('gp_Prot'):	1✔
406	self.hparams['gp_Prot'] = hparams[par]	1✔
407
408	assert len(self.hparams) == 4, """	1✔
409	CeleriteKernel requires exactly 4 hyperparameters with names 'gp_B', 'gp_C', 'gp_L', and 'gp_Prot'.
410	"""
411
412	try:	1✔
413	self.hparams['gp_Prot'].value	1✔
414	self.hparams['gp_C'].value	1✔
415	self.hparams['gp_B'].value	1✔
416	self.hparams['gp_L'].value	1✔
417	except KeyError:	1✔
418	raise KeyError("""	×
419	CeleriteKernel requires hyperparameters 'gp_B', 'gp_C', 'gp_L', and 'gp_Prot*'.
420	""")
421	except AttributeError:	1✔
422	raise AttributeError("CeleriteKernel requires dictionary of radvel.Parameter objects as input.")	1✔
423
424	# get arrays of real and complex parameters
425	def compute_real_and_complex_hparams(self):	1✔
426
427	self.real = np.zeros((1, 4))	1✔
428	self.complex = np.zeros((1, 4))	1✔
429
430	B = self.hparams['gp_B'].value	1✔
431	C = self.hparams['gp_C'].value	1✔
432	L = self.hparams['gp_L'].value	1✔
433	Prot = self.hparams['gp_Prot'].value	1✔
434
435	# Foreman-Mackey et al. (2017) eq 56
436	self.real[0,0] = B*(1+C)/(2+C)	1✔
437	self.real[0,2] = 1/L	1✔
438	self.complex[0,0] = B/(2+C)	1✔
439	self.complex[0,1] = 0.	1✔
440	self.complex[0,2] = 1/L	1✔
441	self.complex[0,3] = 2*np.pi/Prot	1✔
442
443	def __repr__(self):	1✔
444
445	B = self.hparams['gp_B'].value	1✔
446	C = self.hparams['gp_C'].value	1✔
447	L = self.hparams['gp_L'].value	1✔
448	Prot = self.hparams['gp_Prot'].value	1✔
449
450	msg = (	1✔
451	"Celerite Kernel with B = {}, C = {}, L = {}, Prot = {}."
452	).format(B, C, L, Prot)
453	return msg	1✔
454
455	def compute_distances(self, x1, x2):	1✔
456	"""
457	The celerite.solver.CholeskySolver object does
458	not require distances to be precomputed, so
459	this method has been co-opted to define some
460	unchanging variables.
461	"""
462	self.x = x1	1✔
463
464	# blank matrices (corresponding to Cholesky decomp of kernel) needed for celerite solver
465	self.A = np.empty(0)	1✔
466	self.U = np.empty((0,0))	1✔
467	self.V = self.U	1✔
468
469
470	def compute_covmatrix(self, errors):	1✔
471	""" Compute the Cholesky decomposition of a celerite kernel
472
473	Args:
474	errors (array of float): observation errors and jitter added
475	in quadrature
476
477	Returns:
478	celerite.solver.CholeskySolver: the celerite solver object,
479	with Cholesky decomposition computed.
480	"""
481	# initialize celerite solver object
482	solver = CholeskySolver()	1✔
483
484	self.compute_real_and_complex_hparams()	1✔
485	solver.compute(	1✔
486	0., self.real[:,0], self.real[:,2],
487	self.complex[:,0], self.complex[:,1],
488	self.complex[:,2], self.complex[:,3],
489	self.A, self.U, self.V,
490	self.x, errors**2
491	)
492
493	return solver	1✔

California-Planet-Search / radvel / 17445970354

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous