fbartolic / volcano Goto Github PK

I've been trying out a few different ways of representing uncertainty on inferred maps. Here are some examples with a linear fit to mock occultations in emitted light.

One option is to have a plot of the mean map (summed over posterior samples in pixel space) and a separate map for the standard deviation or inverse variance. This is shown below

Another option is to make uncertain parts of the map more transparent, this requires matplotlib 3.2.0 where alpha can be an array in imshow. There a few different ways of mapping uncertainties to alpha values and I've plotted a few I could think of below:

I quite like the second approach, in particular the ones where alpha is scaled by the precision (inverse variance). It makes it obvious that only one side of the occulted object was visible during occultation and all spurious features have high transparency.

All of these plots look a lot prettier in anotebook!

Let me know what you think @dfm @rodluger

Timescales of occultations are on the order of a few minutes so the finite exposure time should be taken into account in the model.

Jupiter has a fuzzy atmosphere with a scale heigh of about 27km which limits the resolution of inferred surface maps which we should take into account.

CI is broken because the figures take too long to run on Azure (>5 hrs or so). I need to think how best to optimize this.

How would we model this in Starry, second light source in addition to the Sun? Is it time variable?

AASTeX supports this, I think.

Mine is: 0000-0002-9328-5652

@rodluger I used the script from planetplanet you mentioned to generate a mock occultation light curve for the super Earth 55 Cancri e with JWST's MIRI instrument, a demo notebook is here. Could you check if I'm doing everything correctly? Here's what it looks assuming an observation cadence of 1min:

The eclipse duration is about ~2min so we only get one data point during an occultation assuming a single eclipse. I guess we could get a lot more with multiple observations? What would be a realistic observing strategy here? Are the potential targets with much longer eclipse times?

• Abstract: "fully Bayesian" remove "fully" or replace this with probabilistic.
• Throughout: see if you can move the figures so that they don't always appear right below a section heading. It becomes hard to follow!
• Bottom of page 3: "Because we have lots of prior information on..." -> "Because of high-resolution resolved imaging of Io, we..."
• Eq 10: Just put the actual Matern kernel here rather than the approximation
• Figure 7 and 8: I think the title of the second map should be "inferred map" or something right?
• Throughout: I don't love the word "superior". I'd prefer something like "better" or (even better!) specifically describe how it's "better". For example, "more rigorous", "more physically motivated", or "less restrictive"...
• Page 31: "the fact that we cannot account for" -> "the fact that we do not yet account for"
• Page 38: "lots of divergences" -> "many divergences"

A group which models mutual occultations between Galilean satellites for the purpose of improving astrometric positions of the satellites uses a reflectance law which is more sophisticated than Lambert's cosine law. Here's a quote from a recent paper:

We successfully solved these issues by adopting a generalisation of
Lambert scattering law given by Oren and Nayar (1994). The Oren-Nayer
model takes into account the direction of radiance and the roughness of
the surface in a natural way, so that the reflectance depends only of the
albedo and in one more parameter that very smoothly tunes a wide range
of surface roughness, and most importantly, regardless of the wavelength.
This model realistically reproduces the illumination of an object
in modern computer graphic scenes for movies and for the full Moon.
(Oren and Nayar, 1994). In Dias-Oliveira et al. (2013) a simplified
version of the model was used. Here, we implemented the complete
version in Oren and Nayar (1994), taking into account the direct illumination
and all inter-reflection components of the radiance.

They don't really describe how they implemented this in that paper but Dias-Oliveira et al. (2013) goes into more details. The relevant equations are 2-9 and A6, A7. The reflected intensity depends on the vector of incident and reflected light light which they get from JPL Horizons. Not sure if these integrals can be computed analytically.

We need to generate mock occultation datasets to test how well can we recover surface features. The simple thing to do is to use Starry to generate the mock datasets and add some noise, this is good for the beginning but Rodrigo is concerned about the circular reasoning issue. An alternative option is to use something healpy and compute the flux numerically over a spotted star.

Either way, we should first generate datasets using a time independent map, the albedo map provided by [USGS] (https://pubs.usgs.gov/sim/3168/) should be good enough for this purpose. We should also generate datasets using a time dependent map to model the volcanically active regions. A simple approach is to place a few bright spots on the surface at locations in (lat, lon) corresponding to major volcanos on Io's surface and model the time dependence as a SHO.

Given x observations of occultations/phase curves in reflected light and a map of fixed degree, what's the size of the nullspace? How does this scale with the degree of the map?

Relevant tutorial.

I have everything I need to compute this (JPL Horizons data), a minor issue is that Starry expects a fixed inclination and obliquity of a given map instance and this is these are time dependent for Io.

Conditioned on fixed orbital parameters and assuming Guassian priors on all SH coefficients the posterior is analytic. If we also want to sample in the nonlinear parameters we can analytically marginalize over the linear parameters to speed up inference as in this but then we get only the conditional posterior over the linear parameters (conditional on the nonlinear parameters).

The issue with both these approaches is that with Gaussian priors on the SH coefficients we can get non physical maps with negative intensity at parts of the map. We thus ignore the specific knowledge we have that a given map must have positive intensity everywhere which constrains the inference to some extent.

One solution might be to do fast analytic inference assuming dumb priors and then reweight samples under sensible priors. Importance sampling might be problematic but I found this recent paper describing an approach for doing exactly that which is superior to naive importance sampling. In the paper they have a linear regression example with similar dimensionality to a typical Starry model.

Even if this works it's still not clear to me how to incorporate this approach with the analytic marginalization over the SH coefficients since the requirement for the algorithm is to have posterior samples over the SH coefficients.

• Change plots to include full model prediction (including GP) and binned data
• Finish conclusions section
• Incorporate feedback from @dfm and @rodluger
• Incorporate feedback from Robert
• Move some of the math into appendices
• Fix build on azure.

@dfm @rodluger anything else?

Currently the paper has two somewhat distinct parts. The first part concerns the question of how to infer static maps from occultation light curves of Io with starry. It introduces pixels/sphixels as a way of ensuring positivity of the maps, it discusses different prior choices, fixing ringing artefacts via Gaussian filtering and the differences between optimization, VI and HMC for fitting the model. The second part generalizes the static model to the case where we the map for each light curve is different. We tackle this problem with NMF in a Bayesian setting. The common thread between both parts is that we use the same Io dataset as a testbed fo our models. The NMF model builds on the first part. Both parts are relevant to stars/exoplanets but the first part more so.

If we decide to split the paper into two, how should we do it? Here's a rough table of contents for both papers:

Paper 1 - Static model:

• background on Io, IRTF dataset
• computing the ephemeris, effective radius of Jupiter, information content of different kinds of light curves (occultations, phase curves, both)
• sphixels/pixels
• influence of different priors on sphixels on the inferred map
• smoothing kernel
• fitting simulated pair of light curves
• fitting multiple pairs of IRTF light curves, recovering Loki
• comparison between different inference algorithms - MAP, VI, HMC
• fitting a simulated JWST dataset of volcanic super Earth?

Possible additions:

• mutual occultations in reflected light
• comparison to parametric models where we fit for location, size and amplitude of a fixed nr. of spots
• simulating uncertainty in orbital parameters and investigating the degeneracies between those parameters and map features

Paper 2 - Dynamic model:

• repeated discussion of Io, dataset
• literature review on NMF, discussion of degeneracies
• priors for NMF
• inferring map and components for simulated data, recovering time variability in some sense, sensitivity to nr. of components, degree of map etc.
• comparison to a simpler model such as Taylor expanding SH coefficients or pixels
• apply model to IRTF data, test sensitivity to fitting different subsets of the data, plot time series of intensities for interesting specific regions on the surface

Possible additions:

• try non Bayesian approach with off the shelf optimizers for NMF such as proxmin

Given that we'd be mostly focused on Io in both papers, I'm leaning towards a paper 1/paper 2 split rather than having two distinct papers. What are your thoughts @rodluger @dfm ?

Simplest thing is to Taylor expand the SH coefficients about some reference time t0 as in this tutorial. The issue with this approach is that the intensity will blow up far from t0 and we have a very long time baseline in our data.

Another approach is to expand the SH coefficients as a Fourier series in time which avoids the issue with the intensity blowing up and it's likely a better model for the time evolution of the volcanos. If we treat the frequencies as free parameters the model is nonlinear.

Third approach is to treat each collection of occultation lightcurves closely separated in time as i.i.d. draws from some latent distribution describing the population of spots (volcanos). Instead of fitting for the SH coefficients here we'd generate a map of a specific degree and place several spots on the map as in this tutorial and infer the hyperparameters describing the spot population. The simplest version of this model would be to place a few spots at fixed (lat, lon) corresponding to locations of known Volcanos such as Loki and fit for the mean and the variance for the brightness of each spot. This could be interesting to planetary science people because we can infer properties of individual volcanos over time, account for uncertainties in a better way and have a longer time baseline. One issue I see here is the situation where we get eruptions of unknown volcanos comparable in intensity to the ones we chose to parametrize. Maybe the solution is to fit a baseline map in addition to the parametrized spots such that the baseline map accounts for all the unaccounted for volcanos?

Finally, Dan's idea is to use a GP to model the temporal and spatial (l-space) evolution of the SH coefficients with a 2D GP. If we also include wavelength it becomes a 3D GP. This is the most elegant approach but I don't know if the inference is tractable.

At the moment I have no idea if the orientation of the Io Starry map at any given time based on JPL Horizons data is correct.

@dfm @rodluger I've tried estimating the errorbars with using the following model: sigma_i_new**2 = (alpha * sigma_i)**2 + (beta * sqrt(F_i))**2 where sigma_i is the errorbar for the i-th data point (which I've previously estimated using a filtering procedure with the Savitzky-Golay filter) and F_i is the modeled flux. I place a bounded normal prior N(1,0.1) for a (a>1), and a bounded normal prior on b N(0,1) (b>0). The numpyro model is implemented here and the functions for plotting the GP predictions are here. The fit of the 1998 pair of occultations is shown below:

The residuals are scaled by their (inferred) errorbar.

Both alpha and beta seem to be well constrained from the data:

The blue histogram corresponds to the ingress light curve and the orange one to the egress light curve.

The inferred GP hyperaparamters are the following:

It seems like the errorbars are overestimated which is why the GP predictions are flat. Any ideas on how to fix this problem? It doesn't seem like the flux dependent independent noise term is working well.

It might be a bug in the code...