@@ -392,6 +392,19 @@ API changes
     with True only if the comparison was true for all elements compared,
     which could lead to confusing circumstances. [#3912]
 
+  - Using ``model.inverse = None`` to reset a model's inverse to its
+    default is deprecated.  In the future this syntax will explicitly make
+    a model not have an inverse (even if it has a default).  Instead, use
+    ``del model.inverse`` to reset a model's inverse to its default (if it
+    has a default, otherwise this just deletes any custom inverse that has
+    been assigned to the model and is still equivalent to setting
+    ``model.inverse = None``). [#4236]
+
+  - Adds a ``model.has_user_inverse`` attribute which indicates whether or not
+    a user has assigned a custom inverse to ``model.inverse``.  This is just
+    for informational purposes, for example, for software that introspects
+    model objects. [#4236]
+
   - Renamed the parameters of ``RotateNative2Celestial`` and
     ``RotateCelestial2Native`` from ``phi``, ``theta``, ``psi`` to
     ``lon``, ``lat`` and ``lon_pole``. [#3578]
 
@@ -246,33 +246,22 @@ def rename(cls, name):
 
     def _create_inverse_property(cls, members):
         inverse = members.get('inverse')
-        if inverse is None:
+        if inverse is None or cls.__bases__[0] is object:
+            # The latter clause is the prevent the below code from running on
+            # the Model base class, which implements the default getter and
+            # setter for .inverse
             return
 
         if isinstance(inverse, property):
-            fget = inverse.fget
-        else:
             # We allow the @property decorator to be omitted entirely from
             # the class definition, though its use should be encouraged for
             # clarity
-            fget = inverse
-
-        def wrapped_fget(self):
-            if self._custom_inverse is not None:
-                return self._custom_inverse
-
-            return fget(self)
-
-        def fset(self, value):
-            if not isinstance(value, (Model, type(None))):
-                raise ValueError(
-                    "The ``inverse`` attribute may be assigned a `Model` "
-                    "instance or `None` (where `None` restores the default "
-                    "inverse for this model if one is defined.")
-
-            self._custom_inverse = value
+            inverse = inverse.fget
 
-        cls.inverse = property(wrapped_fget, fset, doc=inverse.__doc__)
+        # Store the inverse getter internally, then delete the given .inverse
+        # attribute so that cls.inverse resolves to Model.inverse instead
+        cls._inverse = inverse
+        del cls.inverse
 
     def _handle_backwards_compat(cls, name, members):
         # Backwards compatibility check for 'eval' -> 'evaluate'
@@ -579,7 +568,8 @@ class Model(object):
     # optionally; in that case it is of course up to the user to determine
     # whether their inverse is *actually* an inverse to the model they assign
     # it to.
-    _custom_inverse = None
+    _inverse = None
+    _user_inverse = None
 
     # If a bounding_box_default function is defined in the model,
     # then the _bounding_box attribute should be set to 'auto' in the model.
@@ -777,8 +767,6 @@ def ineqcons(self):
 
         return self._constraints['ineqcons']
 
-    # *** Public methods ***
-
     @property
     def inverse(self):
         """
@@ -791,16 +779,73 @@ def inverse(self):
         `~astropy.modeling.polynomial.PolynomialModel`, but not by
         requirement).
 
+        A custom inverse can be deleted with ``del model.inverse``.  In this
+        case the model's inverse is reset to its default, if a default exists
+        (otherwise the default is to raise `NotImplementedError`).
+
         Note to authors of `~astropy.modeling.Model` subclasses:  To define an
         inverse for a model simply override this property to return the
         appropriate model representing the inverse.  The machinery that will
         make the inverse manually-overridable is added automatically by the
         base class.
         """
 
+        if self._user_inverse is not None:
+            return self._user_inverse
+        elif self._inverse is not None:
+            return self._inverse()
+
         raise NotImplementedError("An analytical inverse transform has not "
                                   "been implemented for this model.")
 
+    @inverse.setter
+    def inverse(self, value):
+        if not isinstance(value, (Model, type(None))):
+            raise ValueError(
+                "The ``inverse`` attribute may be assigned a `Model` "
+                "instance or `None` (where `None` restores the default "
+                "inverse for this model if one is defined.")
+
+        if value is None:
+            warnings.warn(
+                "Currently setting `model.inverse = None` resets the inverse "
+                "to the default inverse (if one exists).  However, starting "
+                "in Astropy 1.2, setting `model.inverse = None` explicitly "
+                "forces a model to have no inverse (such that accessing "
+                "`model.inverse` raises a NotImplementedError) even if that "
+                "model's class has a default inverse.\n\n"
+                "Instead, call `del model.inverse` to reset the inverse to "
+                "its default (if a default exists for that model's class--"
+                "otherwise the model is reset to having no inverse.",
+                AstropyDeprecationWarning)
+
+        self._user_inverse = value
+
+    @inverse.deleter
+    def inverse(self):
+        """
+        Resets the model's inverse to its default (if one exists, otherwise
+        the model will have no inverse).
+        """
+
+        del self._user_inverse
+
+    @property
+    def has_user_inverse(self):
+        """
+        A flag indicating whether or not a custom inverse model has been
+        assigned to this model by a user, via assignment to ``model.inverse``.
+        """
+
+        return self._user_inverse is not None
+
+    @property
+    def _custom_inverse(self):
+        # Deprecated alias for _user_inverse--included solely for temporary
+        # backwards compatibility for pyasdf--remove after Astropy v1.1
+        # release.
+        return self._user_inverse
+
     @property
     def bounding_box(self):
         """
@@ -889,6 +934,8 @@ def bounding_box(self, limits):
 
         self._bounding_box = limits
 
+    # *** Public methods ***
+
     @abc.abstractmethod
     def evaluate(self, *args, **kwargs):
         """Evaluate the model on some input variables."""
@@ -1832,7 +1879,7 @@ class members that should be added to the generated
             # intermediate compound models are stored in the tree as well then
             # we can immediately check for custom inverses on sub-models when
             # computing the inverse
-            instance._custom_inverse = mcls._make_custom_inverse(
+            instance._user_inverse = mcls._make_user_inverse(
                     operator, left, right)
 
             return instance
@@ -1851,7 +1898,7 @@ def _handle_backwards_compat(mcls, name, members):
         return
 
     @classmethod
-    def _make_custom_inverse(mcls, operator, left, right):
+    def _make_user_inverse(mcls, operator, left, right):
         """
         Generates an inverse `Model` for this `_CompoundModel` when either
         model in the operation has a *custom inverse* that was manually
@@ -1873,7 +1920,7 @@ def _make_custom_inverse(mcls, operator, left, right):
         """
 
         if not (operator in ('&', '|') and
-                (left._custom_inverse or right._custom_inverse)):
+                (left._user_inverse or right._user_inverse)):
             # These are the only operators that support an inverse right now
             return None
 
 
@@ -688,6 +688,12 @@ def fit_deriv(x, slope, intercept):
         d_intercept = np.ones_like(x)
         return [d_slope, d_intercept]
 
+    @property
+    def inverse(self):
+        new_slope = self.slope ** -1
+        new_intercept = -self.intercept / self.slope
+        return self.__class__(slope=new_slope, intercept=new_intercept)
+
 
 class Lorentz1D(Fittable1DModel):
     """
 
@@ -3,14 +3,15 @@
 from __future__ import (absolute_import, unicode_literals, division,
                         print_function)
 
-import pytest
 import numpy as np
 from numpy.testing.utils import assert_allclose
 from ..core import Model, InputParameterError, custom_model, render_model
 from ..parameters import Parameter
 from .. import models
 
+from ...tests.helper import pytest, catch_warnings
 from ...utils.compat.funcsigs import signature
+from ...utils.exceptions import AstropyDeprecationWarning
 
 
 class NonFittableModel(Model):
@@ -200,12 +201,43 @@ def test_custom_inverse():
     assert_allclose(x, p(p.inverse(x)))
     assert_allclose(x, p.inverse(p(x)))
 
-    p.inverse = None
+    with catch_warnings(AstropyDeprecationWarning) as w:
+        p.inverse = None
+
+    # TODO: This can be removed after Astropy v1.1 or so
+    assert len(w) == 1
 
     with pytest.raises(NotImplementedError):
         p.inverse
 
 
+def test_custom_inverse_reset():
+    """Test resetting a custom inverse to the model's default inverse."""
+
+    class TestModel(Model):
+        inputs = ()
+        outputs = ('y',)
+        @property
+        def inverse(self):
+            return models.Shift()
+
+        @staticmethod
+        def evaluate():
+            return 0
+
+    # The above test model has no meaning, nor does its inverse--this just
+    # tests that setting an inverse and resetting to the default inverse works
+
+    m = TestModel()
+    assert isinstance(m.inverse, models.Shift)
+
+    m.inverse = models.Scale()
+    assert isinstance(m.inverse, models.Scale)
+
+    del m.inverse
+    assert isinstance(m.inverse, models.Shift)
+
+
 def test_render_model_2d():
 
     imshape = (71, 141)
 
@@ -30,6 +30,19 @@ a model may be instantiated with all scalar parameters::
     >>> g
     <Gaussian1D(amplitude=1.0, mean=0.0, stddev=1.0)>
 
+The newly created model instance ``g`` now works like a Gaussian function
+with the given parameters fixed.  It takes a single input::
+
+    >>> g.inputs
+    ('x',)
+    >>> g(x=0)
+    1.0
+
+The model can also be called without explicitly using keyword arguments::
+
+    >>> g(0)
+    1.0
+
 Or it may use all array parameters.  For example if all parameters are 2x2
 arrays the model is computed element-wise using all elements in the arrays::
 
@@ -181,24 +194,62 @@ set.  For scalar inputs like ``g(0)``, ``model_set_axis=False`` is implied
 automatically.  But for array inputs it is necessary to avoid ambiguity.
 
 
-Inputs and Outputs
-==================
-
-Models have an `~astropy.modeling.Model.n_inputs` attribute, which shows how
-many coordinates the model expects as an input. All models expect coordinates
-as separate arguments.  For example a 2-D model expects x and y coordinate
-values to be passed separately, i.e. as two scalars or array-like values.
-
-Models also have an attribute `~astropy.modeling.Model.n_outputs`, which shows
-the number of output coordinates. The `~astropy.modeling.Model.n_inputs` and
-`~astropy.modeling.Model.n_outputs` attributes can be used when chaining
-transforms by adding models in :class:`series
-<astropy.modeling.SerialCompositeModel>` or in :class:`parallel
-<astropy.modeling.SummedCompositeModel>`. Because composite models can be
-nested within other composite models, creating theoretically infinitely complex
-models, a mechanism to map input data to models is needed. In this case the
-input may be wrapped in a `~astropy.modeling.LabeledInput` object-- a dict-like
-object whose items are ``{label: data}`` pairs.
+Model Inverses
+==============
+
+All models have a `Model.inverse <astropy.modeling.Model.inverse>` property
+which may, for some models, return a new model that is the analytic inverse of
+the model it is attached to.  For example::
+
+    >>> from astropy.modeling.models import Linear1D
+    >>> linear = Linear1D(slope=0.8, intercept=1.0)
+    >>> linear.inverse
+    <Linear1D(slope=1.25, intercept=-1.25)>
+
+The inverse of a model will always be a fully instantiated model in its own
+right, and so can be evaluated directly like::
+
+    >>> linear.inverse(2.0)
+    1.25
+
+It is also possible to assign a *custom* inverse to a model.  This may be
+useful, for example, in cases where a model does not have an analytic inverse,
+but may have an approximate inverse that was computed numerically and is
+represented by a polynomial.  This works even if the target model has a
+default analytic inverse--in this case the default is overridden with the
+custom inverse::
+
+    >>> from astropy.modeling.models import Polynomial1D
+    >>> linear.inverse = Polynomial1D(degree=1, c0=-1.25, c1=1.25)
+    >>> linear.inverse
+    <Polynomial1D(1, c0=-1.25, c1=1.25)>
+
+If a custom inverse has been assigned to a model, it can be deleted with
+``del model.inverse``.  This resets the inverse to its default (if one exists).
+If a default does not exist, accessing ``model.inverse`` raises a
+`NotImplementedError`.  For example polynomial models do not have a default
+inverse::
+
+    >>> del linear.inverse
+    >>> linear.inverse
+    <Linear1D(slope=1.25, intercept=-1.25)>
+    >>> p = Polynomial1D(degree=2, c0=1.0, c1=2.0, c2=3.0)
+    >>> p.inverse
+    Traceback (most recent call last):
+      File "<stdin>", line 1, in <module>
+      File "astropy\modeling\core.py", line 796, in inverse
+        raise NotImplementedError("An analytical inverse transform has not "
+    NotImplementedError: An analytical inverse transform has not been
+    implemented for this model. 
+
+One may certainly compute an inverse and assign it to a polynomial model
+though.
+
+.. note::
+
+    When assigning a custom inverse to a model no validation is performed to
+    ensure that it is actually an inverse or even approximate inverse.  So
+    assign custom inverses at your own risk.
 
 
 Further examples
 
@@ -268,6 +268,11 @@ input``.
             new_intercept = -self.intercept / self.slope
             return LineModel(slope=new_slope, intercept=new_intercept)
 
+.. note::
+
+    The above example is essentially equivalent to the built-in
+    `~astropy.modeling.functional_models.Linear1D` model.
+
 
 Defining New Fitter Classes
 ===========================
 
@@ -63,3 +63,4 @@ py:obj list.remove
 py:class classmethod
 py:obj RuntimeError
 py:obj AttributeError
+py:obj NotImplementedError