Tutorial Testing, New EB tutorial "fast" mode doesn't seem to run

[christinahedges]

Hey Rodrigo,

Currently testing your tutorials. FYI version numbers:

theano 1.0.4
starry 1.0.0.dev39+g56020701
pymc3 3.8

The new EB tutorials are awesome and I can't wait to use them in my research!

I'm trying to use this tutorial for the "fast" method of solving for all the parameters in a system including the map. (I'm applying it to a system where I want to find the possible phase curve of the planet). When I try to build a starry model with reasonable parameters for the system I care about I end up with an error when I try to evaluate the flux, and I'm struggling to parse the error.

Input

p1, t01 = 3.0931514, 2.8554633 + 2457000
p2, t02 =  6.7666719, 0.061385 + 2457000

with pm.Model() as model:
    PositiveNormal = pm.Bound(pm.Normal, lower=0.0)

    # Primary
    A_inc = pm.Normal("A_inc", mu=90, sd=5, testval=90)
    A_amp = 1.0
    A_r = PositiveNormal("A_r", mu=0.778, sd=0.005, testval=0.778)
    A_m = PositiveNormal("A_m", mu=0.81, sd=0.03, testval=0.81)
    A_prot = PositiveNormal("A_prot", mu=10, sd=10, testval=10)
    pri = starry.Primary(
        starry.Map(ydeg=5, udeg=2, inc=A_inc, amp=A_amp),
        r=A_r,
        m=A_m,
        prot=A_prot,
    )
    pri.map[1] = 0.6278
    pri.map[2] = 0.1045

    # Secondary
    B_inc = pm.Normal("B_inc", mu=85.05, sd=0.09, testval=85.05)
    B_amp = PositiveNormal('B_amp', mu=0.0001, sd=0.001, testval=0.0001)
    B_r = PositiveNormal("B_r", mu=0.01883983, sd=0.0005, testval=0.01883983)
    B_m = PositiveNormal("B_m", mu=1, sd=0.1, testval=1)

    B_porb = PositiveNormal("B_porb", mu=p1, sd=0.01, testval=p1)
    B_prot = B_porb
    B_t0 = pm.Normal("B_t0", mu=t01, sd=0.0001, testval=t01)

    sec = starry.Secondary(
        starry.Map(ydeg=2, udeg=0, inc=B_inc, amp=B_amp),
        r=B_r,
        m=B_m,
        porb=B_porb,
        prot=B_prot,
        t0=B_t0,
        inc=B_inc,
    )

    # System
    sys = starry.System(pri, sec)
    sys.flux(0).eval()

Error

MissingInputError: Input 0 of the graph (indices start from 0), used to compute Elemwise{exp,no_inplace}(B_r_lowerbound__), was not provided and not given a value. Use the Theano flag exception_verbosity='high', for more information on this error.

I also get the same error when I try to run the full example from the docs:

Input

map = starry.Map(ydeg=5)
map.add_spot(amp=-0.075, sigma=0.1, lat=0, lon=-30)
A_y = np.array(map.y[1:])
map.reset()
map.add_spot(amp=-0.075, sigma=0.1, lat=-30, lon=60)
B_y = np.array(map.y[1:])

A = dict(
    ydeg=5,  # degree of the map
    udeg=2,  # degree of the limb darkening
    inc=80.0,  # inclination in degrees
    amp=1.0,  # amplitude (a value prop. to luminosity)
    r=1.0,  #  radius in R_sun
    m=1.0,  # mass in M_sun
    prot=1.25,  # rotational period in days
    u=[0.40, 0.25],  # limb darkening coefficients
    y=A_y,  # the spherical harmonic coefficients
)

B = dict(
    ydeg=5,  # degree of the map
    udeg=2,  # degree of the limb darkening
    inc=80.0,  # inclination in degrees
    amp=0.1,  # amplitude (a value prop. to luminosity)
    r=0.7,  #  radius in R_sun
    m=0.7,  #  mass in M_sun
    porb=1.00,  # orbital period in days
    prot=0.625,  # rotational period in days
    t0=0.15,  # reference time in days (when it transits A)
    u=[0.20, 0.05],  # limb darkening coefficients
    y=B_y,  # the spherical harmonic coefficients
)

t = np.linspace(-2.5, 2.5, 1000)
flux_true = sys.flux(t)
sigma = 0.0005
flux = flux_true + sigma * np.random.randn(len(t))

with pm.Model() as model:
    PositiveNormal = pm.Bound(pm.Normal, lower=0.0)

    # Primary
    A_inc = pm.Normal("A_inc", mu=80, sd=5, testval=80)
    A_amp = 1.0
    A_r = PositiveNormal("A_r", mu=0.95, sd=0.1, testval=0.95)
    A_m = PositiveNormal("A_m", mu=1.05, sd=0.1, testval=1.05)
    A_prot = PositiveNormal("A_prot", mu=1.25, sd=0.01, testval=1.25)
    pri = starry.Primary(
        starry.Map(ydeg=A["ydeg"], udeg=A["udeg"], inc=A_inc, amp=A_amp),
        r=A_r,
        m=A_m,
        prot=A_prot,
    )
    pri.map[1] = A["u"][0]
    pri.map[2] = A["u"][1]

    # Secondary
    B_inc = pm.Normal("B_inc", mu=80, sd=5, testval=80)
    B_amp = 0.1
    B_r = PositiveNormal("B_r", mu=0.75, sd=0.1, testval=0.75)
    B_m = PositiveNormal("B_m", mu=0.70, sd=0.1, testval=0.70)
    B_prot = PositiveNormal("B_prot", mu=0.625, sd=0.01, testval=0.625)
    B_porb = PositiveNormal("B_porb", mu=1.01, sd=0.01, testval=1.01)
    B_t0 = pm.Normal("B_t0", mu=0.15, sd=0.001, testval=0.15)
    sec = starry.Secondary(
        starry.Map(ydeg=B["ydeg"], udeg=B["udeg"], inc=B_inc, amp=B_amp),
        r=B_r,
        m=B_m,
        porb=B_porb,
        prot=B_prot,
        t0=B_t0,
        inc=B_inc,
    )
    sec.map[1] = B["u"][0]
    sec.map[2] = B["u"][1]

    # System
    sys = starry.System(pri, sec)
    sys.flux(t).eval()

Error

MissingInputError: Input 0 of the graph (indices start from 0), used to compute Elemwise{exp,no_inplace}(B_r_lowerbound__), was not provided and not given a value. Use the Theano flag exception_verbosity='high', for more information on this error.

[rodluger]

Hi @christinahedges, cool! The issue is probably this line

sys.flux(t).eval()

Within a pymc3 model context, none of the Theano variables ever have values ("inputs"), so that's the error you get when you try to evaluate them. pymc3 stores things differently (using test values), so I think what you want to do is

xo.eval_in_model(sys.flux(t))

which is a convenience function DFM wrote to allow you to evaluate things inside of models. If you don't pass the point kwarg to the function, it will evaluate things at the test point, which is what I think you want. Later, once you have a trace from sampling, you can evaluate the ith sample by running

xo.eval_in_model(sys.flux(t), point=trace.point(i))

Let me know if this fixes the issue.

[christinahedges]

Awesome, this definitely does the trick. I think I had not read the earlier tutorials in enough detail to catch this change! I hope you don't mind if I post things I find as I go through and use the tutorials?

A few other small things so far:

1. In the tutorial, we're using known data where we have the "True" flux with the correct normalization. When I get Kepler/TESS data, I tend to have to normalize it to the median value. Could you show in the tutorial for the eclipsing binary how to add a "mean" offset term? (I think the exoplanet tutorials do this)

2. In this tutorial you refer to a MAP variable, but I can't find where it's defined in the tutorial.

map = starry.Map(ydeg=A["ydeg"])
map.inc = map_soln["A_inc"]
map[1:, :] = MAP[0]
map.show(theta=np.linspace(0, 360, 50))

I'm trying to run it so I can see my converged MAPs

[rodluger]

Yes, keep the comments coming!

1. You probably want to solve for the amp property of the two maps, which I'm currently fixing at 1.0 and 0.1 for the primary and secondary, respectively. This used to be called L (for luminosity), but it's unitless, so I changed it to an "amplitude". I'll add a note about this in the tutorial.

2. You should replace MAP with yhat; I must have changed the name of that variable after running that cell -- sorry!

I found some other small issues with the tutorial you're playing with, so I'm re-running the notebook and will update it soon.

[christinahedges]

1. I think you're right, I want to solve for amp, but if my amp_primary=1 and amp_secondary=0.1 then won't my mean flux be 1.1? if I put in a light curve which is normalized to 1, how is starry going to be able to fit it if it is not allowed an offset/mean term?

2. Awesome thanks!

[rodluger]

I think that's fine: you just need to make pri.map.amp and sec.map.amp into pymc3 variables and give them suitable priors. Then you'll get some posterior over them and hopefully the means will be something like 1/1.1 ~ 0.91 and 0.1/1.1 ~ 0.09. Right?

[christinahedges]

Oh I see! I think I was confused because in previous iterations the star luminosity was fixed at "1", now that the star can change I see that relative brightnesses are fine. Ok excellent, I'll keep digging!

[rodluger]

Ah, yes, there's no longer that restriction. Cool!

[rodluger]

@christinahedges Just a heads up, since this affects the notebook you're working on. I just pushed a version to master in which I changed how you sample from the posterior over surface maps in Keplerian systems. Previously, in order to set the map coefficients of the primary and secondary to a random posterior sample, I was running

sys.solve(t)
pri.map.draw()
sec.map.draw()

Now, I added a draw method to the system object itself, so the proper thing to do now is just

sys.solve(t)
sys.draw()

and that will automatically populate all the bodies' spherical harmonic coefficient vectors.

[christinahedges]

I have another quick question that's come up while I've been testing, it may be that this is already documented but I haven't seen it:

When I make a model with two secondaries, starry throws an error if I try to build the model where the secondaries have different ydeg parameters (e.g. if I want to fit for a surface map of just one of the secondaries, this breaks). Is this expected behavior?

Thanks!!! C

[rodluger]

Interesting. No, that's not expected. Can you post a MWE?

On Mon, Dec 16, 2019 at 6:17 PM Christina Hedges notifications@github.com wrote:

I have another quick question that's come up while I've been testing, it may be that this is already documented but I haven't seen it:

When I make a model with two secondaries, starry throws an error if I try to build the model where the secondaries have different ydeg parameters (e.g. if I want to fit for a surface map of just one of the secondaries, this breaks). Is this expected behavior?

Thanks!!! C

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/223?email_source=notifications&email_token=ACHEKKYG64CLZSVLDYL7RATQZAEAHA5CNFSM4JZYLZO2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHAPAGI#issuecomment-566292505, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACHEKKZCRQFWCEFZHFT6KDTQZAEAHANCNFSM4JZYLZOQ .

[christinahedges]

This is a fairly minimal working example

import numpy as np
import starry
import pymc3 as pm
import exoplanet as xo
import corner

starry.config.lazy = True
starry.config.quiet = True
t = np.linspace(0, 10, 10000)

p1, t01 = 3.0931514, 2.8554633 + 2457000
p2, t02 =  6.7666719, 0.061385 + 2457000

with pm.Model() as model:
    PositiveNormal = pm.Bound(pm.Normal, lower=0.0)

    # Primary
#    A_inc = pm.Normal("A_inc", mu=90, sd=5, testval=90)
    A_amp = PositiveNormal('A_amp', mu=1, sd=0.001, testval=1)
    A_r = PositiveNormal("A_r", mu=0.778, sd=0.05, testval=0.778)
    A_m = PositiveNormal("A_m", mu=0.81, sd=0.03, testval=0.81)
    A_prot = PositiveNormal("A_prot", mu=10, sd=10, testval=10)
    pri = starry.Primary(
        starry.Map(ydeg=4, udeg=2, inc=90, amp=A_amp),
        r=A_r,
        m=A_m,
        prot=A_prot,
    )
    pri.map[1] = 0.6278
    pri.map[2] = 0.1045

    # Secondary
    B_inc = pm.Normal("B_inc", mu=90, sd=2, testval=90)
    B_amp = PositiveNormal('B_amp', mu=0.0001, sd=0.001, testval=0.00001)
    B_r = PositiveNormal("B_r", mu=0.01883983, sd=0.005, testval=0.01883983)
    #B_m = PositiveNormal("B_m", mu=1, sd=0.1, testval=1)

    B_porb = PositiveNormal("B_porb", mu=p1, sd=0.0001, testval=p1)
    B_prot = B_porb
    B_t0 = pm.Normal("B_t0", mu=t01, sd=0.000001, testval=t01)

    sec = starry.Secondary(
        starry.Map(ydeg=2, udeg=0, inc=B_inc, amp=B_amp),
        r=B_r,
#        m=B_m,
        porb=B_porb,
        prot=B_prot,
        t0=B_t0,
        inc=B_inc,
    )

    #     # Secondary
    C_inc = pm.Normal("C_inc", mu=90, sd=2, testval=90)
#    C_amp = PositiveNormal('C_amp', mu=0.0001, sd=0.001, testval=0.00001)
    C_r = PositiveNormal("C_r", mu=0.01778585, sd=0.005, testval=0.01777)
#    C_m = PositiveNormal("C_m", mu=1, sd=0.1, testval=1)

    C_porb = PositiveNormal("C_porb", mu=p2, sd=0.0001, testval=p2)
    C_t0 = pm.Normal("C_t0", mu=t02, sd=0.000001, testval=t02)

    other = starry.Secondary(
        starry.Map(ydeg=0, udeg=0, inc=90, amp=0),
        r=C_r,
 #       m=C_m,
        porb=C_porb,
        t0=C_t0,
        inc=C_inc,
    )

    # System
    sys = starry.System(pri, sec, other) 
    mf = xo.eval_in_model(sys.flux(t))

[rodluger]

Awesome, thanks for catching this bug. I believe it's fixed now. The issue was that I was trying to tile the design matrices of each body together, but you can only do that if they are the same shape!

Can you pull again from master and confirm that it's working? Thanks!

[rodluger]

@christinahedges BTW, I'm going to change the linear solve a bit in the next few days (sorry!) It turns out it's easy to linearly solve for/marginalize over the map amplitude. The main difference for the user is that the prior for the map coefficients will now have Ny (instead of Ny - 1) terms, where the first term is the prior on the amplitude. The remaining terms are the prior on the spherical harmonic coefficients weighted by the amplitude. I've fixed all the tutorials to explain this in more detail on the dev branch and will let you know once it's merged.

[christinahedges]

Awesome, I'm away on vacation now but I'll check this out when I have a chance. Looking forward to the update. While I was testing I was having problems with the MCMC getting "stuck" after ~40 iterations. I will wait for the new updates, make a MWE and let you know if I still have the problem!

[rodluger]

@christinahedges Are you still having issues with the sampler getting stuck? I've run into this before, but haven't been able to track down the source. It doesn't seem to be an issue with starry; rather, it has something to do with how pymc3 handles parallelization. If I specify cores=1 in the call to pm.sample, the chain runs smoothly. The issue also went away when I re-installed Anaconda. Can you try running on a single core to see if that helps?

[rodluger]

@christinahedges Any updates on this issue?

[christinahedges]

Hi the original issue is solved, and setting cores to 1 does seem to stop the MCMC getting stuck. 👍