Skip to content

Commit 46132ea

Browse files
hpparvihpparvi
committed
Merge branch 'v2.2' into dev
# Conflicts: # pytransit/models/numba/ma_quadratic_nb.py # pytransit/models/numba/ma_uniform_nb.py # pytransit/version.py # tests/test_qpower2_nb.py
2 parents 9c6f61c + f1cb530 commit 46132ea

17 files changed

Lines changed: 890 additions & 507 deletions

‎CHANGELOG.md‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,11 @@
11
# Changelog
22

3+
## [2.2.0] - 2020-09-13
4+
5+
PyTransit version 2.2 now calculates the normalized planet-star distances using a Taylor series expansion
6+
of the x, and y positions in the sky plane (Parviainen and Korth, 2020, submitted to MNRAS). This gives a
7+
significant speed boost in transit model evaluation that is especially noticeable for eccentric orbits.
8+
39
## [2.1.0] - 2020-07-07
410

511
PyTransit version 2.1 adds a new transit model named *swift* that can use any Python callable to model the stellar

‎README.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -173,7 +173,7 @@ period `p`, scaled semi-major axis `a`, orbital inclination `i`, eccentricity `e
173173
an array of limb darkening coefficients `ldc`
174174

175175
```Python
176-
flux = tm.evaluate_ps(k, ldc, t0, p, a, i, e, w)
176+
flux = tm.evaluate_ps(k,ldc,t0,p,a,i,e,w)
177177
```
178178

179179
or using either a parameter array

‎pytransit/lpf/loglikelihood/multiceleriteloglikelihood.py‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
# PyTransit: fast and easy exoplanet transit modelling in Python.
2+
# Copyright (C) 2010-2020 Hannu Parviainen
3+
#
4+
# This program is free software: you can redistribute it and/or modify
5+
# it under the terms of the GNU General Public License as published by
6+
# the Free Software Foundation, either version 3 of the License, or
7+
# (at your option) any later version.
8+
#
9+
# This program is distributed in the hope that it will be useful,
10+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12+
# GNU General Public License for more details.
13+
#
14+
# You should have received a copy of the GNU General Public License
15+
# along with this program. If not, see <https://www.gnu.org/licenses/>.
16+
#
17+
# This program is free software: you can redistribute it and/or modify
18+
# it under the terms of the GNU General Public License as published by
19+
# the Free Software Foundation, either version 3 of the License, or
20+
# (at your option) any later version.
21+
#
22+
# This program is distributed in the hope that it will be useful,
23+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
24+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25+
# GNU General Public License for more details.
26+
#
27+
28+
from numpy import asarray, unique, zeros, inf, squeeze, zeros_like, isfinite, log10, diff, sqrt
29+
from astropy.stats import mad_std
30+
31+
try:
32+
from celerite import GP
33+
from celerite.terms import Matern32Term
34+
with_celerite = True
35+
except ImportError:
36+
with_celerite = False
37+
38+
from ...param import LParameter, NormalPrior as NP, UniformPrior as UP
39+
40+
41+
class MultiCeleriteLogLikelihood:
42+
def __init__(self, lpf, name: str = 'gp', lcids=None, fixed_hps=None):
43+
if not with_celerite:
44+
raise ImportError("MultiCeleriteLogLikelihood requires celerite.")
45+
46+
self.name = name
47+
self.lpf = lpf
48+
49+
if fixed_hps is None:
50+
self.free = True
51+
else:
52+
self.hps = asarray(fixed_hps)
53+
self.free = False
54+
55+
if lpf.lcids is None:
56+
raise ValueError('The LPF data needs to be initialised before initialising CeleriteLogLikelihood.')
57+
self.lcids = lcids if lcids is not None else unique(lpf.lcids)
58+
self.lcslices = [lpf.lcslices[i] for i in self.lcids]
59+
self.nlc = len(self.lcslices)
60+
61+
self.times, self.fluxes, self.wns = [], [], []
62+
for i in self.lcids:
63+
self.times.append(lpf.times[i])
64+
self.fluxes.append(lpf.fluxes[i])
65+
self.wns.append(diff(lpf.fluxes[i]).std()/sqrt(2))
66+
67+
self.gps = [GP(Matern32Term(0, 0)) for i in range(self.nlc)]
68+
69+
if self.free:
70+
self.init_parameters()
71+
else:
72+
self.compute_gp(None, force=True, hps=self.hps)
73+
74+
def init_parameters(self):
75+
name = self.name
76+
pgp = [LParameter(f'{name}_ln_out', f'{name} ln output scale', '', NP(-6.5, 0.5), bounds=(-inf, inf)),
77+
LParameter(f'{name}_ln_in', f'{name} ln input scale', '', UP(-8, 8), bounds=(-inf, inf))]
78+
self.lpf.ps.thaw()
79+
self.lpf.ps.add_global_block(self.name, pgp)
80+
self.lpf.ps.freeze()
81+
self.pv_slice = self.lpf.ps.blocks[-1].slice
82+
self.pv_start = self.lpf.ps.blocks[-1].start
83+
setattr(self.lpf, f"_sl_{name}", self.pv_slice)
84+
setattr(self.lpf, f"_start_{name}", self.pv_start)
85+
86+
def compute_gp(self, pv, force: bool = False, hps=None):
87+
if self.free or force:
88+
parameters = pv[self.pv_slice] if hps is None else hps
89+
for i, gp in enumerate(self.gps):
90+
gp.set_parameter_vector(parameters)
91+
gp.compute(self.times[i], yerr=self.wns[i])
92+
93+
def compute_gp_lnlikelihood(self, pv, model):
94+
self.compute_gp(pv)
95+
lnlike = 0.0
96+
for i, gp in enumerate(self.gps):
97+
lnlike += gp.log_likelihood(self.fluxes[i] - model[self.lcslices[i]])
98+
return lnlike
99+
100+
def predict_baseline(self, pv):
101+
self.compute_gp(pv)
102+
residuals = self.lpf.ofluxa - squeeze(self.lpf.transit_model(pv))
103+
bl = zeros_like(self.lpf.timea)
104+
for i, gp in enumerate(self.gps):
105+
sl = self.lcslices[i]
106+
bl[sl] = gp.predict(residuals[sl], self.times[i], return_cov=False)
107+
return 1. + bl
108+
109+
def __call__(self, pvp, model):
110+
if pvp.ndim == 1:
111+
lnlike = self.compute_gp_lnlikelihood(pvp, model)
112+
else:
113+
lnlike = zeros(pvp.shape[0])
114+
for ipv, pv in enumerate(pvp):
115+
if all(isfinite(model[ipv])):
116+
lnlike[ipv] = self.compute_gp_lnlikelihood(pv, model[ipv])
117+
else:
118+
lnlike[ipv] = -inf
119+
return lnlike

‎pytransit/models/ma_quadratic.py‎

Lines changed: 18 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -28,11 +28,9 @@
2828
# along with this program. If not, see <https://www.gnu.org/licenses/>.
2929
from typing import Union, Optional, List
3030

31-
from numpy import ndarray, array, squeeze, atleast_2d, atleast_1d, zeros, asarray
31+
from numpy import ndarray, array, squeeze, atleast_2d, atleast_1d, zeros, asarray, inf
3232

33-
from .numba.ma_quadratic_nb import quadratic_model_interpolated_pv, calculate_interpolation_tables, \
34-
quadratic_model_interpolated_v, quadratic_model_interpolated_s, quadratic_model_direct_v, quadratic_model_direct_s, \
35-
quadratic_model_direct_pv
33+
from .numba.ma_quadratic_nb import quadratic_model_pv, calculate_interpolation_tables, quadratic_model_v, quadratic_model_s
3634
from .transitmodel import TransitModel
3735

3836
__all__ = ['QuadraticModel']
@@ -41,7 +39,7 @@ class QuadraticModel(TransitModel):
4139
"""Transit model with quadratic limb darkening (Mandel & Agol, ApJ 580, L171-L175 2002).
4240
"""
4341

44-
def __init__(self, interpolate: bool = True, klims: tuple = (0.005, 0.5), nk: int = 256, nz: int = 256):
42+
def __init__(self, interpolate: bool = True, klims: tuple = (0.01, 0.5), nk: int = 256, nz: int = 256):
4543
"""Transit model with quadratic limb darkening (Mandel & Agol, ApJ 580, L171-L175 2002).
4644
4745
Parameters
@@ -61,13 +59,17 @@ def __init__(self, interpolate: bool = True, klims: tuple = (0.005, 0.5), nk: in
6159
# Interpolation tables for the model components
6260
# ---------------------------------------------
6361
if interpolate:
62+
self._interpolation_initialised = True
6463
self.ed, self.le, self.ld, self.kt, self.zt = calculate_interpolation_tables(klims[0], klims[1], nk, nz)
6564
self.klims = klims
6665
self.nk = nk
6766
self.nz = nz
6867
else:
69-
self.ed, self.le, self.ld, self.kt, self.zt = None, None, None, None, None
70-
self.klims, self.nk, self.nz = None, None, None
68+
self._interpolation_initialised = False
69+
self.ed, self.le, self.ld, self.kt, self.zt = zeros((2,2)), zeros((2,2)), zeros((2,2)), zeros(2), zeros(2)
70+
self.klims = klims
71+
self.nk = 0
72+
self.nz = 0
7173

7274
def evaluate(self, k: Union[float, ndarray], ldc: Union[ndarray, List], t0: Union[float, ndarray], p: Union[float, ndarray],
7375
a: Union[float, ndarray], i: Union[float, ndarray], e: Optional[Union[float, ndarray]] = None,
@@ -125,15 +127,9 @@ def evaluate(self, k: Union[float, ndarray], ldc: Union[ndarray, List], t0: Unio
125127
e = zeros(npv) if e is None else e
126128
w = zeros(npv) if w is None else w
127129

128-
if self.interpolate:
129-
flux = quadratic_model_interpolated_v(self.time, k, t0, p, a, i, e, w, ldc,
130-
self.lcids, self.pbids, self.nsamples, self.exptimes, self.npb,
131-
self._es, self._ms, self._tae, self.ed, self.ld, self.le,
132-
self.kt, self.zt)
133-
else:
134-
flux = quadratic_model_direct_v(self.time, k, t0, p, a, i, e, w, ldc,
135-
self.lcids, self.pbids, self.nsamples, self.exptimes, self.npb,
136-
self._es, self._ms, self._tae)
130+
flux = quadratic_model_v(self.time, k, t0, p, a, i, e, w, ldc, self.lcids, self.pbids, self.nsamples, self.exptimes, self.npb,
131+
self.ed, self.ld, self.le, self.kt, self.zt, self.interpolate)
132+
137133
return squeeze(flux)
138134

139135
def evaluate_ps(self, k: Union[float, ndarray], ldc: ndarray, t0: float, p: float, a: float, i: float,
@@ -153,7 +149,7 @@ def evaluate_ps(self, k: Union[float, ndarray], ldc: ndarray, t0: float, p: floa
153149
a : float
154150
Orbital semi-major axis divided by the stellar radius as a float.
155151
i : float
156-
Orbital inclination(s) as a float.
152+
Orbital inclination as a float.
157153
e : float, optional
158154
Orbital eccentricity as a float.
159155
w : float, optional
@@ -179,15 +175,8 @@ def evaluate_ps(self, k: Union[float, ndarray], ldc: ndarray, t0: float, p: floa
179175
if ldc.size != 2 * self.npb:
180176
raise ValueError("The quadratic model needs two limb darkening coefficients per passband")
181177

182-
if self.interpolate:
183-
flux = quadratic_model_interpolated_s(self.time, k, t0, p, a, i, e, w, ldc,
184-
self.lcids, self.pbids, self.nsamples, self.exptimes, self.npb,
185-
self._es, self._ms, self._tae, self.ed, self.ld, self.le,
186-
self.kt, self.zt)
187-
else:
188-
flux = quadratic_model_direct_s(self.time, k, t0, p, a, i, e, w, ldc,
189-
self.lcids, self.pbids, self.nsamples, self.exptimes, self.npb,
190-
self._es, self._ms, self._tae)
178+
flux = quadratic_model_s(self.time, k, t0, p, a, i, e, w, ldc, self.lcids, self.pbids, self.nsamples, self.exptimes, self.npb,
179+
self.ed, self.ld, self.le, self.kt, self.zt, self.interpolate)
191180
return squeeze(flux)
192181

193182
def evaluate_pv(self, pvp: ndarray, ldc: ndarray, copy: bool = True) -> ndarray:
@@ -198,7 +187,7 @@ def evaluate_pv(self, pvp: ndarray, ldc: ndarray, copy: bool = True) -> ndarray:
198187
pvp: ndarray
199188
Parameter array with a shape `(npv, npar)` where `npv` is the number of parameter vectors, and each row
200189
contains a set of parameters `[k, t0, p, a, i, e, w]`. The radius ratios can also be given per passband,
201-
in which case the row should be structured as `[k_0, k_1, k_2, ..., k_npb, t0, p, a, i, e, w]`.
190+
in which case the row should be structured as `[k_0, k_1, k_2, ..., k_npb, t0, p, a, b, e, w]`.
202191
ldc: ndarray
203192
Limb darkening coefficient array with shape `(npv, 2*npb)`, where `npv` is the number of parameter vectors
204193
and `npb` is the number of passbands.
@@ -220,13 +209,8 @@ def evaluate_pv(self, pvp: ndarray, ldc: ndarray, copy: bool = True) -> ndarray:
220209
if self.time is None:
221210
raise ValueError("Need to set the data before calling the transit model.")
222211

223-
if self.interpolate:
224-
flux = quadratic_model_interpolated_pv(self.time, pvp, ldc, self.lcids, self.pbids, self.nsamples, self.exptimes,
225-
self.npb, self._es, self._ms, self._tae, self.ed, self.ld, self.le,
226-
self.kt, self.zt)
227-
else:
228-
flux = quadratic_model_direct_pv(self.time, pvp, ldc, self.lcids, self.pbids, self.nsamples, self.exptimes,
229-
self.npb, self._es, self._ms, self._tae)
212+
flux = quadratic_model_pv(self.time, pvp, ldc, self.lcids, self.pbids, self.nsamples, self.exptimes,
213+
self.npb, self.ed, self.ld, self.le, self.kt, self.zt, self.interpolate)
230214
return squeeze(flux)
231215

232216
def to_opencl(self):

‎pytransit/models/ma_uniform.py‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ def evaluate(self, k: npfloat, t0: npfloat, p: npfloat, a: npfloat, i: npfloat,
7070
e = 0. if e is None else e
7171
w = 0. if w is None else w
7272
flux = uniform_model_s(self.time, k, t0, p, a, i, e, w, self.lcids, self.pbids, self.nsamples,
73-
self.exptimes, self._es, self._ms, self._tae, zsign=self._zsign)
73+
self.exptimes, zsign=self._zsign)
7474

7575
# Parameter population branch
7676
# ---------------------------
@@ -82,8 +82,7 @@ def evaluate(self, k: npfloat, t0: npfloat, p: npfloat, a: npfloat, i: npfloat,
8282
if k.ndim == 1:
8383
k = k.reshape((k.size,1))
8484

85-
flux = uniform_model_v(self.time, k, t0, p, a, i, e, w, self.lcids, self.pbids, self.nsamples,
86-
self.exptimes, self._es, self._ms, self._tae, zsign=self._zsign)
85+
flux = uniform_model_v(self.time, k, t0, p, a, i, e, w, self.lcids, self.pbids, self.nsamples, self.exptimes, zsign=self._zsign)
8786
return squeeze(flux)
8887

8988
def evaluate_ps(self, k: float, t0: float, p: float, a: float, i: float, e: float = 0., w: float = 0.) -> ndarray:
@@ -121,8 +120,7 @@ def evaluate_ps(self, k: float, t0: float, p: float, a: float, i: float, e: floa
121120
raise ValueError("Need to set the data before calling the transit model.")
122121

123122
k = asarray(k)
124-
flux = uniform_model_s(self.time, k, t0, p, a, i, e, w, self.lcids, self.pbids, self.nsamples, self.exptimes,
125-
self._es, self._ms, self._tae, zsign=self._zsign)
123+
flux = uniform_model_s(self.time, k, t0, p, a, i, e, w, self.lcids, self.pbids, self.nsamples, self.exptimes, zsign=self._zsign)
126124
return squeeze(flux)
127125

128126
def evaluate_pv(self, pvp: ndarray) -> ndarray:
@@ -146,6 +144,5 @@ def evaluate_pv(self, pvp: ndarray) -> ndarray:
146144
Modelled flux either as a 1D or 2D ndarray.
147145
"""
148146
assert self.time is not None, "Need to set the data before calling the transit model."
149-
flux = uniform_model_pv(self.time, pvp, self.lcids, self.pbids, self.nsamples, self.exptimes,
150-
self._es, self._ms, self._tae, zsign=self._zsign)
147+
flux = uniform_model_pv(self.time, pvp, self.lcids, self.pbids, self.nsamples, self.exptimes, zsign=self._zsign)
151148
return squeeze(flux)

0 commit comments

Comments
 (0)