rodluger / starry

Tools for mapping stars and planets.
https://starry.readthedocs.io
MIT License
144 stars 33 forks source link

Issue with starry using pymc3 priors #318

Open ar1326 opened 4 months ago

ar1326 commented 4 months ago

Hi,

I was going through the Starry tutorial for the SpotSolver (https://starry.readthedocs.io/en/latest/notebooks/SpotSolve/) and ran the block of code to determine the posterior constraints.

import numpy as np
import starry 
import pymc3 as pm 
import pymc3_ext as pmx

truth = dict(contrast=0.25, radius=20, lat=30, lon=30)

inc = 60.0
P = 1.0

map0 = starry.Map(15)
map0.inc = inc
map0.spot(
    contrast=truth["contrast"],
    radius=truth["radius"],
    lat=truth["lat"],
    lon=truth["lon"],
)
map0.show()

t = np.linspace(0, 3.0, 500)
flux0 = map0.flux(theta=360.0 / P * t)
np.random.seed(0)
flux_err = 2e-4
flux = flux0 + flux_err * np.random.randn(len(t))

with pm.Model() as model:

    # Priors
    contrast = pm.Uniform(name="contrast", lower=0.0, upper=1.0, testval=0.5)
    radius = pm.Uniform(name="radius", lower=10.0, upper=60.0, testval=15.0)
    lat = pm.Uniform(name="lat", lower=-90.0, upper=90.0, testval=0.1)
    lon = pm.Uniform(name="lon", lower=-180.0, upper=180.0, testval=0.1)

    # Instantiate the map and add the spot
    map0 = starry.Map(ydeg=15)
    map0.inc = inc
    map0.spot(contrast=contrast, lat=lat, lon=lon, radius=radius)

    # Compute the flux model
    flux_model = map0.flux(theta=360.0 / P * t)

    pm.Deterministic("flux_model", flux_model)

    # Save our initial guess
    flux_model_guess = pmx.eval_in_model(flux_model)

    # The likelihood function assuming known Gaussian uncertainty
    pm.Normal("obs", mu=flux_model, sd=flux_err, observed=flux)

However, an exception was raised, specifically TypeError: float() argument must be a string or a number, not 'TransformedRV', and I think the following details indicate that it's an issue with the radius prior:

 File ~/opt/anaconda3/lib/python3.9/site-packages/starry/maps.py:1534 in spot
    radius = self._math.cast(radius) * self._angle_factor

  File ~/opt/anaconda3/lib/python3.9/site-packages/starry/_core/math.py:165 in cast
    return np.array(args[0], dtype=floatX)

ValueError: setting an array element with a sequence.

I was wondering if there is an issue with starry compatibility with pymc3, and how to resolve this issue, please? It's just the standard tutorial code, without any adjustments (as far as I can see), so I'm unsure as to why this exception was raised. Thanks for any help!