LinearLSQFitter does not handle 2D weights properly

When passing 2D weigths to LinearLSQFitter in 2D mode:

Traceback (most recent call last):
  [....]
  File "/Users/michal/sources/astropy/astropy/modeling/fitting.py", line 456, in __call__
    lhs *= weights[:, np.newaxis]
ValueError: operands could not be broadcast together with shapes (441,10) (21,1,21) (441,10)

LevMarLSQFitter works properly.

Test: https://gist.github.com/majkelx/0f105055554ce5061264f18be44b5171

Thanks. I was going to report this but you've saved me the bother. Here's a slightly simpler example that doesn't work FWIW:

import numpy as np
from astropy.modeling import models, fitting

fitter = fitting.LinearLSQFitter()

y, x = np.mgrid[0:5,0:5]
z = x+y

weights = np.ones_like(x)

model = models.Polynomial2D(0)

fit = fitter(model, x, y, z, weights=weights)

 [...]
  File "/home/jturner/repos/astropy/astropy/modeling/fitting.py", line 475, in __call__
    lhs *= weights[:, np.newaxis]
ValueError: non-broadcastable output operand with shape (25,1) doesn't match the broadcast shape (5,25,5)

Related problem: This should probably be fixed at the same time as 2D weights, since they are both a matter of reworking the same logic. LinearLSQFitter expects the dimensionality of weights to match the x array (like the other fitters), but this doesn't allow for model sets to have different weights for different models. An (ideally-optional) model_set_axis dimension should probably be accepted in the argument to weights, since the correct weighting could depend on different data or noise values in different models in the set. There may be some cases where a common set of weights can be used for all this models; this could probably be inferred from the array dimensionality.