docs: 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