Quickstart

Quickstart#

This tutorial walks through the core spotgp workflow in five steps:

  1. Define a spot evolution model — choose an envelope, visibility function, and amplitude

  2. Simulate a lightcurve — generate synthetic photometry from the model

  3. Build the analytic kernel — compute the GP covariance from the physics

  4. Fit the GP to data — find the MAP hyperparameters with GPSolver

  5. Inspect the fit — compare the model prediction, ACF, and PSD to the data

import numpy as np
import matplotlib.pyplot as plt

from spotgp import (
    TrapezoidSymmetricEnvelope,
    VisibilityFunction,
    SpotEvolutionModel,
    LightcurveModel,
    AnalyticKernel,
    GPSolver,
)

np.random.seed(42)

1. Define a spot evolution model#

A SpotEvolutionModel combines three components:

Component

What it controls

Key parameters

Envelope

Spot size vs. time

lspot (plateau duration), tau_spot (rise/decay timescale)

Visibility

Stellar geometry

peq (rotation period), kappa (differential rotation), inc (inclination)

Amplitude

Kernel scaling

sigma_k

 
envelope = TrapezoidSymmetricEnvelope(lspot=14.0, tau_spot=6.0)
visibility = VisibilityFunction(peq=4.0, kappa=0.3, inc=np.pi / 3)
model = SpotEvolutionModel(envelope=envelope, visibility=visibility, sigma_k=0.01)

print(model)

2. Simulate a lightcurve#

LightcurveModel.from_spot_model places nspot spots at random longitudes, latitudes, and emergence times, then sums their flux deficits over the observation window. We add Gaussian noise to mimic a real observation.

lc = LightcurveModel.from_spot_model(
    spot_model=model,
    nspot=30,
    tsim=150,
    tsamp=1.0,
)

# Add observation noise
sigma_n = 0.3 * np.std(lc.flux)
flux_obs = lc.flux + np.random.normal(0, sigma_n, lc.flux.shape)
flux_err = np.full_like(flux_obs, sigma_n)

lc.plot_lightcurve()
plt.show()

3. Build the analytic kernel#

AnalyticKernel computes the GP covariance function directly from the spot evolution physics. Call build_jax() once to pre-compile the XLA kernels — subsequent evaluations are fast.

ak = AnalyticKernel(model).build_jax()

# Evaluate and plot the kernel
tlags = np.linspace(0, 3 * model.peq, 300)
K = ak.kernel(tlags)

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(tlags, K, color="steelblue")
ax.set_xlabel("Time lag [days]")
ax.set_ylabel("K(τ)")
ax.set_title("Analytic GP kernel")
plt.tight_layout()
plt.show()

4. Fit the GP to data#

GPSolver builds the banded covariance matrix from the kernel and finds the maximum a posteriori (MAP) hyperparameters via L-BFGS-B. Define parameter bounds, compile with build_jax(), and run fit_map().

bounds = {
    "peq":         (2.0, 10.0),
    "kappa":       (-0.5, 0.8),
    "inc":         (0.1, np.pi / 2),
    "lspot":       (2.0, 25.0),
    "tau_spot":    (0.5, 12.0),
    "log_sigma_k": (-4.0, 0.0),
}

gp = GPSolver(lc.t, flux_obs, flux_err, model, bounds=bounds).build_jax()

theta_map, opt_result = gp.fit_map(nopt=3)

# Compare MAP estimates to true values
theta_true = dict(peq=4.0, kappa=0.3, inc=np.pi / 3,
                  lspot=14.0, tau_spot=6.0, sigma_k=0.01)

print(f"{'param':>12s}  {'true':>8s}  {'MAP':>8s}")
print("-" * 34)
for k, v_true in theta_true.items():
    v_map = theta_map.get(k, theta_map.get("log_" + k, float("nan")))
    print(f"{k:>12s}  {v_true:8.4f}  {v_map:8.4f}")

5. Inspect the fit#

GPSolver includes built-in diagnostic plots: the GP posterior prediction over the data, the autocorrelation function (ACF), and the power spectral density (PSD).

fig, axes = plt.subplots(3, 1, figsize=(10, 10))

gp.plot_prediction(theta=theta_map, ax=axes[0])
gp.plot_acf(theta=theta_map, ax=axes[1], normalize=True)
gp.plot_psd(theta=theta_map, ax=axes[2])

fig.tight_layout()
plt.show()

Next steps#