Repository: astropy/astropy

Issue:
Introduce `Model.separable`
Add a new boolean attribute, `separable` to simple models. `Model.separable` will indicate whether the inputs to a model are separable or not. This is used in GWCS to compute the separability of axes, e.g. spectral and spatial.

Discussion:
astropy-bot[bot]: Hi there @nden :wave: - thanks for the pull request! I'm just a friendly :robot: that checks for issues related to the changelog and making sure that this pull request is milestoned and labelled correctly. This is mainly intended for the maintainers, so if you are not a maintainer you can ignore this, and a maintainer will let you know if any action is required on your part :smiley:.

I see this is a work in progress pull request. I'll report back on the checks once the PR is ready for review.

*If there are any issues with this message, please report them [here](https://github.com/astropy/astropy-bot/issues)*
bsipocz: (@nden - I've removed the "affects dev" label, as we supposed to use it for bugfixes that are only in the dev version, thus they don't need a changelog entry. 
This one on the other hand seems to be an added feature, so the changelog entry is definitely required.
nden: @bsipocz Oh, right. Thanks.
nden: @bsipocz , @pllim I've addressed your comments, hopefully things are clearer now. Thanks for looking over the PR.
pllim: Need to rebase to resolve conflict.

@astrofrog , do you think this is okay to merge after rebase?

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index e54f01820e..48d0604d2d 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -85,6 +85,12 @@ astropy.modeling
 - A ``deepcopy()`` method was added to models. [#6515]
 
 - Added units support to ``AffineTransformation``. [#6853]
+- Added ``is_separable`` function to modeling totest the
+- Added ``is_separable`` function to modeling to test the
+  separability of a model. [#6746]
+
+- Added ``Model.separable`` property. It returns a boolean value or
+  ``None`` if not set. [#6746]
 
 astropy.nddata
 ^^^^^^^^^^^^^^
diff --git a/astropy/modeling/__init__.py b/astropy/modeling/__init__.py
index 6d21ff3f90..75141f5bbe 100644
--- a/astropy/modeling/__init__.py
+++ b/astropy/modeling/__init__.py
@@ -11,3 +11,4 @@ from . import fitting
 from . import models
 from .core import *
 from .parameters import *
+from .separable import *
diff --git a/astropy/modeling/core.py b/astropy/modeling/core.py
index 5aee030750..0d2745a1fb 100644
--- a/astropy/modeling/core.py
+++ b/astropy/modeling/core.py
@@ -633,6 +633,9 @@ class Model(metaclass=_ModelMeta):
     fittable = False
     linear = True
 
+    _separable = None
+    """ A boolean flag to indicate whether a model is separable."""
+
     meta = metadata.MetaData()
     """A dict-like object to store optional information."""
 
@@ -1111,6 +1114,17 @@ class Model(metaclass=_ModelMeta):
 
         return self._user_bounding_box is not None
 
+    @property
+    def separable(self):
+        """ A flag indicating whether a model is separable."""
+
+        if self._separable is not None:
+            return self._separable
+        else:
+            raise NotImplementedError(
+                'The "separable" property is not defined for '
+                'model {}'.format(self.__class__.__name__))
+
     # *** Public methods ***
 
     def without_units_for_data(self, **kwargs):
@@ -2080,6 +2094,7 @@ class Fittable1DModel(FittableModel):
 
     inputs = ('x',)
     outputs = ('y',)
+    _separable = True
 
 
 class Fittable2DModel(FittableModel):
diff --git a/astropy/modeling/polynomial.py b/astropy/modeling/polynomial.py
index 933fed2ad6..07d81fda64 100644
--- a/astropy/modeling/polynomial.py
+++ b/astropy/modeling/polynomial.py
@@ -4,7 +4,6 @@
 This module contains models representing polynomials and polynomial series.
 """
 
-
 from collections import OrderedDict
 
 import numpy as np
@@ -352,6 +351,7 @@ class Chebyshev1D(PolynomialModel):
 
     inputs = ('x',)
     outputs = ('y',)
+    _separable = True
 
     def __init__(self, degree, domain=None, window=[-1, 1], n_models=None,
                  model_set_axis=None, name=None, meta=None, **params):
@@ -458,6 +458,7 @@ class Hermite1D(PolynomialModel):
 
     inputs = ('x')
     outputs = ('y')
+    _separable = True
 
     def __init__(self, degree, domain=None, window=[-1, 1], n_models=None,
                  model_set_axis=None, name=None, meta=None, **params):
@@ -566,6 +567,7 @@ class Hermite2D(OrthoPolynomialBase):
     example, the third Hermite polynomial (H2) is 4x^2-2, but if x was
     specified with units, 4x^2 and -2 would have incompatible units.
     """
+    _separable = False
 
     def __init__(self, x_degree, y_degree, x_domain=None, x_window=[-1, 1],
                  y_domain=None, y_window=[-1, 1], n_models=None,
@@ -688,6 +690,7 @@ class Legendre1D(PolynomialModel):
 
     inputs = ('x',)
     outputs = ('y',)
+    _separable = False
 
     def __init__(self, degree, domain=None, window=[-1, 1], n_models=None,
                  model_set_axis=None, name=None, meta=None, **params):
@@ -781,6 +784,7 @@ class Polynomial1D(PolynomialModel):
 
     inputs = ('x',)
     outputs = ('y',)
+    _separable = True
 
     def __init__(self, degree, domain=[-1, 1], window=[-1, 1], n_models=None,
                  model_set_axis=None, name=None, meta=None, **params):
@@ -881,6 +885,7 @@ class Polynomial2D(PolynomialModel):
 
     inputs = ('x', 'y')
     outputs = ('z',)
+    _separable = False
 
     def __init__(self, degree, x_domain=[-1, 1], y_domain=[-1, 1],
                  x_window=[-1, 1], y_window=[-1, 1], n_models=None,
@@ -1057,6 +1062,7 @@ class Chebyshev2D(OrthoPolynomialBase):
     example, the third Chebyshev polynomial (T2) is 2x^2-1, but if x was
     specified with units, 2x^2 and -1 would have incompatible units.
     """
+    _separable = False
 
     def __init__(self, x_degree, y_degree, x_domain=None, x_window=[-1, 1],
                  y_domain=None, y_window=[-1, 1], n_models=None,
@@ -1187,6 +1193,7 @@ class Legendre2D(OrthoPolynomialBase):
     third Legendre polynomial (P2) is 1.5x^2-0.5, but if x was specified with
     units, 1.5x^2 and -0.5 would have incompatible units.
     """
+    _separable = False
 
     def __init__(self, x_degree, y_degree, x_domain=None, x_window=[-1, 1],
                  y_domain=None, y_window=[-1, 1], n_models=None,
@@ -1276,6 +1283,8 @@ class _SIP1D(PolynomialBase):
 
     inputs = ('u', 'v')
     outputs = ('w',)
+    _separable = False
+
 
     def __init__(self, order, coeff_prefix, n_models=None,
                  model_set_axis=None, name=None, meta=None, **params):
@@ -1390,6 +1399,7 @@ class SIP(Model):
 
     inputs = ('u', 'v')
     outputs = ('x', 'y')
+    _separable = False
 
     def __init__(self, crpix, a_order, b_order, a_coeff={}, b_coeff={},
                  ap_order=None, bp_order=None, ap_coeff={}, bp_coeff={},
@@ -1459,6 +1469,7 @@ class InverseSIP(Model):
 
     inputs = ('x', 'y')
     outputs = ('u', 'v')
+    _separable = False
 
     def __init__(self, ap_order, bp_order, ap_coeff={}, bp_coeff={},
                  n_models=None, model_set_axis=None, name=None, meta=None):
diff --git a/astropy/modeling/projections.py b/astropy/modeling/projections.py
index cb8b180771..63ddf5d2af 100644
--- a/astropy/modeling/projections.py
+++ b/astropy/modeling/projections.py
@@ -106,6 +106,8 @@ class Projection(Model):
     # This sets the circumference to 360 deg so that arc length is measured in deg.
     r0 = 180 * u.deg / np.pi
 
+    _separable = False
+
     @property
     @abc.abstractmethod
     def inverse(self):
@@ -170,6 +172,8 @@ class Zenithal(Projection):
         y &= R_\theta \cos \phi
     """
 
+    _separable = False
+
 
 class Pix2Sky_ZenithalPerspective(Pix2SkyProjection, Zenithal):
     r"""
@@ -702,6 +706,7 @@ class Cylindrical(Projection):
     Cylindrical projections are so-named because the surface of
     projection is a cylinder.
     """
+    _separable = True
 
 
 class Pix2Sky_CylindricalPerspective(Pix2SkyProjection, Cylindrical):
@@ -972,6 +977,8 @@ class PseudoCylindrical(Projection):
     distortion there.  Consequently, the meridians are curved.
     """
 
+    _separable = True
+
 
 class Pix2Sky_SansonFlamsteed(Pix2SkyProjection, PseudoCylindrical):
     r"""
@@ -1030,6 +1037,8 @@ class Pix2Sky_Parabolic(Pix2SkyProjection, PseudoCylindrical):
         \theta &= 3 \sin^{-1}\left(\frac{y}{180^\circ}\right)
     """
 
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_Parabolic()
@@ -1053,6 +1062,8 @@ class Sky2Pix_Parabolic(Sky2PixProjection, PseudoCylindrical):
         y &= 180^\circ \sin \frac{\theta}{3}
     """
 
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_Parabolic()
@@ -1076,6 +1087,8 @@ class Pix2Sky_Molleweide(Pix2SkyProjection, PseudoCylindrical):
         \theta &= \sin^{-1}\left(\frac{1}{90^\circ}\sin^{-1}\left(\frac{\pi}{180^\circ}\frac{y}{\sqrt{2}}\right) + \frac{y}{180^\circ}\sqrt{2 - \left(\frac{\pi}{180^\circ}y\right)^2}\right)
     """
 
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_Molleweide()
@@ -1106,6 +1119,8 @@ class Sky2Pix_Molleweide(Sky2PixProjection, PseudoCylindrical):
         \sin \theta = \frac{\gamma}{90^\circ} + \frac{\sin 2 \gamma}{\pi}
     """
 
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_Molleweide()
@@ -1129,6 +1144,8 @@ class Pix2Sky_HammerAitoff(Pix2SkyProjection, PseudoCylindrical):
         \theta &= \sin^{-1}\left(\frac{\pi}{180^\circ}yZ\right)
     """
 
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_HammerAitoff()
@@ -1157,6 +1174,8 @@ class Sky2Pix_HammerAitoff(Sky2PixProjection, PseudoCylindrical):
         \gamma = \frac{180^\circ}{\pi} \sqrt{\frac{2}{1 + \cos \theta \cos(\phi / 2)}}
     """
 
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_HammerAitoff()
@@ -1196,6 +1215,8 @@ class Conic(Projection):
     sigma = Parameter(default=90.0, getter=_to_orig_unit, setter=_to_radian)
     delta = Parameter(default=0.0, getter=_to_orig_unit, setter=_to_radian)
 
+    _separable = False
+
 
 class Pix2Sky_ConicPerspective(Pix2SkyProjection, Conic):
     r"""
@@ -1576,6 +1597,7 @@ class Pix2Sky_BonneEqualArea(Pix2SkyProjection, PseudoConic):
         Bonne conformal latitude, in degrees.
     """
     theta1 = Parameter(default=0.0, getter=_to_orig_unit, setter=_to_radian)
+    _separable = True
 
     @property
     def inverse(self):
@@ -1612,6 +1634,7 @@ class Sky2Pix_BonneEqualArea(Sky2PixProjection, PseudoConic):
         Bonne conformal latitude, in degrees.
     """
     theta1 = Parameter(default=0.0, getter=_to_orig_unit, setter=_to_radian)
+    _separable = True
 
     @property
     def inverse(self):
@@ -1632,6 +1655,8 @@ class Pix2Sky_Polyconic(Pix2SkyProjection, PseudoConic):
 
     Corresponds to the ``PCO`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_Polyconic()
@@ -1650,6 +1675,8 @@ class Sky2Pix_Polyconic(Sky2PixProjection, PseudoConic):
 
     Corresponds to the ``PCO`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_Polyconic()
@@ -1685,6 +1712,8 @@ class Pix2Sky_TangentialSphericalCube(Pix2SkyProjection, QuadCube):
 
     Corresponds to the ``TSC`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_TangentialSphericalCube()
@@ -1703,6 +1732,8 @@ class Sky2Pix_TangentialSphericalCube(Sky2PixProjection, QuadCube):
 
     Corresponds to the ``PCO`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_TangentialSphericalCube()
@@ -1721,6 +1752,8 @@ class Pix2Sky_COBEQuadSphericalCube(Pix2SkyProjection, QuadCube):
 
     Corresponds to the ``CSC`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_COBEQuadSphericalCube()
@@ -1739,6 +1772,8 @@ class Sky2Pix_COBEQuadSphericalCube(Sky2PixProjection, QuadCube):
 
     Corresponds to the ``CSC`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_COBEQuadSphericalCube()
@@ -1757,6 +1792,8 @@ class Pix2Sky_QuadSphericalCube(Pix2SkyProjection, QuadCube):
 
     Corresponds to the ``QSC`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_QuadSphericalCube()
@@ -1775,6 +1812,8 @@ class Sky2Pix_QuadSphericalCube(Sky2PixProjection, QuadCube):
 
     Corresponds to the ``QSC`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_QuadSphericalCube()
@@ -1806,6 +1845,8 @@ class Pix2Sky_HEALPix(Pix2SkyProjection, HEALPix):
     X : float
         The number of facets in latitude direction.
     """
+    _separable = True
+
     H = Parameter(default=4.0)
     X = Parameter(default=3.0)
 
@@ -1835,6 +1876,8 @@ class Sky2Pix_HEALPix(Sky2PixProjection, HEALPix):
     X : float
         The number of facets in latitude direction.
     """
+    _separable = True
+
     H = Parameter(default=4.0)
     X = Parameter(default=3.0)
 
@@ -1856,6 +1899,8 @@ class Pix2Sky_HEALPixPolar(Pix2SkyProjection, HEALPix):
 
     Corresponds to the ``XPH`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Sky2Pix_HEALPix()
@@ -1874,6 +1919,8 @@ class Sky2Pix_HEALPixPolar(Sky2PixProjection, HEALPix):
 
     Corresponds to the ``XPH`` projection in FITS WCS.
     """
+    _separable = False
+
     @property
     def inverse(self):
         return Pix2Sky_HEALPix()
@@ -1906,6 +1953,8 @@ class AffineTransformation2D(Model):
 
     standard_broadcasting = False
 
+    _separable = False
+
     matrix = Parameter(default=[[1.0, 0.0], [0.0, 1.0]])
     translation = Parameter(default=[0.0, 0.0])
 
diff --git a/astropy/modeling/rotations.py b/astropy/modeling/rotations.py
index eb5be04ee2..81678d5679 100644
--- a/astropy/modeling/rotations.py
+++ b/astropy/modeling/rotations.py
@@ -19,7 +19,6 @@ References
 .. [1] Calabretta, M.R., Greisen, E.W., 2002, A&A, 395, 1077 (Paper II)
 """
 
-
 import math
 
 import numpy as np
@@ -40,6 +39,8 @@ class _EulerRotation:
     Base class which does the actual computation.
     """
 
+    _separable = False
+
     def _create_matrix(self, phi, theta, psi, axes_order):
         matrices = []
         for angle, axis in zip([phi, theta, psi], axes_order):
@@ -342,6 +343,7 @@ class Rotation2D(Model):
 
     inputs = ('x', 'y')
     outputs = ('x', 'y')
+    _separable = False
 
     angle = Parameter(default=0.0, getter=_to_orig_unit, setter=_to_radian)
 
diff --git a/astropy/modeling/separable.py b/astropy/modeling/separable.py
new file mode 100644
index 0000000000..687d04241e
--- /dev/null
+++ b/astropy/modeling/separable.py
@@ -0,0 +1,277 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+Functions to determine if a model is separable, i.e.
+if the model outputs are independent.
+
+It analyzes ``n_inputs``, ``n_outputs`` and the operators
+in a compound model by stepping through the transforms
+and creating a ``coord_matrix`` of shape (``n_outputs``, ``n_inputs``).
+
+
+Each modeling operator is represented by a function which
+takes two simple models (or two ``coord_matrix`` arrays) and
+returns an array of shape (``n_outputs``, ``n_inputs``).
+
+"""
+
+import numpy as np
+
+from .core import Model, _CompoundModel, ModelDefinitionError
+from .mappings import Mapping
+
+
+__all__ = ["is_separable"]
+
+
+def is_separable(transform):
+    """
+    A separability test for the outputs of a transform.
+
+    Parameters
+    ----------
+    transform : `~astropy.modeling.core.Model`
+        A (compound) model.
+
+    Returns
+    -------
+    is_separable : ndarray
+        A boolean array with size ``transform.n_outputs`` where
+        each element indicates whether the output is independent
+        and the result of a separable transform.
+
+    Examples
+    --------
+    >>> from astropy.modeling.models import Shift, Scale, Rotation2D, Polynomial2D
+    >>> is_separable(Shift(1) & Shift(2) | Scale(1) & Scale(2))
+        array([ True,  True]...)
+    >>> is_separable(Shift(1) & Shift(2) | Rotation2D(2))
+        array([False, False]...)
+    >>> is_separable(Shift(1) & Shift(2) | Mapping([0, 1, 0, 1]) | \
+        Polynomial2D(1) & Polynomial2D(2))
+        array([False, False]...)
+    >>> is_separable(Shift(1) & Shift(2) | Mapping([0, 1, 0, 1]))
+        array([ True,  True,  True,  True]...)
+
+    """
+    if transform.n_inputs == 1 and transform.n_outputs > 1:
+        is_separable = np.array([False] * transform.n_outputs).T
+        return is_separable
+    separable_matrix = _separable(transform)
+    is_separable = separable_matrix.sum(1)
+    is_separable = np.where(is_separable != 1, False, True)
+    return is_separable
+
+
+def _compute_n_outputs(left, right):
+    """
+    Compute the number of outputs of two models.
+
+    The two models are the left and right model to an operation in
+    the expression tree of a compound model.
+
+    Parameters
+    ----------
+    left, right : `astropy.modeling.Model` or ndarray
+        If input is of an array, it is the output of `coord_matrix`.
+
+    """
+    if isinstance(left, Model):
+        lnout = left.n_outputs
+    else:
+        lnout = left.shape[0]
+    if isinstance(right, Model):
+        rnout = right.n_outputs
+    else:
+        rnout = right.shape[0]
+    noutp = lnout + rnout
+    return noutp
+
+
+def _arith_oper(left, right):
+    """
+    Function corresponding to one of the arithmetic operators
+    ['+', '-'. '*', '/', '**'].
+
+    This always returns a nonseparable output.
+
+
+    Parameters
+    ----------
+    left, right : `astropy.modeling.Model` or ndarray
+        If input is of an array, it is the output of `coord_matrix`.
+
+    Returns
+    -------
+    result : ndarray
+        Result from this operation.
+    """
+    # models have the same number of inputs and outputs
+    def _n_inputs_outputs(input):
+        if isinstance(input, Model):
+            n_outputs, n_inputs = input.n_outputs, input.n_inputs
+        else:
+            n_outputs, n_inputs = input.shape
+        return n_inputs, n_outputs
+
+
+    left_inputs, left_outputs = _n_inputs_outputs(left)
+    right_inputs, right_outputs = _n_inputs_outputs(right)
+
+    if left_inputs != right_inputs or left_outputs != right_outputs:
+        raise ModelDefinitionError(
+            "Unsupported operands for arithmetic operator: left (n_inputs={0}, "
+            "n_outputs={1}) and right (n_inputs={2}, n_outputs={3}); "
+            "models must have the same n_inputs and the same "
+            "n_outputs for this operator.".format(
+                left_inputs, left_outputs, right_inputs, right_outputs))
+
+    result = np.ones((left_outputs, left_inputs))
+    return result
+
+
+def _coord_matrix(model, pos, noutp):
+    """
+    Create an array representing inputs and outputs of a simple model.
+
+    The array has a shape (noutp, model.n_inputs).
+
+    Parameters
+    ----------
+    model : `astropy.modeling.Model`
+        model
+    pos : str
+        Position of this model in the expression tree.
+        One of ['left', 'right'].
+    noutp : int
+        Number of outputs of the compound model of which the input model
+        is a left or right child.
+
+    """
+    if isinstance(model, Mapping):
+        axes = []
+        for i in model.mapping:
+            axis = np.zeros((model.n_inputs,))
+            axis[i] = 1
+            axes.append(axis)
+        m = np.vstack(axes)
+        mat = np.zeros((noutp, model.n_inputs))
+        if pos == 'left':
+            mat[: model.n_outputs, :model.n_inputs] = m
+        else:
+            mat[-model.n_outputs:, -model.n_inputs:] = m
+        return mat
+    if not model.separable:
+        # this does not work for more than 2 coordinates
+        mat = np.zeros((noutp, model.n_inputs))
+        if pos == 'left':
+            mat[:model.n_outputs, : model.n_inputs] = 1
+        else:
+            mat[-model.n_outputs:, -model.n_inputs:] = 1
+    else:
+        mat = np.zeros((noutp, model.n_inputs))
+
+        for i in range(model.n_inputs):
+            mat[i, i] = 1
+        if pos == 'right':
+            mat = np.roll(mat, (noutp - model.n_outputs))
+    return mat
+
+
+def _cstack(left, right):
+    """
+    Function corresponding to '&' operation.
+
+    Parameters
+    ----------
+    left, right : `astropy.modeling.Model` or ndarray
+        If input is of an array, it is the output of `coord_matrix`.
+
+    Returns
+    -------
+    result : ndarray
+        Result from this operation.
+
+    """
+    noutp = _compute_n_outputs(left, right)
+
+    if isinstance(left, Model):
+        cleft = _coord_matrix(left, 'left', noutp)
+    else:
+        cleft = np.zeros((noutp, left.shape[1]))
+        cleft[: left.shape[0], : left.shape[1]] = left
+    if isinstance(right, Model):
+        cright = _coord_matrix(right, 'right', noutp)
+    else:
+        cright = np.zeros((noutp, right.shape[1]))
+        cright[-right.shape[0]:, -right.shape[1]:] = 1
+
+    return np.hstack([cleft, cright])
+
+
+def _cdot(left, right):
+    """
+    Function corresponding to "|" operation.
+
+    Parameters
+    ----------
+    left, right : `astropy.modeling.Model` or ndarray
+        If input is of an array, it is the output of `coord_matrix`.
+
+    Returns
+    -------
+    result : ndarray
+        Result from this operation.
+    """
+
+    left, right = right, left
+
+    def _n_inputs_outputs(input, position):
+        """
+        Return ``n_inputs``, ``n_outputs`` for a model or coord_matrix.
+        """
+        if isinstance(input, Model):
+            coords = _coord_matrix(input, position, input.n_outputs)
+        else:
+            coords = input
+        return coords
+
+    cleft = _n_inputs_outputs(left, 'left')
+    cright = _n_inputs_outputs(right, 'right')
+
+    try:
+        result = np.dot(cleft, cright)
+    except ValueError:
+        raise ModelDefinitionError(
+            'Models cannot be combined with the "|" operator; '
+            'left coord_matrix is {0}, right coord_matrix is {1}'.format(
+                cright, cleft))
+    return result
+
+
+def _separable(transform):
+    """
+    Calculate the separability of outputs.
+
+    Parameters
+    ----------
+    transform : `astropy.modeling.Model`
+        A transform (usually a compound model).
+
+    Returns
+    -------
+    is_separable : ndarray of dtype np.bool
+        An array of shape (transform.n_outputs,) of boolean type
+        Each element represents the separablity of the corresponding output.
+    """
+    if isinstance(transform, _CompoundModel):
+        is_separable = transform._tree.evaluate(_operators)
+    elif isinstance(transform, Model):
+        is_separable = _coord_matrix(transform, 'left', transform.n_outputs)
+    return is_separable
+
+
+# Maps modeling operators to a function computing and represents the
+# relationship of axes as an array of 0-es and 1-s
+_operators = {'&': _cstack, '|': _cdot, '+': _arith_oper, '-': _arith_oper,
+              '*': _arith_oper, '/': _arith_oper, '**': _arith_oper}
diff --git a/astropy/modeling/tabular.py b/astropy/modeling/tabular.py
index 5dba0323ff..152d88cccb 100644
--- a/astropy/modeling/tabular.py
+++ b/astropy/modeling/tabular.py
@@ -277,6 +277,12 @@ def tabular_model(dim, name=None):
     table = np.zeros([2] * dim)
     inputs = tuple('x{0}'.format(idx) for idx in range(table.ndim))
     members = {'lookup_table': table, 'inputs': inputs}
+
+    if dim == 1:
+        members['_separable'] = True
+    else:
+        members['_separable'] = False
+
     if name is None:
         model_id = _Tabular._id
         _Tabular._id += 1
diff --git a/docs/modeling/index.rst b/docs/modeling/index.rst
index 9aaa62609c..52cc84e2af 100644
--- a/docs/modeling/index.rst
+++ b/docs/modeling/index.rst
@@ -381,3 +381,4 @@ Reference/API
 .. automodapi:: astropy.modeling.fitting
 .. automodapi:: astropy.modeling.optimizers
 .. automodapi:: astropy.modeling.statistic
+.. automodapi:: astropy.modeling.separable
diff --git a/docs/modeling/models.rst b/docs/modeling/models.rst
index 2d9b5e4dd1..15c34b1da1 100644
--- a/docs/modeling/models.rst
+++ b/docs/modeling/models.rst
@@ -251,6 +251,21 @@ though.
     ensure that it is actually an inverse or even approximate inverse.  So
     assign custom inverses at your own risk.
 
+.. _separability:
+
+Model Separability
+==================
+
+Simple models have a boolean `Model.separable <astropy.modeling.Model.separable>` property.
+It indicates whether the outputs are independent and is essential for computing the
+separability of compound models using the :func:`~astropy.modeling.is_separable` function.
+Having a separable compound model means that it can be decomposed into independent models,
+which in turn is useful in many applications.
+For example, it may be easier to define inverses using the independent parts of a model
+than the entire model.
+In other cases, tools using `GWCS <https://gwcs.readthedocs.io/en/latest/>`_,
+can be more flexible and take advantage of separable spectral and spatial transforms.
+
 
 Further examples
 ================
diff --git a/docs/modeling/new.rst b/docs/modeling/new.rst
index 8e99657ae3..cabe84f464 100644
--- a/docs/modeling/new.rst
+++ b/docs/modeling/new.rst
@@ -124,6 +124,11 @@ value of the `~astropy.modeling.Model.linear` attribute is ``False``.  Linear
 models should define the ``linear`` class attribute as ``True``.  Because this
 model is non-linear we can stick with the default.
 
+Models which inherit from `~astropy.modeling.Fittable1DModel` have the
+``Model._separable`` property already set to ``True``.
+All other models should define this property to indicate the
+:ref:`separability`.
+
 Next, provide methods called ``evaluate`` to evaluate the model and
 ``fit_deriv``, to compute its derivatives with respect to parameters.  These
 may be normal methods, `classmethod`, or `staticmethod`, though the convention


Test patch:
diff --git a/astropy/modeling/tests/test_separable.py b/astropy/modeling/tests/test_separable.py
new file mode 100644
index 0000000000..a8a9c93413
--- /dev/null
+++ b/astropy/modeling/tests/test_separable.py
@@ -0,0 +1,113 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+Test separability of models.
+
+"""
+import pytest
+import numpy as np
+from numpy.testing import assert_allclose
+
+from .. import models
+from ..models import Mapping
+from .. separable import (_coord_matrix, is_separable, _cdot,
+                          _cstack, _arith_oper)
+
+
+sh1 = models.Shift(1, name='shift1')
+sh2 = models.Shift(2, name='sh2')
+scl1 = models.Scale(1, name='scl1')
+scl2 = models.Scale(2, name='scl2')
+map1 = Mapping((0, 1, 0, 1), name='map1')
+map2 = Mapping((0, 0, 1), name='map2')
+map3 = Mapping((0, 0), name='map3')
+rot = models.Rotation2D(2, name='rotation')
+p2 = models.Polynomial2D(1, name='p2')
+p22 = models.Polynomial2D(2, name='p22')
+p1 = models.Polynomial1D(1, name='p1')
+
+
+compound_models = {
+    'cm1': (map3 & sh1 | rot & sh1 | sh1 & sh2 & sh1,
+            np.array([False, False, True])
+            ),
+    'cm2': (sh1 & sh2 | rot | map1 | p2 & p22,
+            np.array([False, False])
+            ),
+    'cm3': (map2 | rot & scl1,
+            np.array([False, False, True])
+            ),
+    'cm4': (sh1 & sh2 | map2 | rot & scl1,
+            np.array([False, False, True])
+            ),
+    'cm5': (map3 | sh1 & sh2 | scl1 & scl2,
+            np.array([False, False])
+            ),
+    'cm7': (map2 | p2 & sh1,
+            np. array([False, True])
+            )
+}
+
+
+def test_coord_matrix():
+    c = _coord_matrix(p2, 'left', 2)
+    assert_allclose(np.array([[1, 1], [0, 0]]), c)
+    c = _coord_matrix(p2, 'right', 2)
+    assert_allclose(np.array([[0, 0], [1, 1]]), c)
+    c = _coord_matrix(p1, 'left', 2)
+    assert_allclose(np.array([[1], [0]]), c)
+    c = _coord_matrix(p1, 'left', 1)
+    assert_allclose(np.array([[1]]), c)
+    c = _coord_matrix(sh1, 'left', 2)
+    assert_allclose(np.array([[1], [0]]), c)
+    c = _coord_matrix(sh1, 'right', 2)
+    assert_allclose(np.array([[0], [1]]), c)
+    c = _coord_matrix(sh1, 'right', 3)
+    assert_allclose(np.array([[0], [0], [1]]), c)
+    c = _coord_matrix(map3, 'left', 2)
+    assert_allclose(np.array([[1], [1]]), c)
+    c = _coord_matrix(map3, 'left', 3)
+    assert_allclose(np.array([[1], [1], [0]]), c)
+
+
+def test_cdot():
+    result = _cdot(sh1, scl1)
+    assert_allclose(result, np.array([[1]]))
+
+    result = _cdot(rot, p2)
+    assert_allclose(result, np.array([[2, 2]]))
+
+    result = _cdot(rot, rot)
+    assert_allclose(result, np.array([[2, 2], [2, 2]]))
+
+    result = _cdot(Mapping((0, 0)), rot)
+    assert_allclose(result, np.array([[2], [2]]))
+
+
+def test_cstack():
+    result = _cstack(sh1, scl1)
+    assert_allclose(result, np.array([[1, 0], [0, 1]]))
+
+    result = _cstack(sh1, rot)
+    assert_allclose(result,
+                    np.array([[1, 0, 0],
+                              [0, 1, 1],
+                              [0, 1, 1]])
+                    )
+    result = _cstack(rot, sh1)
+    assert_allclose(result,
+                    np.array([[1, 1, 0],
+                              [1, 1, 0],
+                              [0, 0, 1]])
+                    )
+
+
+def test_arith_oper():
+    result = _arith_oper(sh1, scl1)
+    assert_allclose(result, np.array([[1]]))
+    result = _arith_oper(rot, rot)
+    assert_allclose(result, np.array([[1, 1], [1, 1]]))
+
+
+@pytest.mark.parametrize(('compound_model', 'result'), compound_models.values())
+def test_separable(compound_model, result):
+    assert_allclose(is_separable(compound_model), result)