Slides for my plenary talk at the 241st American Astronomical Society meeting.
OPENSOFTWAREFORASTROPHYSICSDan Foreman-Mackey
View Slide
case study:Gaussian Processes
AAS 225 / 2015 / Seattle AAS 231 / 2018 / National Harbor
°0.6°0.30.00.30.6raw [ppt]0 5 10 15 20 25time [days]°0.30°0.150.00de-trended [ppt]N = 1000reference: DFM+ (2017)
reference: Aigrain & DFM (2022)
reference: Aigrain & DFM (2022)ignoringcorrelatednoiseaccountingforcorrelatednoise
a Gaussian Process is adrop-in replacement forchi-squared
more details:Aigrain & Foreman-Mackey (2023)arXiv:2209.08940
7[1] model building[2] computational cost
k(tn, tm; θ)“kernel” or “covariance”
import georgeimport celeriteimport tinygp
my first try: george1
import numpy as npdef 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)
from george.kernels import *k1 = 1.5 * ExpSquaredKernel(2.3)k2 = 5.5 * Matern32Kernel(0.1)kernel = 0.5 * (k1 + k2)
from george import GPgp = GP(kernel)gp.compute(x, yerr)gp.log_likelihood(y)
from george import GPgp = GP(kernel)gp.compute(x, yerr)gp.log_likelihood(y)gp.fit(y) ???
the astronomicalPython ecosystem+ MANY MORE!
* API design (library vs scripts)* don’t reinvent the wheel
faster: celerite*2* yes, that truly is how you pronounce it…
“semi/quasi-separable” matrices
102 103 104 105number of data points [N]10 510 410 310 210 1100computational cost [seconds]1248163264128256directO(N)100 101number oreference: DFM, Agol, Ambikasaran, Angus (2017)
102 103 104 105number of data points [N]10 410 310 210 1100computational cost [seconds]1248163264128256O(N)100 101number oreference: DFM, Agol, Ambikasaran, Angus (2017)
+
+ +vs
* interdisciplinary collaboration* importance of implementation
7[1] 1 (ish) dimensional input[2] specific type of kernelrestrictions:
modern infrastructure: tinygp3
what’s missing fromthe astronomicalPython ecosystem?
7[1] differentiable programming[2] hardware acceleration
the broadernumerical computingPython ecosystem+ SO MANY MORE!
jax.readthedocs.io
import numpy as npdef linear_least_squares(x, y): A = np.vander(x, 2)return np.linalg.lstsq(A, y)[0]
import jax.numpy as jnpdef linear_least_squares(x, y): A = jnp.vander(x, 2)return jnp.linalg.lstsq(A, y)[0]
import jax.numpy as jnp@jax.jitdef linear_least_squares(x, y): A = jnp.vander(x, 2)return jnp.linalg.lstsq(A, y)[0]
tinygp.readthedocs.io
* I <3 JAX* don’t reinvent the wheel
the why & how ofopen softwarein astrophysics
credit: Adrian Price-Whelan/ /data: SAO/NASA ADS
takeaways
open software is foundational toastrophysics researchlet’s consider & discuss interfacedesign and user interactionleverage existing infrastructure &learn when to start fresh
get in touch!dfm.iogithub.com/dfm