Project the maximum likelihood values and sampled posterior points as 
    spectra.

    :param model:
        The model employed.

    :type model:
        :class:`sick.models.Model`

    :param data:
        The observed spectra.

    :type data:
        iterable of :class:`sick.specutils.Spectrum1D` objects

    :param theta: [optional]
        The optimised model parameters given the data. Either theta
        or chain should be given.

    :type theta:
        dict

    :param chain: [optional]
        The chain of sampled parameters.

    :type chain:
        :class:`numpy.ndarray`    

    :param extents: [optional]
        The wavelength extents to plot for each channel in the form of 
        [(min_chan_1, max_chan_1), ..., (min_chan_N, max_chan_N)]
    
    :type extents:
        tuple or None

    :param uncertainties: [optional]
        Show uncertainty of the data points.

    :type uncertainties:
        bool

    :param title: [optional]
        Title to set for the top axes.

    :type title:
        str

    :param fig: [optional]
        Overplot onto the provided figure object.

    :type fig:
        :class:`matplotlib.Figure` or None
    
    :param figsize: [optional]
        The figure size (x-dimension, y-dimension) in inches.

    :type figsize:
        tuple or None

    :raises ValueError:
        If a ``fig`` is provided with the incorrect number of axes.

    :raise TypeError:
        If the ``data`` are not provided in the correct type.

    :returns:
        The projection figure.

    :rtype:
        :class:`maplotlib.Figure`

def old_projection(model, data, theta=None, chain=None, n=100, extents=None, 
    uncertainties=True, title=None, fig=None, figsize=None):
    """
    Project the maximum likelihood values and sampled posterior points as 
    spectra.

    :param model:
        The model employed.

    :type model:
        :class:`sick.models.Model`

    :param data:
        The observed spectra.

    :type data:
        iterable of :class:`sick.specutils.Spectrum1D` objects

    :param theta: [optional]
        The optimised model parameters given the data. Either theta
        or chain should be given.

    :type theta:
        dict

    :param chain: [optional]
        The chain of sampled parameters.

    :type chain:
        :class:`numpy.ndarray`    

    :param extents: [optional]
        The wavelength extents to plot for each channel in the form of 
        [(min_chan_1, max_chan_1), ..., (min_chan_N, max_chan_N)]
    
    :type extents:
        tuple or None

    :param uncertainties: [optional]
        Show uncertainty of the data points.

    :type uncertainties:
        bool

    :param title: [optional]
        Title to set for the top axes.

    :type title:
        str

    :param fig: [optional]
        Overplot onto the provided figure object.

    :type fig:
        :class:`matplotlib.Figure` or None
    
    :param figsize: [optional]
        The figure size (x-dimension, y-dimension) in inches.

    :type figsize:
        tuple or None

    :raises ValueError:
        If a ``fig`` is provided with the incorrect number of axes.

    :raise TypeError:
        If the ``data`` are not provided in the correct type.

    :returns:
        The projection figure.

    :rtype:
        :class:`maplotlib.Figure`
    """
    if not isinstance(data, (tuple, list)) or \
    any([not isinstance(each, specutils.Spectrum1D) for each in data]):
        raise TypeError("Data must be a list-type of Spectrum1D objects.")

    K = len(data)

    factor = 2.0
    lbdim = 0.5 * factor
    trdim = 0.2 * factor
    whspace = 0.10
    width = np.max([len(each.disp) for each in data])/500.
    height = factor*K + factor * (K - 1.) * whspace
    dimy = lbdim + height + trdim
    dimx = lbdim + width + trdim

    if figsize is None:
        figsize = (dimx, dimy)
    if fig is None:
        fig, axes = plt.subplots(K, 1, figsize=figsize)

    else:
        try:
            axes = np.array(fig.axes).reshape((K, 1))
        except:
            raise ValueError("Provided figure has {0} axes, but data has "
                "parameters K={1}".format(len(fig.axes), K))

    if chain is not None:

        flat_chain = chain.reshape(-1, len(model.parameters))
        map_theta = np.mean(flat_chain, axis=0)

        try:
            map_fluxes = model(data=data, **model._dictify_theta(map_theta))
        except:
            logger.warn("Could not draw MAP fluxes from posterior")

        if n > 0:
            # Draw samples from sampler.chain and compute spectra for them
            sampled_fluxes = []
            n_samples = len(flat_chain)

            for i in range(n):
                sampled_theta = dict(zip(model.parameters,
                    flat_chain[np.random.randint(0, n_samples)]))
                try:
                    sampler_flux = model(data=data, **sampled_theta)
                except:
                    logger.warn("Could not draw sample flux from posterior")
                    continue
                else:
                    sampled_fluxes.append(sampler_flux)
        
    elif theta is not None:

        sampled_fluxes = []
        map_fluxes = model(data=data, **model._dictify_theta(theta))
        
    else:
        raise ValueError("either theta or chain should be given")


    if len(data) == 1:
        axes = [axes]

    for k, (map_flux, observed_spectrum) in enumerate(zip(map_fluxes, data)):

        ax = axes[k]

        # Draw the random samples from the chain
        if n > 0:
            for sampled_flux in sampled_fluxes:
                ax.plot(observed_spectrum.disp, sampled_flux[k], color="r",
                    alpha=0.5, zorder=90)

        # Draw the ML spectra
        ax.plot(observed_spectrum.disp, map_flux, color="r", lw=2, zorder=100)

        # Plot the data
        if uncertainties:
            ax.fill_between(observed_spectrum.disp,
                observed_spectrum.flux - observed_spectrum.variance**0.5,
                observed_spectrum.flux + observed_spectrum.variance**0.5,
                facecolor="#cccccc", edgecolor="#666666", zorder=-1)
        ax.plot(observed_spectrum.disp, observed_spectrum.flux, color="k",
            zorder=10)

        # By default only show common overlap between the model and data
        if extents is None:
            finite_data = np.isfinite(observed_spectrum.flux)
            finite_model = np.isfinite(map_flux)
            finite_points = (finite_model, finite_data)

            x_extent = [
                np.max([observed_spectrum.disp[s][0]  for s in finite_points]),
                np.min([observed_spectrum.disp[s][-1] for s in finite_points]),
            ]

            indices = observed_spectrum.disp.searchsorted(x_extent)
            finite_flux = observed_spectrum.flux[indices[0]:indices[1]]

            if len(finite_flux) > 0:
                #y_extent = [
                #    0.9 * np.min(finite_flux[np.isfinite(finite_flux)]),
                #    1.1 * np.max(finite_flux[np.isfinite(finite_flux)])
                #]
                ax.set_ylim([0.9, 1.1] * np.percentile(finite_flux[np.isfinite(finite_flux)], [0.5, 99.5]))

            ax.set_xlim(x_extent)

        else:
            ax.set_xlim(extents[k][0])
            ax.set_ylim(extents[k][1])

        # Labels and ticks
        if not (k < K - 1):
            ax.set_xlabel("Wavelength, $\lambda$ ($\AA$)")

        ax.set_ylabel("Flux, $F_\lambda$")
        ax.yaxis.set_label_coords(-0.05, 0.5)

        ax.xaxis.set_major_locator(MaxNLocator(5))
        ax.yaxis.set_major_locator(MaxNLocator(5))
        [l.set_rotation(45) for l in ax.get_yticklabels()]


    if title is not None and isinstance(title, (str, unicode)):
        axes[0].set_title(title)

    return fig