BoxLeastSquares.jl

Installation

To install use Pkg. From the REPL, press ] to enter Pkg-mode

pkg> add BoxLeastSquares

If you want to use the most up-to-date version of the code, check it out from main

pkg> add BoxLeastSquares#main

Usage

First, import the package

julia> using BoxLeastSquares

you can optionally alias the package name, too

julia> import BoxLeastSquares as BLS

now, load some data. If you don't have an estimate of the y error it will default to 1.

julia> t, y, yerr = load_data(); # load data somehow

The primary interface is through the BLS method

julia> result = BLS(t, y, yerr; duration=0.16)
BLSPeriodogram
==============
input dim: 1000
output dim: 1820
period range: 0.32 - 5.014724142709022
duration range: 0.16 - 0.16
objective: likelihood

parameters
----------
index: 1633
period: 1.99930396919953
duration: 0.16
t0: 0.5001330656464655
depth: 0.19594118110109113 ± 0.0008688097746093883
snr: 225.52828804117118
log-likelihood: 27396.365214805144

to extract the parameters in a convenient named tuple use BoxLeastSquares.params

julia> BoxLeastSquares.params(result)
(index = 1633, power = 27396.365214805144, period = 1.99930396919953, duration = 0.16, t0 = 0.5001330656464655, depth = 0.19594118110109113, depth_err = 0.0008688097746093883, snr = 225.52828804117118, loglike = 27396.365214805144)

The period grid was automatically determined using autoperiod, but you can supply your own, too:

julia> periods = exp.(range(log(2) - 0.1, log(2) + 0.1, length=1000));

julia> result_fine = BLS(t, y, yerr; duration=0.12:0.01:0.20, periods=periods)
BLSPeriodogram
==============
input dim: 1000
output dim: 1000
period range: 1.809674836071919 - 2.210341836151295
duration range: 0.12 - 0.2
objective: likelihood

parameters
----------
index: 503
period: 2.001001251543549
duration: 0.168
t0: 0.4961330656464656
depth: 0.19466955969052016 ± 0.0008627202098527317
snr: 225.64622628204188
log-likelihood: 27457.6383039924

Unitful.jl

BoxLeastSquares.jl is fully compatible with Unitful.jl (although it is not a dependency of the library). For example

julia> using Unitful

julia> tu = t * u"d";

julia> results_units = BLS(tu, y, yerr; duration=(2:0.1:4)u"hr")
BLSPeriodogram
==============
input dim: 1000
output dim: 3343
period range: 0.3333333333333333 d - 4.988348864592586 d
duration range: 2.0 hr - 4.0 hr
objective: likelihood

parameters
----------
index: 2986
period: 2.0019235780121827 d
duration: 3.8000000000000003 hr
t0: 0.4916330656464656 d
depth: 0.19445716575012517 ± 0.0008692454825826517
snr: 223.70799693127577
log-likelihood: 26953.643422397385

Plotting

BoxLeastSquares.BLSPeriodogram has plotting shorthands built right in- by default it will plot the period grid and the computed power

using Plots, UnitfulRecipes

plot(results_units, label="")

now let's look at how the transit model compares to the data

pars = BoxLeastSquares.params(results_units)
wrap = 0.5 * pars.period
phases = @. (mod(t - pars.t0 + wrap, pars.period) - wrap) / pars.period
inds = sortperm(phases)
model = BoxLeastSquares.model(results_units)

scatter(phases[inds], y[inds], yerr=yerr[inds],
    label="data", xlabel="phase", xlim=(-0.2, 0.2), leg=:bottomright)
plot!(phases[inds], model[inds], lw=3, label="BLS model")

This code has been benchmarked against the C implementation in astropy.timeseries.bls. The C version uses OpenMP to multi-thread some parts of the core BLS algorithm, but BoxLeastSquares.jl has no threading support currently. For a fair comparison, we set OMP_NUM_THREADS to 1 for the following tests.

This first benchmark is simply the time it takes to evaluate the BLS periodogram. Periods are pre-computed using autoperiod. We simulate different sizes of data sets (x-axis) as well as different sizes of period grids (shape). This benchmark does not use units. The code can be found in bench/benchmark.jl. Here is the information for my system-

Julia Version 1.6.0
Commit f9720dc2eb* (2021-03-24 12:55 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin20.3.0)
  CPU: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  OMP_NUM_THREADS = 1
  JULIA_NUM_THREADS = 1

Contributing and Support