Open Software for Astrophysics, AAS241

Slide 1

Slide 1 text

OPEN SOFTWARE FOR ASTROPHYSICS Dan Foreman-Mackey

Slide 2

Slide 2 text

No content

Slide 3

Slide 3 text

case study: Gaussian Processes

Slide 4

Slide 4 text

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

Slide 5

Slide 5 text

°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)

Slide 6

Slide 6 text

°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)

Slide 7

Slide 7 text

reference: Aigrain & DFM (2022)

Slide 8

Slide 8 text

reference: Aigrain & DFM (2022)

Slide 9

Slide 9 text

reference: Aigrain & DFM (2022) ignoring correlated noise accounting for correlated noise

Slide 10

Slide 10 text

reference: Aigrain & DFM (2022)

Slide 11

Slide 11 text

a Gaussian Process is a drop - in replacement for chi - squared

Slide 12

Slide 12 text

more details: Aigrain & Foreman-Mackey (2023) arXiv:2209.08940

Slide 13

Slide 13 text

7 [1] model building [2] computational cost

Slide 14

Slide 14 text

k(tn , tm ; θ) “kernel” or “covariance”

Slide 15

Slide 15 text

No content

Slide 16

Slide 16 text

import george import celerite import tinygp

Slide 17

Slide 17 text

my f i rst try: george 1

Slide 18

Slide 18 text

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)

Slide 19

Slide 19 text

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)

Slide 20

Slide 20 text

k(tn , tm ; θ) “kernel” or “covariance”

Slide 21

Slide 21 text

from george.kernels import * k1 = 1.5 * ExpSquaredKernel(2.3) k2 = 5.5 * Matern32Kernel(0.1) kernel = 0.5 * (k1 + k2)

Slide 22

Slide 22 text

from george import GP gp = GP(kernel) gp.compute(x, yerr) gp.log_likelihood(y)

Slide 23

Slide 23 text

from george import GP gp = GP(kernel) gp.compute(x, yerr) gp.log_likelihood(y) gp.f i t(y) ???

Slide 24

Slide 24 text

the astronomical Python ecosystem + MANY MORE!

Slide 25

Slide 25 text

* API design (library vs scripts) * don’t reinvent the wheel

Slide 26

Slide 26 text

No content

Slide 27

Slide 27 text

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

Slide 28

Slide 28 text

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)

Slide 29

Slide 29 text

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)

Slide 30

Slide 30 text

No content

Slide 31

Slide 31 text

“semi/quasi - separable” matrices

Slide 32

Slide 32 text

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)

Slide 33

Slide 33 text

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)

Slide 34

Slide 34 text

No content

Slide 35

Slide 35 text

Slide 36

Slide 36 text

+ + vs

Slide 37

Slide 37 text

* interdisciplinary collaboration * importance of implementation

Slide 38

Slide 38 text

7 [1] 1 (ish) dimensional input [2] specif i c type of kernel restrictions:

Slide 39

Slide 39 text

modern infrastructure: tinygp 3

Slide 40

Slide 40 text

what’s missing from the astronomical Python ecosystem?

Slide 41

Slide 41 text

7 [1] differentiable programming [2] hardware acceleration

Slide 42

Slide 42 text

the broader numerical computing Python ecosystem + SO MANY MORE!

Slide 43

Slide 43 text

jax.readthedocs.io

Slide 44

Slide 44 text

import numpy as np def linear_least_squares(x, y) : A = np.vander(x, 2) return np.linalg.lstsq(A, y)[0]

Slide 45

Slide 45 text

import jax.numpy as jnp def linear_least_squares(x, y) : A = jnp.vander(x, 2) return jnp.linalg.lstsq(A, y)[0]

Slide 46

Slide 46 text

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]

Slide 47

Slide 47 text

No content

Slide 48

Slide 48 text

tinygp.readthedocs.io

Slide 49

Slide 49 text

the broader numerical computing Python ecosystem + SO MANY MORE!

Slide 50

Slide 50 text

* I <3 JAX * don’t reinvent the wheel

Slide 51

Slide 51 text

the why & how of open software in astrophysics

Slide 52

Slide 52 text

credit: Adrian Price-Whelan / / data: SAO/NASA ADS

Slide 53

Slide 53 text

No content

Slide 54

Slide 54 text

No content

Slide 55

Slide 55 text

No content

Slide 56

Slide 56 text

No content

Slide 57

Slide 57 text

takeaways

Slide 58

Slide 58 text

open software is foundational to astrophysics research let’s consider & discuss interface design and user interaction leverage existing infrastructure & learn when to start fresh