rodluger / starry

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

map.flux() Normalization problem #316

Closed Shailen9901 closed 4 months ago

Shailen9901 commented 4 months ago

I am currently trying to model a spot to some data that I have. However, the data is demeaned relative flux:

demeaned relative flux = (relative flux / mean relative flux)

So my data has a mean of 1 with the other values representing the fractional difference from this mean.

When using map.flux(), 1 is normalized so that it is the unspotted star, so currently the plots using map.flux() look like this:

Normal Plot

But by dividing by the map.flux() by its mean I can get the plot to be the same as my data:

Demeaned Plot

However when trying to apply this to define the spot model I get an error.

with pm.Model() as model:

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

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

    # Compute the flux model
    flux_model = map.flux(theta=360.0 / P * t)
    flux_model = flux_model/np.mean(flux_model) # error here

    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)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-3f5c185e6187> in <module>
     14     # Compute the flux model
     15     flux_model = map.flux(theta=360.0 / P * t)
---> 16     flux_model = flux_model/np.mean(flux_model)
     17     print(flux_model)
     18     pm.Deterministic("flux_model", flux_model)

<__array_function__ internals> in mean(*args, **kwargs)

~\anaconda3\lib\site-packages\numpy\core\fromnumeric.py in mean(a, axis, dtype, out, keepdims, where)
   3436             pass
   3437         else:
-> 3438             return mean(axis=axis, dtype=dtype, out=out, **kwargs)
   3439 
   3440     return _methods._mean(a, axis=axis, dtype=dtype,

TypeError: mean() got an unexpected keyword argument 'out'

Is there a way to change the flux model so that it is a plot around 1 as the mean (demeaned)?

dfm commented 4 months ago

This error is happening because it's not possible to use numpy functions on Theano tensor objects. In this case, you should be able to add the following import:

import theano.tensor as tt

then update the line with the error to:

-     flux_model = flux_model/np.mean(flux_model)
+     flux_model = flux_model/tt.mean(flux_model)
Shailen9901 commented 4 months ago

Thank you so much!!! The code is now producing spots that fit the new flux models (demeaned relative flux) to my data.