Optimise the model parameters, given the data.

def optimise(self, data, initial_theta=None, full_output=False, **kwargs):
        """
        Optimise the model parameters, given the data.
        """

        data = self._format_data(data)

        if initial_theta is None:
            initial_theta = self.estimate(data)

        # Which parameters will be optimised, and which will be fixed?
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)

        # For the Cannon model, the normalised fluxes are a multiplicative
        # operation, so we don't actually need to optimise the astrophysical
        # parameters, we only need to optimise the other ones. At each set of
        # non-astrophysical parameters we can determine the optimal
        # astrophysical parameters by a matrix inversion.
        parameters = [p for p in self.parameters if p not in ignore_parameters \
            and p not in self.grid_points.dtype.names]

        # What model wavelength ranges will be required?
        wavelengths_required = []
        for channel, spectrum in zip(matched_channels, data):
            if channel is None: continue
            z = initial_theta.get("z", initial_theta.get("z_{}".format(channel),
                0))
            wavelengths_required.append(
                [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])

        # (If necessary,..) train the Cannon model around the closest point.
        subset_bounds = self._initialise_approximator(closest_point=\
            [initial_theta[p] for p in self.grid_points.dtype.names], 
            wavelengths_required=wavelengths_required, **kwargs)
        
        # Get the optimisation keyword arguments.
        op_kwargs = self._configuration.get("optimise", {}).copy()
        op_kwargs.update(kwargs)

        # Get fixed keywords.
        fixed = op_kwargs.pop("fixed", {})
        if fixed is None: fixed = {}
        if fixed:
            # Remove non-parameters from the 'fixed' keywords.
            keys = set(fixed).intersection(self.parameters)
            ignored_keys = set(fixed).difference(parameters)
            if ignored_keys:
                logger.warn("Ignoring optimise.fixed keys because they are not"
                    " model parameters: {0}".format(", ".join(ignored_keys)))
            # If the 'fixed' value is provided, use that. Otherwise if it is
            # None then use the initial_theta value.
            fixed = dict(zip(keys, 
                [(fixed[k], initial_theta.get(k, None))[fixed[k] is None] \
                    for k in keys]))

            logger.info("Fixing keyword arguments (these will not be optimised)"
                ": {}".format(fixed))

        # Remove fixed parameters from the parameters to be optimised
        #parameters = list(set(parameters).difference(fixed))
        parameters = [p for p in parameters if p not in fixed]

        # Translate input bounds.
        nbs = (None, None) # No boundaries.
        input_bounds = op_kwargs.pop("bounds", {})
        # Because astrophysical parameters are not being 'optimised', we may not
        # actually need the subset_bounds.get call, but there may be fringe
        # cases that I cannot immediately think of.
        op_kwargs["bounds"] = [input_bounds.get(p, subset_bounds.get(p, nbs)) \
            for p in parameters]

        # Apply data masks now so we don't have to do it on the fly.
        masked_data, pixels_affected = self._apply_data_mask(data)

        # Prepare the convolution functions.
        self._create_convolution_functions(matched_channels, data, parameters)

        logger.info("Optimising parameters: {0}".format(", ".join(parameters)))
        logger.info("Optimisation keywords: {0}".format(op_kwargs))

        # Create the objective function. The objective function needs to
        # continuum-normalise the observed spectra, put it at rest, then solve
        # for the best astrophysical parameters, then *generate* that spectrum.
        debug = kwargs.get("debug", False)
        def nlp(t, return_labels=False):
            # Apply fixed keywords.
            theta = dict(zip(parameters, t))
            theta.update(fixed)

            observed_variances = []
            observed_intensities = []
            for channel, spectrum in zip(matched_channels, data):

                # For each channel:
                # 1) Put the observed spectrum at rest.
                z = theta.get("z", theta.get("z_{}".format(channel), 0))
                rest_observed_disp = spectrum.disp * (1 - z)

                # 2) Correct for continuum.
                j, coeffs = 0, []
                while theta.get("continuum_{0}_{1}".format(channel, j), None):
                    coeffs.append(theta["continuum_{0}_{1}".format(channel, j)])
                    j += 1

                # Remember: model.__call__ calculates continuum based on the
                # *observed* wavelength points, so here we do the same (e.g.,
                # not those that have potentially been corrected for redshift)
                if not coeffs: continuum = 1.0
                else:
                    continuum = np.polyval(coeffs[::-1], spectrum.disp) 
                
                # 3) Put the observed data onto the self.wavelengths scale.
                # [TODO] do the resampling correctly.
                rebinned_observed_intensities = np.interp(self.wavelengths,
                    rest_observed_disp, spectrum.flux / continuum,
                    left=np.nan, right=np.nan)

                rebinned_observed_variances = np.interp(self.wavelengths,
                    rest_observed_disp, spectrum.variance / continuum,
                    left=np.nan, right=np.nan)

                observed_variances.append(rebinned_observed_variances)
                observed_intensities.append(rebinned_observed_intensities)

            # [TODO] This may be the wrong thing to do.
            observed_variances \
                = np.nanmean(np.vstack(observed_variances), axis=0)
            observed_intensities \
                = np.nanmean(np.vstack(observed_intensities), axis=0)

            # Solve for the astrophysical parameters.
            try:
                labels = self._solve_labels(observed_intensities,
                    observed_variances)

            except:
                logger.exception("Could not determine labels:")
                if debug: raise
                return np.inf

            if return_labels:
                return labels

            # Put the optimal astrophysical parameters into our theta.
            theta.update(labels)

            # Make the log-probability call.
            # Note: The inference._ln_probability call takes a dictionary.
            return -inference._ln_probability(theta, self, data, debug,
                matched_channels=matched_channels)

        # Do the optimisation.
        p0 = np.array([initial_theta[p] for p in parameters])
        x_opt = optimise.minimise(nlp, p0, **op_kwargs)

        # Put the result into a usable form.
        result = dict(zip(parameters, x_opt))
        result.update(nlp(x_opt, return_labels=True))
        result.update(fixed)
        
        result = OrderedDict([(k, result[k]) \
            for k in self.parameters if k in result])

        if full_output:
            # Create model fluxes and calculate some metric.
            chi_sq, dof, model_fluxes = self._chi_sq(result, data)
            result = (result, chi_sq, dof, model_fluxes)

        # Remove any prepared convolution functions.
        self._destroy_convolution_functions()

        return result