Fitting of CompoundModel is slow

[saimn]

I'm playing with some photometry computation with photutils, which is horribly slow, and this appears to be because of astropy.modeling trying to know if the model has some units. From what I understand with the traceback below (after interrupting the process), accessing self.input_units triggers a recursive lookup that takes a lot of time, exponentially increasing with the number of sources in the CompoundModel. And this for a model that has no units, and for each iteration if I get it correctly.

https://github.com/astropy/astropy/blob/master/astropy/modeling/core.py#L2861

I cannot provide reliable timing for now as it takes just too long, but I measured one call to prepare_units to 17sec, with a model with 24 (gaussian) sources. Replacing ._validate_input_units by a no-op method makes thing work, fitting the 24 sources group in 5 seconds (I will try to get the timing without this tomorrow).

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-43-b96405ebf871> in <module>
----> 1 photometry_results = basic_photometry(subim)
      2 photometry_results

~/dev/photutils/photutils/psf/photometry.py in __call__(self, image, init_guesses)
    185         including the `__call__` signature.
    186         """
--> 187 
    188         return self.do_photometry(image, init_guesses)
    189 

~/dev/photutils/photutils/psf/photometry.py in do_photometry(self, image, init_guesses)
    288                                     [getattr(self.psf_model, param).value])
    289 
--> 290         star_groups = self.group_maker(init_guesses)
    291         output_tab, self._residual_image = self.nstar(image, star_groups)
    292 

~/dev/photutils/photutils/psf/photometry.py in nstar(self, image, star_groups)
    348                                         mode='trim')[0]] = True
    349 
--> 350             fit_model = self.fitter(group_psf, x[usepixel], y[usepixel],
    351                                     image[usepixel])
    352             param_table = self._model_params2table(fit_model,

~/dev/astropy/astropy/modeling/fitting.py in wrapper(self, model, x, y, z, **kwargs)
    176         else:
    177 
--> 178             return func(self, model, x, y, z=z, **kwargs)
    179 
    180     return wrapper

~/dev/astropy/astropy/modeling/fitting.py in __call__(self, model, x, y, z, weights, maxiter, acc, epsilon, estimate_jacobian)
    885             self.objective_function, init_values, args=farg, Dfun=dfunc,
    886             col_deriv=model_copy.col_fit_deriv, maxfev=maxiter, epsfcn=epsilon,
--> 887             xtol=acc, full_output=True)
    888         _fitter_to_model_params(model_copy, fitparams)
    889         self.fit_info.update(dinfo)

~/.pyenv/versions/3.7.0/lib/python3.7/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    392             maxfev = 200*(n + 1)
    393         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 394                                  gtol, maxfev, epsfcn, factor, diag)
    395     else:
    396         if col_deriv:

~/dev/astropy/astropy/modeling/fitting.py in objective_function(self, fps, *args)
    823         meas = args[-1]
    824         if weights is None:
--> 825             return np.ravel(model(*args[2: -1]) - meas)
    826         else:
    827             return np.ravel(weights * (model(*args[2: -1]) - meas))

~/dev/astropy/astropy/modeling/core.py in __call__(self, x, y, model_set_axis, with_bounding_box, fill_value, equivalencies)
    380                                      ('with_bounding_box', False),
    381                                      ('fill_value', np.nan),
--> 382                                      ('equivalencies', None)])
    383 
    384             # The following makes it look like __call__ was defined in the class

~/dev/astropy/astropy/modeling/core.py in __call__(self, *inputs, **kwargs)
    361             def __call__(self, *inputs, **kwargs):
    362                 """Evaluate this model on the supplied inputs."""
--> 363                 return super(cls, self).__call__(*inputs, **kwargs)
    364 
    365             # When called, models can take two optional keyword arguments:

~/dev/astropy/astropy/modeling/core.py in __call__(self, *inputs, **kwargs)
    755         that were specified when the model was instantiated.
    756         """
--> 757         inputs, format_info = self.prepare_inputs(*inputs, **kwargs)
    758 
    759         parameters = self._param_sets(raw=True, units=True)

~/dev/astropy/astropy/modeling/core.py in prepare_inputs(self, model_set_axis, equivalencies, *inputs, **kwargs)
   1477                                model_set_axis, self.standard_broadcasting)
   1478 
-> 1479         inputs = self._validate_input_units(inputs, equivalencies)
   1480 
   1481         # The input formatting required for single models versus a multiple

~/dev/astropy/astropy/modeling/core.py in _validate_input_units(self, inputs, equivalencies)
   1495         # Check that the units are correct, if applicable
   1496 
-> 1497         if self.input_units is not None:
   1498 
   1499             # We combine any instance-level input equivalencies with user

~/dev/astropy/astropy/modeling/core.py in input_units(self)
   2860     @property
   2861     def input_units(self):
-> 2862         return self._generate_input_output_units_dict(self._tree.inputs_map, 'input_units')
   2863 
   2864     @property

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

~/dev/astropy/astropy/modeling/utils.py in inputs_map(self)
     83         else:
     84             for inp in self.left.inputs:
---> 85                 m, inp2 = self._recursive_lookup(self.left, self.left.inputs_map, inp)
     86                 inputs_map[inp] = m, inp2
     87 

KeyboardInterrupt: 

[pllim]

cc @nden @perrygreenfield

[saimn]

So, the model with 24 sources took 2210 sec., compared to 3.45 sec. without the units check.

[perrygreenfield]

Can you supply the model you used? I've been looking at this particular call very recently in my rework of modeling, just out of chance. There isn't that much to it so I'm surprised that recursive call takes that much time, but I'd like to look into it.

[saimn]

The model is created by photutils, it is a sum (https://github.com/astropy/photutils/blob/master/photutils/psf/utils.py#L109) of IntegratedGaussianPRF (https://github.com/astropy/photutils/blob/master/photutils/psf/models.py#L795). I wanted to try if I can reproduce with a sum of astropy Gaussian2D, but did not have the time yet.

[saimn]

Here is a reproduction directly with astropy, just adding Gaussian2D models: https://gist.github.com/saimn/cea2d1dff0ed933525efb743553a80ad
36 sec. for 25 gaussians, 1324 sec. for 30 ! Just for a .render().

[saimn]

Ping @larrybradley in case you missed this, since this makes photutils' psf photometry nearly impossible to use, unless stars have no close neighbors ;).

[perrygreenfield]

I just confirmed the slowness with just Gaussian2D models a few minutes before you posted that. If it is any consolation, the new version of modeling (won't be out until the next major version of astropy unfortunately) doesn't suffer this problem. The time for 25 Gaussian2D components in that only takes 4.5 ms. I'm still waiting for the other to finish...

[perrygreenfield]

Just finished. 36 seconds (I used timeit so it ran it 7 times). What, a factor of 8,000 is unacceptable?

[perrygreenfield]

I'll see if I can understand why the current one behaves this way and if there is any simple fix.

[larrybradley]

@saimn Yes, I noticed. :) I'm hope the new modeling updates will be a big improvement.

[larrybradley]

Fitting ePSFs to ~1000 sources took me about 20 minutes. A factor of 8000 improvement would bring that down to 0.15 seconds, which is certainly acceptable. 😄

[saimn]

Great for the new version!
Meanwhile it would be great if we can find a fix, there must be another way to know if a model uses unit or not ?

[perrygreenfield]

The problem was tracked down to util.py ExpressionTree property inputs_map. This only shows up when there are more than one input coordinate and results in a doubling of calls to this property at every level of the expression tree. At 24 levels that basically becomes something on the order of 2**24 calls, a moderately substantial number, even for a simple routine ;-) I have a fix for this that I will later submit. I'll need to test it more thoroughly.

[nden]

@perrygreenfield I assume you didn't mean to close the issue.

[perrygreenfield]

Right. They should make that button harder to get to :-)