Dan Foreman-Mackey
January 12, 2023
310

# Open Software for Astrophysics, AAS241

Slides for my plenary talk at the 241st American Astronomical Society meeting.

January 12, 2023

## Transcript

1. OPEN

SOFTWARE

FOR

ASTROPHYSICS
Dan Foreman-Mackey

2. case study:

Gaussian Processes

3. AAS 225 / 2015 / Seattle AAS 231 / 2018 / National Harbor

4. °0.6
°0.3
0.0
0.3
0.6
raw [ppt]
0 5 10 15 20 25
time [days]
°0.30
°0.15
0.00
de-trended [ppt]
N = 1000
reference: DFM+ (2017)

5. °0.6
°0.3
0.0
0.3
0.6
raw [ppt]
0 5 10 15 20 25
time [days]
°0.30
°0.15
0.00
de-trended [ppt]
N = 1000
reference: DFM+ (2017)

6. reference: Aigrain & DFM (2022)

7. reference: Aigrain & DFM (2022)

8. reference: Aigrain & DFM (2022)
ignoring

correlated

noise
accounting

for

correlated

noise

9. reference: Aigrain & DFM (2022)

10. a Gaussian Process is a

drop
-
in replacement for

chi
-
squared

11. more details:

Aigrain & Foreman-Mackey (2023)

arXiv:2209.08940

12. 7
[1] model building

[2] computational cost

13. k(tn
, tm
; θ)

“kernel” or “covariance”

14. import george

import celerite

import tinygp

15. my f
i
rst try: george
1

16. import numpy as np

def log_likelihood(params, x, diag, r)
:

K = build_kernel_matrix(params, x, diag)

gof = r.T @ np.linalg.solve(K, r)

norm = np.linalg.slogdet(K)[1]

return -0.5 * (gof + norm)

17. import numpy as np

def log_likelihood(params, x, diag, r)
:

K = build_kernel_matrix(params, x, diag)

gof = r.T @ np.linalg.solve(K, r)

norm = np.linalg.slogdet(K)[1]

return -0.5 * (gof + norm)

18. k(tn
, tm
; θ)

“kernel” or “covariance”

19. from george.kernels import *

k1 = 1.5 * ExpSquaredKernel(2.3)

k2 = 5.5 * Matern32Kernel(0.1)

kernel = 0.5 * (k1 + k2)

20. from george import GP

gp = GP(kernel)

gp.compute(x, yerr)

gp.log_likelihood(y)

21. from george import GP

gp = GP(kernel)

gp.compute(x, yerr)

gp.log_likelihood(y)

gp.f
i
t(y) ???

22. the astronomical

Python ecosystem
+ MANY MORE!

23. * API design (library vs scripts)

* don’t reinvent the wheel

24. faster: celerite*
2
* yes, that truly is how you pronounce it…

25. import numpy as np

def log_likelihood(params, x, diag, r)
:

K = build_kernel_matrix(params, x, diag)

gof = r.T @ np.linalg.solve(K, r)

norm = np.linalg.slogdet(K)[1]

return -0.5 * (gof + norm)

26. import numpy as np

def log_likelihood(params, x, diag, r)
:

K = build_kernel_matrix(params, x, diag)

gof = r.T @ np.linalg.solve(K, r)

norm = np.linalg.slogdet(K)[1]

return -0.5 * (gof + norm)

27. “semi/quasi
-
separable” matrices

28. 102 103 104 105
number of data points [N]
10 5
10 4
10 3
10 2
10 1
100
computational cost [seconds]
1
2
4
8
16
32
64
128
256
direct
O(N)
100 101
number o
reference: DFM, Agol, Ambikasaran, Angus (2017)

29. 102 103 104 105
number of data points [N]
10 4
10 3
10 2
10 1
100
computational cost [seconds]
1
2
4
8
16
32
64
128
256
O(N)
100 101
number o
reference: DFM, Agol, Ambikasaran, Angus (2017)

30. +

31. + +
vs

32. * interdisciplinary collaboration

* importance of implementation

33. 7
[1] 1 (ish) dimensional input

[2] specif
i
c type of kernel
restrictions:

34. modern infrastructure: tinygp
3

35. what’s missing from

the astronomical

Python ecosystem?

36. 7
[1] differentiable programming

[2] hardware acceleration

numerical computing

Python ecosystem
+ SO MANY MORE!

39. import numpy as np

def linear_least_squares(x, y)
:

A = np.vander(x, 2)

return np.linalg.lstsq(A, y)[0]

40. import jax.numpy as jnp

def linear_least_squares(x, y)
:

A = jnp.vander(x, 2)

return jnp.linalg.lstsq(A, y)[0]

41. import jax.numpy as jnp

@jax.jit

def linear_least_squares(x, y)
:

A = jnp.vander(x, 2)

return jnp.linalg.lstsq(A, y)[0]

numerical computing

Python ecosystem
+ SO MANY MORE!

44. * I <3 JAX

* don’t reinvent the wheel

45. the why & how of

open software

in astrophysics

/ /

47. takeaways

48. open software is foundational to
astrophysics research

let’s consider & discuss interface
design and user interaction

leverage existing infrastructure &
learn when to start fresh

49. get in touch!

dfm.io

github.com/dfm