Models to include analytic integral, if possible

[pllim]

I propose a new method property for Astropy model that returns the model's analytic integral as another model (if available) or None (if not). Pseudo-code as follows:

In core.py:

class _ModelMeta(...):
    @property
    def integral(self):
        return None

In the module of relevant model(s):

class SomeModel1D(Fittable1DModel):
    @property
    def integral(self):
        return ErrorFunction(...)

class SomeOtherModel(...):
    @property
    def integral(self):
        return ErfModel(...)

Alternatively, models can also fallback to some kind of basic integration (e.g., trapezoid integration) if there is no analytic one. But that can be a separate discussion. Here is an example by @mhvk on how to do such a thing:

if model.integral:
    result = model.integral(stop) - model.integral(start)
else:
    # code to do trapezoidal integration
    ....

[mcara]

@pllim @perrygreenfield While I do not know the specific needs that prompted this issue, It seems that adding analytical integration to models to be the beginning of an unnecessary over-complication: why stop at integration? Why not add differentiation, Jacobians, Laplacians, gradients, etc.?

However, it does appear that integration is a special kind of operation used quite often. So, let's say models will get this kind of method. However, I wonder, what would be the integral of the Shift and Scale or, even more exciting, of the Rotation2D or RotateNative2Celestial models? What would be the meaning of these integrals? What about analytic integrals of compound models? Are we going to implement some sort of symbolic integration algorithms?

[mcara]

@pllim @perrygreenfield Another issue is that in your example you havestart and stop arguments. This may be fine for 1D case but for 2D models one should be able to provide some complex domains of integration (not just rectangular) if these methods are to be useful for a wider set of problems.

[mcara]

I would think that it may be much more powerful to return an "integrated model" instead of a number. For example, Gaussian1D.integrated_model~Erf_model? But then, we need to implement models for Erf, and its integral, and the integral of the integral of Erf, etc.

[pllim]

The motivation on my part is to support flux integration in pysynphot, which currently only has trapezoid integration (no matter what the model is). Certainly not all models have analytic integral, which is why it returns None in the base model, so that models that don't have that defined would just return a None. Also, @perrygreenfield mentioned some time ago via private communications that the analytic integral argument(s) can be different for different models.

My proposal is very basic, and I certainly welcome other ideas. integrated_model sounds interesting... Can every analytic integral be represented as a model? Perhaps yes...

[mhvk]

@mcara - The derivative of a model would be very valuable for model fitting, so that would certainly be useful! Do agree that it needs some thought for composite models, though at least in principle it is not difficult (but yes I do worry about reinventing symbolic algebra...).

Like @pllim, I rather like your idea of any such method to return a new model instance; indeed, that should avoid a lot of duplicate code (e.g., suppose we did have an ErrorFunction model (may we do, I didn't check!), then it one wanted its derivatives for fitting one would need to add a method that calculates gaussians, which is silly to re-implement if we already have a Gaussian model. Of course, one would need to be very careful with proper conversion of parameters.

For this specific purpose, I also like that returning the integral as a model removes the need for the modeling class to care about how these integrals are being used, i.e., pysynphot can integrate between some start and stop itself knowing the number of dimensions, etc., but the code doesn't have to care.

Finally, maybe good to note that this echoes what is already done for inverse (though I was a bit puzzled to see that e.g., in EulerAngleRotation it is a method whereas in all other classes in rotations it is a property... ... but looking at the core, I see "We allow the @property decorator to be omitted entirely from the class definition, though its use should be encouraged for clarity" -- not sure I think that is a great idea...).

By analogy with inverse, I think one can keep the name simple: derivative and integral would seem fine to me.

[mhvk]

@mhvk I think you just added @property to the list of recipients of your message 😆

Which does turn out to be a user on github.... Though not one that has contributed anything.

[pllim]

I think Erik B. must have done some magic to the model properties. Because even if inverse is not explicitly stated to be a property in subclass, it is actually still a property.

>>> from astropy.modeling.models import EulerAngleRotation
>>> m = EulerAngleRotation(0, 0, 0, ['x', 'y', 'z'])
>>> m.inverse
<EulerAngleRotation(phi=-0.0, theta=-0.0, psi=-0.0)>

I think so far we agree that @mcara's suggestion to have the integral itself to be a model is a good one. Therefore, I'll edit my proposal above to indicate it as such. But I still prefer to call it analytic_integral, since we have inverse and not inverse_model.

[pllim]

However... The question that remains is if analytic_integral returns a model, how do I actually integrate it, say, to get integrated flux?

[mhvk]

If we are returning a model like inverse, then isn't it OK to just call the properties integral and derivative? (i.e., why add analytic -- by analogy with inverse, I think this should be the default; we can worry about other ways to do it later.).

As for evaluating, isn't the idea that you would simply do

if model.integral:
    result = model.integral(stop) - model.integral(start)
else:
    # code to do trapezoidal integration
    ....

[pllim]

Yeah, I am okay with just integral. Will update the proposal again.

[mcara]

I am inclined to think that analytic be dropped. On the other hand, if integration is performed numerically (when analytically is not feasible), it would be useful to specify accuracy, method, etc. Also, if model integrals are evaluated numerically, then for computing definite integral such as in @mhvk example above (result = model.integral(stop) - model.integral(start)), things are inefficient.

Also, for 1D case @mhvk example is what needs to be done. For a 2D case - I am not sure it can de done in a simple way (boundaries, order of integration issues, etc.)

So, after @pllim question, I am not sure if this approach of returning an "integrated model" is even doable in general [Edit: Maybe it is doable but would it be that useful?]

For particular needs of @pllim maybe it would make more sense to have a method/@property that returns "integrated flux" with some smart way of providing boundaries for >=2D cases.

[pllim]

Let's start with something simple like Box1D, the integrated result (which I think is the same as its area) is amplitude * width. But what would be its integral model? And how can such a model be used to evaluate something equivalent to amplitude * width?

[mhvk]

The integral would be a piecewise model, 0 before the start, linear in slope over the box, and constant thereafter.

[pllim]

Looking at this particular case for Box1D, having integral() to be a method that returns the integrated result (not the model) might be more desirable. I mean, simply multiplying the height and width is a lot more efficient here. Don't you think?

[mhvk]

@pllim - but that is only useful if your start and stop are on both sides of the box. What if one of them is in the middle?

[mcara]

What about having two methods/properties such as [rename as you wish]: integrated_flux(boundaries), antiderivative_model (for 1D models)? It seems that the integrated_flux is more urgent than antiderivative.

[pllim]

Along the line of @mcara's thought above... Maybe something like:

class SomeModel(...):
    @property
    def integral(self):
        return SomeOtherModel(...)

    def integrate(self, start=None, stop=None, ...):
        # Here we can opt to use self.integral or do simple math
        integrated_result = do_something
        return integrated_result

For the Box1D, it would be like:

class Box1D(...):
    @property
    def integral(self):
        return SomeFancyModel(...)

    # Here, the actual integration does not use the model
    def integrate(self, start=None, stop=None):
        x1, x2 = self.bounding_box

        if start is None:
            start = x1
        if stop is None:
            stop = x2
        if stop < start:
            raise ValueError(...)  # Or we can be fancy here and switch them?

        w = stop - start
        return w * self.amplitude

EDIT: Corrected some typo.

[mcara]

@pllim 👍 for succinct naming (Edit: I am not being sarcastic)

[mcara]

If more "differentiation" between attribute names is desired, the method that returns a model (in @pllim example - integral) could be named (see https://en.wikipedia.org/wiki/Antiderivative) primitive, primitive_model, antiderivative, antiderivative_model, integral_model [again, I do not have strong preferences and these are just suggestions].

[mhvk]

The integral is not particularly fancy: it just needs a new Ramp1D model, with an amplitude, center-of-ramp, and a width:

class Ramp1D(...):
    @staticmethod
    def evaluate(x, amplitude, x_0, width):
        return amplitude * np.clip((x - (x_0 - width / 2.)) / width, 0., 1.)

(Of course, in evaluating Box1D.integral, it should set the amplitude to amplitude * width.)

And now just do

integrated_flux = box1d.integral(stop) - box1d.integral(start)

p.s. Your example is not quite right: you need to check for the case where start < x1 and/or stop > x2.

[pllim]

@mhvk , but performance-wise, using a model in place of simple math would incur unnecessary overhead?

[mhvk]

My main worry is that for integrate, I'm not sure the proposed API will be sufficient (e.g., the 2-D case mentioned by @mcara; how to deal width gridding in non-analytically integrable functions; etc.), while I'm fairly confident that an integral property will be useful no matter what. So, my proposal would be to add that first, and see if it does in fact help you enough in pysynphot.

As for overhead: it depends a bit on what you need. For modelling spectra, e.g., one would have a whole slew of pixels to integrate over, and for that case I would guess the overhead of creating the new model instance would be not very relevant. Anyway, all of this may be premature optimization...

p.s. But just for fun, since evaluate is a staticmethod, one could quite easily do something like:

class Box1D(...):
    def integrate(...):
        ramp_amp = amplitude * width
        return Ramp1D.evaluate(stop, ramp_amp, x0, width) - Ramp1D.evaluate(start, ramp_amp, x0, width)

[pllim]

Okay, putting my boxy concerns aside, there are some things that need some thinking, which are voiced above but I am putting them in a list here so we don't get lost in the long thread:

1. Is it even possible to define a integral model for 2D cases? If not, are we going to just focus on 1D and ignore 2D?
2. Should we just implement integrate for now (returns the integrated flux, not model) and worry about model later? Also, my focus right now is 1D, so some thinking need to be done for 2D as well.
3. Should we rename integral and integrate?

Did I miss anything?

[mhvk]

Your comment about 2-D is well taken. In the analytic sense, the only way I can see to solve it would be to pass on which axis the integral would be over, which would then return a model that would need to be integrated over the other axis.

How an API that makes this possible in principle and also solves your wish:

def get_integral(self, start=None, stop=None, over='x'):
    """Return integral over ``over``.

    If start and stop are not given, the indefinite integral over ``over`` is returned,
    otherwise the integral is evaluated over the range ``start`` to ``stop``.
    """
    if over == 'x':
        integral_model = IntegralOverXModel(...)
    elif over == 'y':
        ...
    if start is None and stop is None:
        return integral_model
    elif start is not None and stop is not None:
        return integral_model(stop) - integral_model(start)
    else:
        raise ValueError("Start and stop should either both be absent or both given")

integral = property(get_integral)

For 2-D integral_model(stop) - integral_model(start) would need to return a new model, which would depend only on y (and parameters partially determined by start and stop).

If one gets really fancy this could even be extended to integrating over parameters rather than an axis... Also, for different models one might need to special case, e.g., start=-np.inf, etc.

[pllim]

Thanks, @mhvk and @mcara ! I should turn this into a prototype PR when I get the chance (unfortunately, not for another week or two).

[pllim]

See #5108 for a possible implementation. I thought about adding start and stop but decided not to, as it overly complicate things.

[pllim]

With the 4.0 refactoring in modeling and no clear consensus for support beyond 1D, I am giving up on a generic integral support here and probably just going to roll my own in my package as needed.

[pllim]

FWIW, I implemented integrals for select 1D models over at spacetelescope/synphot_refactor#252 if there is interest to revisit this in the future. I did not end up implementing them as models themselves to not over-complicate things for my simple use cases.

[mcara]

I think this is doable for the 1D case. If you had the following inheritance: Model->Model1D (and Model2D,...)->Fittable1DModel->SpecificModel, then integration interface could be made available in Model1D. Higher order models would not have methods such as integrate. Model1D would be just like the Fittable1DModel.

[pllim]

Feel free to open a PR over at synphot_refactor if you really want them as models. I think I am done with this for a while...