astropy: FittingWithOutlierRemoval doesn't work for model sets

Attempting to use FittingWithOutlierRemoval in place of a simple LinearLSQFitter instance produces a ValueError: Input argument u'x' does not have the correct dimensions in model_set_axis=0 for a model set with n_models=4176. This blocks the very common use case of fitting each row of an image with rejection, as in IRAF’s fit1d. One could, of course, fit each row with a separate model instance, but that would be much slower – and even with model sets, fitting already seems unpleasantly slow compared with IRAF, so it probably wouldn’t be viable.

Unfortunately, I think fixing this might be non-trivial: Within FittingWithOutlierRemoval, L557 (also L567-568) selects only good points from an input array of x co-ordinates, like so:

https://github.com/astropy/astropy/blob/667b7ab621dfa80290fbe36cd711ad0a427ca272/astropy/modeling/fitting.py#L556-L559

Obviously the rejected points will differ from one row/model to another, which means that FittingWithOutlierRemoval cannot simply pass the underlying Fitter a single array of x co-ordinates (in N-1 dimensions) to use for every model. Unlike model evaluation, I think fitting of model sets in general only works with a common x array for all the models, so LinearLSQFitter and friends may need substantial changes for this to work. Moreover, I believe the matrix equation solved by the underlying np.linalg routine assumes a 2D Vandermonde matrix, corresponding to a single set of points in x – I’m not even sure (off the top of my head) that the problem can be expressed in a form that the back-end solver can work with, which would be rather unfortunate.

Would you mind confirming that this is a correct analysis, @nden (assuming you know off hand)? I’ve always had a bit of a hard time discerning from the documentation exactly what dimensionality of inputs the models and fitters accept under different circumstances.