Reflected light occultations - "Point not on the unit disk."

[fbartolic]

I'm getting an error when calling map.design_matrix for a reflected light occultation. Here's a reproducible example:

xo, yo, zo, ro = 31.03953239062832, 23.892679948795926, 1.0, 39.10406741663172
xs, ys, zs = -0.1023801975400136, 0.02539410349374292, 5.423840491698775

map = starry.Map(ydeg=6, reflected = True)
map.design_matrix(xo=xo, yo=yo, zo=zo, ro=ro, xs=xs, ys=ys, zs=zs)

Traceback:

Pre-computing some matrices... Done.
Compiling `X`... /Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/theano/gof/cc.py:968: UserWarning: Your g++ compiler fails to compile OpenMP code. We know this happen with some version of the EPD mingw compiler and LLVM compiler on Mac OS X. We disable openmp everywhere in Theano. To remove this warning set the theano flags `openmp` to False.
  ret += x.c_compile_args()
Done.
/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
~/anaconda3/envs/io/lib/python3.8/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    902             outputs =\
--> 903                 self.fn() if output_subset is None else\
    904                 self.fn(output_subset=output_subset)

~/anaconda3/envs/io/lib/python3.8/site-packages/theano/gof/op.py in rval(p, i, o, n)
    891             def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 892                 r = p(n, [x[0] for x in i], o)
    893                 for o in node.outputs:

~/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/ops/integration.py in perform(self, node, inputs, outputs)
    125         dummy = np.zeros((len(b), self.N))
--> 126         outputs[0][0], _, _, _, _ = self.func(b, theta, bo, ro, dummy)
    127 

RuntimeError: Point not on the unit disk.

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
<ipython-input-3-65f7bca0e77b> in <module>
      1 map = starry.Map(ydeg=6, reflected = True)
----> 2 map.design_matrix(xo=xo, yo=yo, zo=zo, ro=ro, xs=xs, ys=ys, zs=zs)

~/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/maps.py in design_matrix(self, **kwargs)
   1780 
   1781         # Compute & return
-> 1782         return self.ops.X(
   1783             theta,
   1784             xs,

~/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/utils.py in wrapper(instance, *args)
    130 
    131             # Return the compiled version
--> 132             return getattr(instance, cname)(*args)
    133 
    134     return wrapper

~/anaconda3/envs/io/lib/python3.8/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    912                 if hasattr(self.fn, 'thunks'):
    913                     thunk = self.fn.thunks[self.fn.position_of_error]
--> 914                 gof.link.raise_with_op(
    915                     node=self.fn.nodes[self.fn.position_of_error],
    916                     thunk=thunk,

~/anaconda3/envs/io/lib/python3.8/site-packages/theano/gof/link.py in raise_with_op(node, thunk, exc_info, storage_map)
    323         # extra long error message in that case.
    324         pass
--> 325     reraise(exc_type, exc_value, exc_trace)
    326 
    327 

~/anaconda3/envs/io/lib/python3.8/site-packages/six.py in reraise(tp, value, tb)
    700                 value = tp()
    701             if value.__traceback__ is not tb:
--> 702                 raise value.with_traceback(tb)
    703             raise value
    704         finally:

~/anaconda3/envs/io/lib/python3.8/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    901         try:
    902             outputs =\
--> 903                 self.fn() if output_subset is None else\
    904                 self.fn(output_subset=output_subset)
    905         except Exception:

~/anaconda3/envs/io/lib/python3.8/site-packages/theano/gof/op.py in rval(p, i, o, n)
    890             # default arguments are stored in the closure of `rval`
    891             def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 892                 r = p(n, [x[0] for x in i], o)
    893                 for o in node.outputs:
    894                     compute_map[o][0] = True

~/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/ops/integration.py in perform(self, node, inputs, outputs)
    124         b, theta, bo, ro = inputs
    125         dummy = np.zeros((len(b), self.N))
--> 126         outputs[0][0], _, _, _, _ = self.func(b, theta, bo, ro, dummy)
    127 
    128     def grad(self, inputs, gradients):

RuntimeError: Point not on the unit disk.
Apply node that caused the error: <starry._core.ops.integration.sTReflectedOp object at 0x7fdf78bd2ac0>(AdvancedSubtensor1.0, AdvancedSubtensor1.0, AdvancedSubtensor1.0, <TensorType(float64, scalar)>)
Toposort index: 33
Inputs types: [TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, scalar)]
Inputs shapes: [(1,), (1,), (1,), ()]
Inputs strides: [(8,), (8,), (8,), ()]
Inputs values: [array([-0.99981094]), array([2.24244038]), array([39.17030414]), array(39.10406742)]
Outputs clients: [[SparseDot(<starry._core.ops.integration.sTReflectedOp object at 0x7fdf78bd2ac0>.0, SparseConstant{csc,float64,shape=(49, 49),nnz=186})]]

Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
  File "<ipython-input-3-65f7bca0e77b>", line 2, in <module>
    map.design_matrix(xo=xo, yo=yo, zo=zo, ro=ro, xs=xs, ys=ys, zs=zs)
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/maps.py", line 1782, in design_matrix
    return self.ops.X(
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/utils.py", line 125, in wrapper
    func(instance, *dummy_args),
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/core.py", line 875, in X
    sT = self.sT(b_term[i_occ], theta_term[i_occ], bo[i_occ], ro)
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/utils.py", line 105, in wrapper
    return func(instance, *args)
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/core.py", line 762, in sT
    return self._sT(b, theta, bo, ro)
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/theano/gof/op.py", line 615, in __call__
    node = self.make_node(*inputs, **kwargs)
  File "/Users/fbartolic/anaconda3/envs/io/lib/python3.8/site-packages/starry-1.0.1.dev14+gf585d6a0-py3.8-macosx-10.9-x86_64.egg/starry/_core/ops/integration.py", line 111, in make_node
    outputs = [tt.TensorType(inputs[-1].dtype, (False, False))()]

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

Starry version is '1.0.1.dev14+gf585d6a0'

[rodluger]

I'll look into it. Thanks!

[rodluger]

@fbartolic Quick update: I fixed the issue but found a new instability in this limit. Working on resolving that, should have this working by tomorrow morning.

[rodluger]

OK I fixed the instability, I think.

xo = np.linspace(29.03953239062832, 33.03953239062832, 1000)
yo, zo, ro = 23.892679948795926, 1.0, 39.10406741663172
xs, ys, zs = -0.1023801975400136, 0.02539410349374292, 5.423840491698775
map = starry.Map(ydeg=6, reflected=True)
map[1:, :] = 1
plt.plot(xo, map.flux(xo=xo, yo=yo, zo=zo, ro=ro, xs=xs, ys=ys, zs=zs))
plt.show()

@fbartolic Let me know if this works. There is still an instability when the sub-stellar point gets very close to the center of the disk -- I'm going to tackle that tomorrow.

[rodluger]

@fbartolic

OK, it looks like I've fixed all the instabilities I'm aware of. Below is a (really unrealistic) mock example of the Earth in occultation, where the radius ratio is the same as that of Jupiter/Io. I made the orbital, rotational, and illumination periods all different to break most of the degeneracies in the problem and verify that we recover the truth when doing inference (and it works!) In the bottom plot, left is truth and right is a posterior draw.

So, do play around with this when you get a chance and let me know if you run into trouble. Thanks!

import starry
import numpy as np
import matplotlib.pyplot as plt

starry.config.lazy = False
np.random.seed(3)

# Orbital/geometric parameters
npts = 50000
t = np.linspace(0, 1, npts)
porb = 0.19
prot = 0.12
rorb = 50
ro = 38.0
yo = np.sin(2 * np.pi / porb * t + 0.5)
xo = np.cos(2 * np.pi / porb * t)
zo = np.sin(2 * np.pi / porb * t)
amp = rorb / np.sqrt(xo ** 2 + yo ** 2 + zo ** 2)
xo *= amp
yo *= amp
zo *= amp
theta = 360.0 / prot * t
xs = np.sin(7 * np.pi * t)
ys = np.cos(5 * np.pi * t)
zs = 5
kwargs = dict(xs=xs, ys=ys, zs=zs, theta=theta, xo=xo, yo=yo, zo=zo, ro=ro)

# Generate a synthetic dataset
map = starry.Map(ydeg=10, reflected=True)
map.load("earth")
img0 = map.render(projection="rect", illuminate=False)
flux0 = map.flux(**kwargs)
err = 3e-4
flux = flux0 + np.random.randn(npts) * err

# Solve the linear problem & draw a sample
map.set_data(flux, C=err ** 2)
map.set_prior(L=1e-4)
map.solve(**kwargs)
map.draw()
model = map.flux(**kwargs)
img = map.render(projection="rect", illuminate=False)

# Plot the light curve model
plt.plot(t, flux)
plt.plot(t, model)

# Plot the posterior draw for the map
fig, ax = plt.subplots(1, 2)
vmin = img0.min()
vmax = img0.max()
ax[0].imshow(
    img0,
    origin="lower",
    cmap="plasma",
    extent=(-180, 180, -90, 90),
    vmin=vmin,
    vmax=vmax,
)
ax[1].imshow(
    img,
    origin="lower",
    cmap="plasma",
    extent=(-180, 180, -90, 90),
    vmin=vmin,
    vmax=vmax,
)
plt.show()