Implement bounding_box model property and render_model function - #3909

I will review this later today.

@patti – Thanks! Can't wait to use this for my analyses!

I'll let @adonath do the review of the API / implementation.
One thing you could do is translate your example notebook into tests and docs.

@embray @adonath For high-level docs, what do you think is best?

• Briefly mention this in http://astropy.readthedocs.org/en/latest/modeling/index.html#getting-started ?
• A separate page in http://astropy.readthedocs.org/en/latest/modeling/index.html#using-astropy-modeling with explanations and a 1d / 2d example?
• Or something else, e.g. have some high-level docs in the nddata section?

@patti – Sorry, I forgot to ask you explicitly ... what do you think would be best for the high-level docs for this PR (ideally so that they can grow to cover what you'll implement in the coming weeks)?

Concerning tests: start with one simple test that proves that this works at all (once it does). Then maybe parametrise it and include cases where the source is very small and very large and overlaps with the image boundaries. Gaussians and box models are great for tests, because they have analytical solutions on how much flux is within rectangular boxes, i.e. it's possible to verify the reference values.

Thanks for working on this! I might make sense to add 'WIP: ' to the title of the PR, just to make clear that it's work in progress. And could you add a to do list in your first post? We can work that list out later during the hangout...

I'll review all this in a bit, but I'll just add that it's fine if it's only working, for now, on a handful of models. I think it's expected that these features will be unimplemented by default, and would be implemented on individual models on a case-by-case basis.

Thanks for providing the notebook! I hope large parts of it will transform into nice plots and usage examples for the docs:-) I noticed, that the difference image doesn't look correct. It seems like the positions with and without the bounding box are shifted with respect to each other. We'll have to check that in detail. Could you add a very simple test, that checks that the evaluation with and without the bounding box gives the same result for let's say a Gaussian or a box? At best this is checked for a few
positions and shapes of the bbox.

@embray In particular I'd be interesting to know, how you think the implementation of the bounding box feature would fit best in the current framework, concerning e.g. n-dimensional models and compound models.

@adonath regarding "It seems like the positions with and without the bounding box are shifted with respect to each other."

It just looks that way because I didn't add 1 to the upper bounds to make it inclusive. I'll fix that.

I'm having a hard time trying to diagnose a bug where I get a row of zeros in the Gaussian2D bounded image for certain values of theta. It's reproduced in a notebook, here. Any ideas?

I'll go over this in more detail in a sec, but there's something in your example notebook that got me thinking, where you created parameters for 1000 sources, and then created 1000 Gaussian2D instances to model them:

imshape = (800,1250)
sky_image = np.ones(imshape)
nsrc = 1000
data = {}
data['x_mean'] = np.random.rand(nsrc)*imshape[1]
data['y_mean'] = np.random.rand(nsrc)*imshape[0]
data['r_amp'] = np.abs(np.random.randn(nsrc)+10.)
data['g_amp'] = np.abs(np.random.randn(nsrc)+10.)
data['b_amp'] = np.abs(np.random.randn(nsrc)+10.)
data['x_stddev'] = np.abs(np.random.rand(nsrc)+3.)
data['y_stddev'] = data['x_stddev']*(np.random.rand(nsrc)+1)
data['theta'] = np.random.rand(nsrc)*np.pi*2.


mod_list = [Gaussian2D(x_mean=data['x_mean'][i], y_mean=data['y_mean'][i], 
                        x_stddev=data['x_stddev'][i], y_stddev=data['y_stddev'][i],
                        amplitude=data['g_amp'][i], theta=data['theta'][i])
            for i in range(nsrc)]

Ideally this should be unnecessary. The Model interface has a useful, but poorly understood feature called "model sets" (introduced here and described in more detail here. It is for exactly this purpose--when you have many different parameterizations of the same model. You could make a model set of 1000 Gaussian2Ds like:

models = Gaussian2D(x_mean=data['x_mean'], y_mean=data['y_mean'], 
                                     x_stddev=data['x_stddev'], y_stddev=data['y_stddev'], 
                                     amplitude=data['g_amp'], theta=data['theta'], n_models=nsrc)

Then you can evaluate all of these models simultaneously against the image like:

sources = models(x, y, model_set_axis=False)

where model_set_axis=False can be read as "evaluate each model in the set on the same input".

You can then combine them into a single image just be calling image = np.sum(sources, axis=0). Unfortunately, for a case like this, this is extremely prohibitive on a typical desktop system, because you end up with 1000 800x1250 images, one for each source. However, if each model has a bounding box, we could change the behavior of model set evaluation to evaluate each model on its bounding box (in this case the .bbox property would probably return an array, with one entry for each model in the set). Since each bbox isn't the same size I guess it would end up having to be an array of different-sized arrays, but that would still be a lot better. Then one could build up the image like:

image = np.zeros(imgsize, dtype='float')
sources = models(x, y, model_set_axis=False)  # Returns bounded results
for bbox, source in zip(models.bbox, sources):
    image[bbox] += source

I think being able to do something akin to that should be able to work, but what do you think? (I don't personally claim to know exactly what the priorities are for this work myself :)

I'm having a hard time trying to diagnose a bug where I get a row of zeros in the Gaussian2D bounded image

I haven't fully investigated this yet, but I think there are still some rounding issues with the bounding boxes. They are not correctly rounded up or down. For example, when you call np.meshgrid you just call it with arguments like: np.arange(bbox[0],bbox[1]).astype(int). But astype(int) effectively takes the floor, where in this case it looks like you might be wanting to round to the closest integer instead.

I also noticed you were using np.meshgrid directly, but you can more simply use

y, x = np.mgrid[bbox[2]:bbox[3], bbox[0]:bbox[1]]

and similarly use the bbox itself and slicing to index the image, instead of using large array indices.
In fact, it might be nice if the bbox attribute just returned a pair of slice objects instead of a 4-tuple, so it can be used to directly slice an array.

@embray Thanks for looking into this! Confirmed on the rounding issue. The plan is to use extract_array and add_array in nddata.utils to save myself some hassle down the line anyway. I just haven't got that working quite yet.

Just to make sure I understand what you mean by a slice object, this would look like a tuple of arrays bbox = (x_slice, y_slice) corresponding to the coordinate indices of the bbox? Effectively, y_slice, x_slice = np.mgrid[bbox[2]:bbox[3], bbox[0]:bbox[1]], or do you mean something different?

Regarding model sets, this was discussed and there were reasons we opted against using them that I don't recall exactly. @adonath or @cdeil?

In general, I agree it would make sense to implement this in model.call. @adonath didn't think it belonged in the super class- maybe because it's just more involved? Would this call some version of render_models2sky() that does live in the modeling.utils?

@embray @nden -- friendly bump on this...any chance you could approve and merge in the next few days before @patti's GSOC ends? I'd offer to review but know nothing about the modeling package :)

Sorry, I haven't been following this. I'll try to review it very soon (unless @embray beats me to it).

I can use this kind of functionality in my work as well. A natural extension would be to have this integrated in the __call__ method, and also working with fitters and compound models. I guess for compound models the default would be the intersection of the all default bbox ranges but something similar to the custom inverse would be useful as well. This is a very good start on this topic. Thanks @patti !

@nden – Thank you!

Sorry, I fell behind on being able to review this--I think would have preferred not to merge as-is, unless there was some agreement to make changes in a follow-up PR. But I'll try to catch up on what the latest conversation was.

Update: Ah, I understand now GSoC was coming to an end and we wanted to get this merged first. In that case that's fine, sorry for whining. And sorry especially for falling behind on code review.

(Also please don't merge PRs without first assigning a milestone).

Some nitpicks above, but otherwise this is looking great. Love the example image.
Sorry I didn't get to review sooner. I'll try to catch up on the past discussion now.

@adonath @patti Thanks for the explanation of the argument against using slices. This makes sense to me and is nice and more generalizable.

I do think it would be nice to be able to get the bounding box as a slice/tuple of slices where applicable though (such as when the coordinates are pixel coordinates). This, in addition to some other issues, is leading me to think that there should be some kind of tuple subclass for representing bounding boxes. This subclass could include a .slice attribute or something of the sort that automatically converts to a slice when possible. Users wouldn't have to use this subclass directly since the implementation can automatically convert user-supplied tuples to this special class.

Not a high priority but something to consider.

@nden

I can use this kind of functionality in my work as well. A natural extension would be to have this integrated in the __call__ method

This was my first thought as well. Though after some (much) earlier discussion (now buried by the terrible GitHub discussion interface), there was some agreement that maybe it would be too confusing to lump this functionality into the normal __call__ interface. Instead the render_model method would be added to models for this purpose.

I could go either way to be honest, but there were some good points (that don't recall exactly) made in favor of using a separate method. (Edit: Ah right, that discussion was here: #3909 (diff))

I guess for compound models the default would be the intersection of the all default bbox ranges but something similar to the custom inverse would be useful as well.

Yes, that's what I was thinking as well. I think for a compound model the default_bounding_box would be the intersection. But @patti's API should also already allow for setting a custom bounding box on a compound model.

@embray Thanks for the review!
@patti I think it's best to address these comments in the follow up PR #4040

What @adonath said. I could probably try to address a few of my own comments too if it will help move things along (since most of them were minor nitpicks).