Point source fitting fails with image PSF?

[hcferguson]

I tried the quick-start guide with just the point source. It works if I use the Gaussian model PSF, but if I swap in the PSFs that were read in earlier in the example (by assigning model_psf = psfs it fails. Traceback below.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<timed eval> in <module>

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/blend.py in fit(self, max_iter, e_rel, **alg_kwargs)
     97             V=V,
     98             Vhat=Vhat,
---> 99             **alg_kwargs
    100         )
    101 

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/proxmin/algorithms.py in adaprox(X, grad, step, prox, scheme, b1, b2, eps, check_convergence, p, e_rel, max_iter, prox_max_iter, M, V, Vhat, callback)
    331         try:
    332             callback(*X, it=it)
--> 333             G = _as_tuple(grad(*X))
    334             Alpha = _as_tuple(step(*X, it=it))
    335             if check_convergence:

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/blend.py in <lambda>(*X)
     64             x.prior(x.view(np.ndarray)) if x.prior is not None else 0 for x in X
     65         )
---> 66         _grad = lambda *X: tuple(l + p for l, p in zip(grad_logL(*X), grad_logP(*X)))
     67         _step = lambda *X, it: tuple(
     68             x.step(x, it=it) if hasattr(x.step, "__call__") else x.step for x in X

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/autograd/wrap_util.py in nary_f(*args, **kwargs)
     18             else:
     19                 x = tuple(args[i] for i in argnum)
---> 20             return unary_operator(unary_f, x, *nary_op_args, **nary_op_kwargs)
     21         return nary_f
     22     return nary_operator

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/autograd/differential_operators.py in grad(fun, x)
     23     arguments as `fun`, but returns the gradient instead. The function `fun`
     24     should be scalar-valued. The gradient has the same type as the argument."""
---> 25     vjp, ans = _make_vjp(fun, x)
     26     if not vspace(ans).size == 1:
     27         raise TypeError("Grad only applies to real scalar-output functions. "

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/autograd/core.py in make_vjp(fun, x)
      8 def make_vjp(fun, x):
      9     start_node = VJPNode.new_root()
---> 10     end_value, end_node =  trace(start_node, fun, x)
     11     if end_node is None:
     12         def vjp(g): return vspace(x).zeros()

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/autograd/tracer.py in trace(start_node, fun, x)
      8     with trace_stack.new_trace() as t:
      9         start_box = new_box(x, t, start_node)
---> 10         end_box = fun(start_box)
     11         if isbox(end_box) and end_box._trace == start_box._trace:
     12             return end_box._value, end_box._node

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/autograd/wrap_util.py in unary_f(x)
     13                 else:
     14                     subargs = subvals(args, zip(argnum, x))
---> 15                 return fun(*subargs, **kwargs)
     16             if isinstance(argnum, int):
     17                 x = args[argnum]

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/blend.py in _loss(self, *parameters)
    117         parameter
    118         """
--> 119         model = self.get_model(*parameters)
    120         # Caculate the total loss function from all of the observations
    121         total_loss = 0

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/component.py in get_model(self, *params)
    533                 p = params[i : i + j]
    534                 i += j
--> 535                 model = model + c.get_model(*p)
    536         else:
    537             for c in self.components:

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/component.py in get_model(self, *parameters)
    329             morph = self.morph
    330         else:
--> 331             morph = self._func(*fparams)
    332             self._morph = morph._value
    333             morph = self._pad_morph(morph)

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/component.py in _func(self, *parameters)
    300 
    301     def _func(self, *parameters):
--> 302         return self.kwargs["func"](*parameters)
    303 
    304     def get_model(self, *parameters):

~/anaconda3/envs/scarlet37/lib/python3.7/site-packages/scarlet-1.0670feaa-py3.7-macosx-10.9-x86_64.egg/scarlet/source.py in _psf_wrapper(self, *parameters)
    394 
    395     def _psf_wrapper(self, *parameters):
--> 396         return self.frame.psf.__call__(*parameters, bbox=self.bbox)[0]
    397 
    398 

TypeError: 'NoneType' object is not subscriptable

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-18-767a06276956> in <module>
      1 blend = scarlet.Blend(sources, observation)
      2 get_ipython().run_line_magic('time', 'blend.fit(200)')
----> 3 print("scarlet ran for {0} iterations to logL = {1}".format(len(blend.loss), -blend.loss[-1]))
      4 plt.plot(-np.array(blend.loss))
      5 plt.xlabel('Iteration')

IndexError: list index out of range

[pmelchior]

That's maybe a little opaque but not a bug. We split the observed psfs into a minimal part and a difference kernel. The minimal part prevents undersampling. It's the effective PSF of the model, which we convolved with the diff kernel to match the observations. PointSource renders a shifted model_psf from its analytic form, a small Gaussian. That's what self.frame.psf.__call__ is supposed to do. You attempted to use the images of observed PSF instead. That fails because we don't include code to shift a PSF image to an arbitrary location.

[hcferguson]

Thanks @pmelchior. If I have an analytic ePSF that is separate for each channel, can I pass it a function that takes the channel as an argument, as well as x,y?

[fred3m]

That is possible, just not with the PointSource class included in scarlet. I think the best way to do this would be to initialize both your model frame and observation without passing either a PSF. That prevents the code from doing a convolution during fitting. For your point sources you'll need to either subclass PointSource or more likely scarlet.component.FunctionComponent and take ePSF as an argument. Then you can set your parameters (x, y, amplitude) similar to the way it is done in PointSource, where you can make your model as simple or complex as you like.

[hcferguson]

Thanks. I'll give that a try.

[pmelchior]

@hcferguson Can I ask you why you want a model PSF that is different in each band. We don't support this mode for a good reason, namely that you can't guarantee consistent colors that way.

[hcferguson]

I have created an oversampled ePSF "model" using the machinery in Photutils. For my current tests, it's not essential, I didn't instantly see that you could feed Scarlet an oversampled PSF, and in any case the photutils machinery already does the interpolation, so I figured I could just pickle the models.