Make the right decisions for prediction with a model
    based on the standard arguments of plotting.

    This is quite complex and will take a while to understand,
    so do not change anything in here lightly!!!

def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, which_data_ycols, predict_kw, samples=0):
    """
    Make the right decisions for prediction with a model
    based on the standard arguments of plotting.

    This is quite complex and will take a while to understand,
    so do not change anything in here lightly!!!
    """
    # Put some standards into the predict_kw so that prediction is done automatically:
    if predict_kw is None:
        predict_kw = {}
    if 'likelihood' not in predict_kw:
        if plot_raw:
            from ...likelihoods import Gaussian
            from ...likelihoods.link_functions import Identity
            lik = Gaussian(Identity(), 1e-9) # Make the likelihood not add any noise
        else:
            lik = None
        predict_kw['likelihood'] = lik
    if 'Y_metadata' not in predict_kw:
        predict_kw['Y_metadata'] = {}
    if 'output_index' not in predict_kw['Y_metadata']:
        predict_kw['Y_metadata']['output_index'] = Xgrid[:,-1:].astype(np.int)

    mu, _ = self.predict(Xgrid, **predict_kw)

    if percentiles is not None:
        percentiles = self.predict_quantiles(Xgrid, quantiles=percentiles, **predict_kw)
    else: percentiles = []

    if samples > 0:
        fsamples = self.posterior_samples(Xgrid, size=samples, **predict_kw)
        fsamples = fsamples[:, which_data_ycols, :]
    else:
        fsamples = None

    # Filter out the ycolums which we want to plot:
    retmu = mu[:, which_data_ycols]
    percs = [p[:, which_data_ycols] for p in percentiles]

    if plot_raw and apply_link:
        for i in range(len(which_data_ycols)):
            retmu[:, [i]] = self.likelihood.gp_link.transf(mu[:, [i]])
            for perc in percs:
                perc[:, [i]] = self.likelihood.gp_link.transf(perc[:, [i]])
            if fsamples is not None:
                for s in range(fsamples.shape[-1]):
                    fsamples[:, i, s] = self.likelihood.gp_link.transf(fsamples[:, i, s])
    return retmu, percs, fsamples