Symbolic Equations with Sympy#
Every major class in spotgp has a get_sympy() method that renders its defining equations as LaTeX. This is useful for:
Inspecting the analytic form of your model
Verifying parameter choices
Including equations in papers or presentations
This tutorial demonstrates get_sympy() for each component class.
import sys
sys.path.append("../..")
import numpy as np
from spotgp import (
TrapezoidSymmetricEnvelope,
TrapezoidAsymmetricEnvelope,
SkewedGaussianEnvelope,
ExponentialEnvelope,
VisibilityFunction,
LatitudeDistributionFunction,
SpotEvolutionModel,
DeltaDistribution,
UniformDistribution,
GaussianDistribution,
LogNormalDistribution,
)
1. Envelope functions#
Each envelope displays \(\Gamma(t)\), \(\hat{\Gamma}(\omega)\), and \(R_\Gamma(\tau)\).
Symmetric trapezoid#
env = TrapezoidSymmetricEnvelope(lspot=10.0, tau_spot=3.0)
env.get_sympy();
\[\displaystyle \textbf{TrapezoidSymmetricEnvelope}\]
\[\begin{split}\displaystyle \Gamma(t) = \begin{cases} 0 & \text{for}\: t < - \frac{\ell}{2} - \tau_{\rm spot} \\\frac{\left(\frac{\ell}{2} + \tau_{\rm spot} + t\right)^{2}}{\tau_{\rm spot}^{2}} & \text{for}\: t < - \frac{\ell}{2} \\1 & \text{for}\: t \leq \frac{\ell}{2} \\\frac{\left(\frac{\ell}{2} + \tau_{\rm spot} - t\right)^{2}}{\tau_{\rm spot}^{2}} & \text{for}\: t < \frac{\ell}{2} + \tau_{\rm spot} \\0 & \text{otherwise} \end{cases}\end{split}\]
\[\displaystyle \hat{\Gamma}(\omega) = \frac{4 \left(\omega \tau_{\rm spot} \cos{\left(\frac{\ell \omega}{2} \right)} + \sin{\left(\frac{\ell \omega}{2} \right)} - \sin{\left(\frac{\ell \omega}{2} + \omega \tau_{\rm spot} \right)}\right)}{\omega^{3} \tau_{\rm spot}^{2}}\]
\[\displaystyle R_{\Gamma}(\tau) = \int\limits_{0}^{\infty} \Gamma{\left(t \right)} \Gamma{\left(\tau + t \right)}\, dt \quad \text{[numerical]}\]
Asymmetric trapezoid#
env_asym = TrapezoidAsymmetricEnvelope(lspot=10.0, tau_em=2.0, tau_dec=5.0)
env_asym.get_sympy();
\[\displaystyle \textbf{TrapezoidAsymmetricEnvelope}\]
\[\begin{split}\displaystyle \Gamma(t) = \begin{cases} 0 & \text{for}\: t < - \frac{\ell}{2} - \tau_{\rm em} \\\frac{\left(\frac{\ell}{2} + \tau_{\rm em} + t\right)^{2}}{\tau_{\rm em}^{2}} & \text{for}\: t < - \frac{\ell}{2} \\1 & \text{for}\: t \leq \frac{\ell}{2} \\\frac{\left(\frac{\ell}{2} + \tau_{\rm dec} - t\right)^{2}}{\tau_{\rm dec}^{2}} & \text{for}\: t < \frac{\ell}{2} + \tau_{\rm dec} \\0 & \text{otherwise} \end{cases}\end{split}\]
\[\displaystyle \hat{\Gamma}(\omega) = \int\limits_{-\infty}^{\infty} \Gamma{\left(t \right)} e^{- i \omega t}\, dt \quad \text{[numerical]}\]
\[\displaystyle R_{\Gamma}(\tau) = \int\limits_{0}^{\infty} \Gamma{\left(t \right)} \Gamma{\left(\tau + t \right)}\, dt \quad \text{[numerical]}\]
Skewed Gaussian#
env_sn = SkewedGaussianEnvelope(sigma_sn=3.0, n_sn=-2.0)
env_sn.get_sympy();
\[\displaystyle \textbf{SkewedGaussianEnvelope}\]
\[\displaystyle \Gamma(t) = \text{[numerical]}\]
\[\displaystyle \hat{\Gamma}(\omega) = \int\limits_{-\infty}^{\infty} \Gamma{\left(t \right)} e^{- i \omega t}\, dt \quad \text{[numerical]}\]
\[\displaystyle R_{\Gamma}(\tau) = \int\limits_{0}^{\infty} \Gamma{\left(t \right)} \Gamma{\left(\tau + t \right)}\, dt \quad \text{[numerical]}\]
Exponential#
env_exp = ExponentialEnvelope(tau_spot=3.0)
env_exp.get_sympy();
\[\displaystyle \textbf{ExponentialEnvelope}\]
\[\displaystyle \Gamma(t) = e^{- \frac{\left|{t}\right|}{\tau_{\rm spot}}}\]
\[\displaystyle \hat{\Gamma}(\omega) = \frac{2 \tau_{\rm spot}}{\omega^{2} \tau_{\rm spot}^{2} + 1}\]
\[\displaystyle R_{\Gamma}(\tau) = \left(\tau_{\rm spot} + \left|{\tau}\right|\right) e^{- \frac{\left|{\tau}\right|}{\tau_{\rm spot}}}\]
Deriving symbolic expressions#
Pass compute_symbolic=True to attempt symbolic derivation of \(\hat{\Gamma}(\omega)\) and \(R_\Gamma(\tau)\) from \(\Gamma(t)\) via sympy integration (can be slow for complex envelopes).
See also
This can be useful for custom envelopes where you only have specified the analytic equation for \(\Gamma(t)\). See: Custom Envelope Functions
env_exp.get_sympy(compute_symbolic=True);
\[\displaystyle \textbf{ExponentialEnvelope}\]
\[\displaystyle \Gamma(t) = e^{- \frac{\left|{t}\right|}{\tau_{\rm spot}}}\]
\[\displaystyle \hat{\Gamma}(\omega) = \frac{2 \tau_{\rm spot}}{\omega^{2} \tau_{\rm spot}^{2} + 1}\]
\[\displaystyle R_{\Gamma}(\tau) = \left(\tau_{\rm spot} + \left|{\tau}\right|\right) e^{- \frac{\left|{\tau}\right|}{\tau_{\rm spot}}}\]
2. Visibility function#
Displays the rotation frequency \(\omega_0(\phi)\), intermediate geometry variables (\(a_0\), \(a_1\), \(\theta_v\)), and Fourier coefficients (\(c_0\), \(c_1\), \(c_n\)).
vis = VisibilityFunction(peq=8.0, kappa=0.3, inc=np.pi / 3)
vis.get_sympy();
\[\displaystyle \textbf{VisibilityFunction}\]
\[\displaystyle \omega_0(\phi) = \frac{2 \pi \left(- \kappa \sin^{2}{\left(\phi \right)} + 1\right)}{P_{\rm eq}}\]
\[\displaystyle a_0 = \sin{\left(\phi \right)} \cos{\left(i \right)}\]
\[\displaystyle a_1 = \sin{\left(i \right)} \cos{\left(\phi \right)}\]
\[\displaystyle \theta_v = \operatorname{acos}{\left(- \frac{a_{0}}{a_{1}} \right)}\]
\[\displaystyle c_0 = \frac{\theta_{v} a_{0} + a_{1} \sin{\left(\theta_{v} \right)}}{\pi}\]
\[\displaystyle c_1 = \frac{a_{0} \sin{\left(\theta_{v} \right)} + \frac{a_{1} \left(\theta_{v} + \sin{\left(\theta_{v} \right)} \cos{\left(\theta_{v} \right)}\right)}{2}}{\pi}\]
\[\displaystyle c_n \; (n \geq 2) = \frac{\frac{a_{0} \sin{\left(\theta_{v} n \right)}}{n} + \frac{a_{1} \left(\frac{\sin{\left(\theta_{v} \left(n + 1\right) \right)}}{n + 1} + \frac{\sin{\left(\theta_{v} \left(n - 1\right) \right)}}{n - 1}\right)}{2}}{\pi}\]
3. Latitude distribution#
Displays the latitude PDF \(p(\phi)\).
# Default: uniform
lat = LatitudeDistributionFunction()
lat.get_sympy();
\[\displaystyle \textbf{LatitudeDistributionFunction}\]
\[\displaystyle p(\phi) = 1\]
import sympy as sp
class GaussianLatitude(LatitudeDistributionFunction):
def __init__(self, sigma_deg=20.0):
self.sigma = np.radians(sigma_deg)
def __call__(self, phi):
return np.exp(-0.5 * (phi / self.sigma) ** 2)
def sympy_pdf(self):
phi = sp.Symbol(r'\phi', real=True)
sigma = sp.Float(self.sigma)
return sp.exp(-sp.Rational(1, 2) * (phi / sigma) ** 2)
lat_gauss = GaussianLatitude(sigma_deg=20.0)
lat_gauss.get_sympy(status="Gaussian, sigma=20 deg");
\[\displaystyle \textbf{GaussianLatitude} \text{[Gaussian, sigma=20 deg]}\]
\[\displaystyle p(\phi) = e^{- 4.10350793751468 \phi^{2}}\]
4. Parameter distributions#
Each ParameterDistribution subclass can display its PDF. The var_name keyword controls the variable symbol.
DeltaDistribution(0.01).get_sympy(var_name=r"\sigma_k");
\[\displaystyle \textbf{DeltaDistribution}\]
\[\displaystyle p(\sigma_k) = \delta\left(x - 0.01\right)\]
UniformDistribution(3.0, 12.0).get_sympy(var_name=r"\ell_{\rm spot}");
\[\displaystyle \textbf{UniformDistribution}\]
\[\displaystyle p(\ell_{\rm spot}) = 0.111111111111111\]
GaussianDistribution(mu=5.0, sigma=1.0).get_sympy(var_name=r"\tau_{\rm spot}");
\[\displaystyle \textbf{GaussianDistribution}\]
\[\displaystyle p(\tau_{\rm spot}) = \frac{0.5 \sqrt{2} e^{- \frac{\left(1.0 x - 5.0\right)^{2}}{2}}}{\sqrt{\pi}}\]
LogNormalDistribution(mu=-4.6, sigma=0.3).get_sympy(var_name=r"\sigma_k");
\[\displaystyle \textbf{LogNormalDistribution}\]
\[\displaystyle p(\sigma_k) = \frac{1.66666666666667 \sqrt{2} e^{- \frac{\left(3.33333333333333 \log{\left(x \right)} + 15.3333333333333\right)^{2}}{2}}}{\sqrt{\pi} x}\]
Custom distribution with sympy#
Override sympy_pdf() in your subclass to get LaTeX rendering:
from spotgp import ParameterDistribution
from scipy.special import beta as beta_func
class BetaDistribution(ParameterDistribution):
"""Beta distribution scaled to [lo, hi]."""
def __init__(self, a, b, lo=0.0, hi=1.0):
self.a, self.b = float(a), float(b)
self.lo, self.hi = float(lo), float(hi)
self._B = beta_func(self.a, self.b)
@property
def support(self):
return (self.lo, self.hi)
def __call__(self, x):
t = (x - self.lo) / (self.hi - self.lo)
if t <= 0 or t >= 1:
return 0.0
return t ** (self.a - 1) * (1 - t) ** (self.b - 1) / self._B
def sympy_pdf(self):
x = sp.Symbol("x")
a, b = sp.Float(self.a), sp.Float(self.b)
lo, hi = sp.Float(self.lo), sp.Float(self.hi)
t = (x - lo) / (hi - lo)
return t ** (a - 1) * (1 - t) ** (b - 1) / (
(hi - lo) * sp.beta(a, b))
BetaDistribution(a=2, b=5, lo=2.0, hi=15.0).get_sympy(
var_name=r"\ell_{\rm spot}");
\[\displaystyle \textbf{BetaDistribution}\]
\[\displaystyle p(\ell_{\rm spot}) = 2.30769230769231 \left(1.15384615384615 - 0.0769230769230769 x\right)^{4.0} \left(0.0769230769230769 x - 0.153846153846154\right)^{1.0}\]
5. Full model#
SpotEvolutionModel.get_sympy() renders all components in sequence: envelope, visibility, and latitude distribution.
model = SpotEvolutionModel(
envelope=TrapezoidSymmetricEnvelope(lspot=10.0, tau_spot=3.0),
visibility=VisibilityFunction(peq=8.0, kappa=0.3, inc=np.pi / 3),
sigma_k=0.01,
latitude_distribution=GaussianLatitude(sigma_deg=20.0),
)
model.get_sympy();
\[\displaystyle \textbf{TrapezoidSymmetricEnvelope} \text{[user defined]}\]
\[\begin{split}\displaystyle \Gamma(t) = \begin{cases} 0 & \text{for}\: t < - \frac{\ell}{2} - \tau_{\rm spot} \\\frac{\left(\frac{\ell}{2} + \tau_{\rm spot} + t\right)^{2}}{\tau_{\rm spot}^{2}} & \text{for}\: t < - \frac{\ell}{2} \\1 & \text{for}\: t \leq \frac{\ell}{2} \\\frac{\left(\frac{\ell}{2} + \tau_{\rm spot} - t\right)^{2}}{\tau_{\rm spot}^{2}} & \text{for}\: t < \frac{\ell}{2} + \tau_{\rm spot} \\0 & \text{otherwise} \end{cases}\end{split}\]
\[\displaystyle \hat{\Gamma}(\omega) = \frac{4 \left(\omega \tau_{\rm spot} \cos{\left(\frac{\ell \omega}{2} \right)} + \sin{\left(\frac{\ell \omega}{2} \right)} - \sin{\left(\frac{\ell \omega}{2} + \omega \tau_{\rm spot} \right)}\right)}{\omega^{3} \tau_{\rm spot}^{2}}\]
\[\displaystyle R_{\Gamma}(\tau) = \int\limits_{0}^{\infty} \Gamma{\left(t \right)} \Gamma{\left(\tau + t \right)}\, dt \quad \text{[numerical]}\]
\[\displaystyle \textbf{VisibilityFunction} \text{[user defined]}\]
\[\displaystyle \omega_0(\phi) = \frac{2 \pi \left(- \kappa \sin^{2}{\left(\phi \right)} + 1\right)}{P_{\rm eq}}\]
\[\displaystyle a_0 = \sin{\left(\phi \right)} \cos{\left(i \right)}\]
\[\displaystyle a_1 = \sin{\left(i \right)} \cos{\left(\phi \right)}\]
\[\displaystyle \theta_v = \operatorname{acos}{\left(- \frac{a_{0}}{a_{1}} \right)}\]
\[\displaystyle c_0 = \frac{\theta_{v} a_{0} + a_{1} \sin{\left(\theta_{v} \right)}}{\pi}\]
\[\displaystyle c_1 = \frac{a_{0} \sin{\left(\theta_{v} \right)} + \frac{a_{1} \left(\theta_{v} + \sin{\left(\theta_{v} \right)} \cos{\left(\theta_{v} \right)}\right)}{2}}{\pi}\]
\[\displaystyle c_n \; (n \geq 2) = \frac{\frac{a_{0} \sin{\left(\theta_{v} n \right)}}{n} + \frac{a_{1} \left(\frac{\sin{\left(\theta_{v} \left(n + 1\right) \right)}}{n + 1} + \frac{\sin{\left(\theta_{v} \left(n - 1\right) \right)}}{n - 1}\right)}{2}}{\pi}\]
\[\displaystyle \textbf{GaussianLatitude} \text{[user defined]}\]
\[\displaystyle p(\phi) = e^{- 4.10350793751468 \phi^{2}}\]
Summary#
Class |
|
|---|---|
|
\(\Gamma(t)\), \(\hat{\Gamma}(\omega)\), \(R_\Gamma(\tau)\) |
|
\(\omega_0(\phi)\), \(a_0\), \(a_1\), \(\theta_v\), \(c_0\), \(c_1\), \(c_n\) |
|
\(p(\phi)\) |
|
\(p(x)\) (with custom |
|
All of the above combined |
To add sympy support to a custom subclass, override sympy_pdf() (for distributions/latitude) or sympy_Gamma() / sympy_Gamma_hat() / sympy_R_Gamma() (for envelopes).