Open Software for Astrophysics, AAS241

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

Transcript

1. OPEN
   
   
   SOFTWARE
   
   
   FOR
   
   
   ASTROPHYSICS
   Dan Foreman-Mackey
   

   View Slide

2. View Slide

3. case study:
   
   
   Gaussian Processes
   

   View Slide

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

   View Slide

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)
   

   View Slide

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

   View Slide

7. reference: Aigrain & DFM (2022)
   

   View Slide

8. reference: Aigrain & DFM (2022)
   

   View Slide

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

   View Slide

10. reference: Aigrain & DFM (2022)
    

    View Slide

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

    View Slide

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

    View Slide

13. 7
    [1] model building
    
    
    [2] computational cost
    

    View Slide

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

    View Slide

15. View Slide

16. import george
    
    
    import celerite
    
    
    import tinygp
    

    View Slide

17. my f
    i
    rst try: george
    1
    

    View Slide

18. 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)
    

    View Slide

19. 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)
    

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

24. the astronomical
    
    
    Python ecosystem
    + MANY MORE!
    

    View Slide

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

    View Slide

26. View Slide

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

    View Slide

28. 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)
    

    View Slide

29. 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)
    

    View Slide

30. View Slide

31. “semi/quasi
    -
    separable” matrices
    

    View Slide

32. 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)
    

    View Slide

33. 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)
    

    View Slide

34. View Slide

35. +
    

    View Slide

36. + +
    vs
    

    View Slide

37. * interdisciplinary collaboration
    
    
    * importance of implementation
    

    View Slide

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

    View Slide

39. modern infrastructure: tinygp
    3
    

    View Slide

40. what’s missing from
    
    
    the astronomical
    
    
    Python ecosystem?
    

    View Slide

41. 7
    [1] differentiable programming
    
    
    [2] hardware acceleration
    

    View Slide

42. the broader
    
    
    numerical computing
    
    
    Python ecosystem
    + SO MANY MORE!
    

    View Slide

43. jax.readthedocs.io
    

    View Slide

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

    View Slide

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

    View Slide

46. 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]
    

    View Slide

47. View Slide

48. tinygp.readthedocs.io
    

    View Slide

49. the broader
    
    
    numerical computing
    
    
    Python ecosystem
    + SO MANY MORE!
    

    View Slide

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

    View Slide

51. the why & how of
    
    
    open software
    
    
    in astrophysics
    

    View Slide

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

    View Slide

53. View Slide

54. View Slide

55. View Slide

56. View Slide

57. takeaways
    

    View Slide

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

    View Slide

59. get in touch!
    
    
    dfm.io
    
    
    github.com/dfm
    

    View Slide

60. View Slide