docs: def estimate(self, data, full_output=False, **kwargs):
          """
          Estimate the model parameters, given the data.
          """
  
          # Number of model comparisons can be specified in the configuration.
          num_model_comparisons = self._configuration.get("estimate", {}).get(
              "num_model_comparisons", self.grid_points.size)
          # If it's a fraction, we need to convert that to an integer.
          if 1 > num_model_comparisons > 0:
              num_model_comparisons *= self.grid_points.size
  
          # If the num_model_comparison is provided as a keyword argument, use it.
          num_model_comparisons = kwargs.pop("num_model_comparisons",
              int(num_model_comparisons))
  
          logger.debug("Number of model comparisons to make for initial estimate:"
              " {0}".format(num_model_comparisons))
          
          # Match the data to the model channels.
          matched_channels, missing_channels, ignore_parameters \
              = self._match_channels_to_data(data)
  
          logger.debug("Matched channels: {0}, missing channels: {1}, ignore "
              "parameters: {2}".format(matched_channels, missing_channels,
                  ignore_parameters))
  
          # Load the intensities
          t = time()
          s = self.grid_points.size/num_model_comparisons # step size
          grid_points = self.grid_points[::s]
          intensities = np.memmap(
              self._configuration["model_grid"]["intensities"], dtype="float32",
              mode="r", shape=(self.grid_points.size, self.wavelengths.size))[::s]
          logger.debug("Took {:.0f} seconds to load and slice intensities".format(
              time() - t))
          # Which matched, data channel has the highest S/N?
          # (This channel will be used to estimate astrophysical parameters)
          data, pixels_affected = self._apply_data_mask(data)
          median_snr = dict(zip(matched_channels,
              [np.nanmedian(spec.flux/(spec.variance**0.5)) for spec in data]))
          median_snr.pop(None, None) # Remove unmatched data spectra
  
          ccf_channel = self._configuration.get("settings", {}).get("ccf_channel",
              max(median_snr, key=median_snr.get))
          if ccf_channel not in matched_channels:
              logger.warn("Ignoring CCF channel {0} because it was not a matched"
                  " channel".format(ccf_channel))
              ccf_channel = max(median_snr, key=median_snr.get)
  
          logger.debug("Channel with peak SNR is {0}".format(ccf_channel))
  
          # Are there *any* continuum parameters in any matched channel?
          any_continuum_parameters = any(map(lambda s: s.startswith("continuum_"),
              set(self.parameters).difference(ignore_parameters)))
  
          # [TODO]: CCF MASK
          # [TODO]: Don't require CCF if we have only continuum parameters.
  
          z_limits = self._configuration["settings"].get("ccf_z_limits", None)
  
          theta = {} # Dictionary for the estimated model parameters.
          best_grid_index = None
          c = speed_of_light.to("km/s").value
          for matched_channel, spectrum in zip(matched_channels, data):
              if matched_channel is None: continue
  
              # Do we need todo cross-correlation for this channel?
              # We do if there are redshift parameters for this channel,
              # or if there is a global redshift or global continuum parameters
              #   and this channel is the highest S/N.
              if "z_{}".format(matched_channel) in self.parameters \
              or ((any_continuum_parameters or "z" in self.parameters) \
              and matched_channel == ccf_channel):
  
                  # Get the continuum degree for this channel.
                  continuum_degree = self._configuration["model"].get("continuum",
                      { matched_channel: -1 })[matched_channel]
  
                  logger.debug("Perfoming CCF on {0} channel with a continuum "
                      "degree of {1}".format(matched_channel, continuum_degree))
  
                  # Get model wavelength indices that match the data.
                  # get the points that are in the mask, and within the spectrum
                  # limits
  
                  # TODO: Make this CCF not model mask.
                  idx = np.where(self._model_mask() \
                      * (self.wavelengths >= spectrum.disp[0]) \
                      * (spectrum.disp[-1] >= self.wavelengths))[0]
  
                  v, v_err, R = spectrum.cross_correlate(
                      (self.wavelengths[idx], intensities[:, idx]),
                      #(self.wavelengths, intensities),
                      continuum_degree=continuum_degree, z_limits=z_limits)
  
                  # Identify the best point by the CCF peak.
                  best = np.nanargmax(R)
  
                  # Now, why did we do CCF in this channel? Which model parameters
                  # should be updated?
                  if "z_{}".format(matched_channel) in self.parameters:
                      theta["z_{}".format(matched_channel)] = v[best] / c
                  elif "z" in self.parameters: 
                      # If there is a global redshift, update it.
                      theta["z"] = v[best] / c
  
                      # Continuum parameters will be updated later, so that each
                      # channel is checked to see if it has the highest S/N,
                      # otherwise we might be trying to calculate continuum
                      # parameters when we haven't done CCF on the highest S/N
                      # spectra yet.
  
                  if matched_channel == ccf_channel:
                      # Update astrophysical parameters.
                      theta.update(dict(zip(grid_points.dtype.names,
                          grid_points[best])))
                      best_grid_index = best
  
          # If there are continuum parameters, calculate them from the best point.
          if any_continuum_parameters:
              for matched_channel, spectrum in zip(matched_channels, data):
                  if matched_channel is None: continue
  
                  # The template spectra at the best point needs to be
                  # redshifted to the data, and then continuum coefficients
                  # calculated from that.
  
                  # Get the continuum degree for this channel.
                  continuum_degree = self._configuration["model"].get("continuum",
                      { matched_channel: -1 })[matched_channel]
  
              
                  # Get model wavelength indices that match the data.
                  idx = np.clip(self.wavelengths.searchsorted(
                      [spectrum.disp[0], spectrum.disp[-1]]) + [0, 1],
                      0, self.wavelengths.size)
  
                  # Redshift and bin the spectrum.
                  z = theta.get("z_{}".format(matched_channel), theta.get("z", 0))
  
                  best_intensities \
                      = np.copy(intensities[best_grid_index, idx[0]:idx[1]]).flatten()
  
                  # Apply model mask.
                  model_mask = self._model_mask(self.wavelengths[idx[0]:idx[1]])
                  best_intensities[~model_mask] = np.nan
  
                  best_intensities = best_intensities * specutils.sample.resample(
                      self.wavelengths[idx[0]:idx[1]] * (1 + z), spectrum.disp)
                  
                  # Calculate the continuum coefficients for this channel.
                  continuum = spectrum.flux/best_intensities
                  finite = np.isfinite(continuum)
  
                  try:
                      coefficients = np.polyfit(
                          spectrum.disp[finite], continuum[finite], continuum_degree,
                          )#w=spectrum.ivariance[finite])
                  except np.linalg.linalg.LinAlgError:
                      logger.exception("Exception in initial polynomial fit")
                      coefficients = np.polyfit(spectrum.disp[finite], continuum[finite],
                          continuum_degree)
  
                  # They go into theta backwards. such that coefficients[-1] is
                  # continuum_{name}_0
                  theta.update(dict(zip(
                      ["continuum_{0}_{1}".format(matched_channel, i) \
                          for i in range(continuum_degree + 1)],
                      coefficients[::-1]
                  )))
  
          # Remaining parameters could be: resolving power, outlier pixels,
          # underestimated variance.
          remaining_parameters = set(self.parameters)\
              .difference(ignore_parameters)\
              .difference(theta)
  
          if remaining_parameters:
              logger.debug("Remaining parameters to estimate: {0}. For these we "
                  "will just assume reasonable initial values.".format(
                  remaining_parameters))
  
              for parameter in remaining_parameters:
                  if parameter == "resolution" \
                  or parameter.startswith("resolution_"):
  
                      if parameter.startswith("resolution_"):
                          spectra = [data[matched_channels.index(
                              parameter.split("_")[1])]]
                      else:
                          spectra = [s for s in data if s is not None]
  
                      R = [s.disp.mean()/np.diff(s.disp).mean() for s in spectra]
  
                      # Assume oversampling rate of ~5.
                      theta.update({ parameter: np.median(R)/5.})
  
                  elif parameter == "ln_f" or parameter.startswith("ln_f_"):
                      theta.update({ parameter: 0.5 }) # Not overestimated.
  
                  elif parameter in ("Po", "Vo"):
                      theta.update({
                          "Po": 0.01, # 1% outlier pixels.
                          "Vo": np.mean([np.nanmedian(s.variance) for s in data]),
                      })
  
          logger.info("Initial estimate: {}".format(theta))
          # Having full_output = True means return the best spectra estimate.
          if full_output:
  
              # Create model fluxes and calculate some metric.
              __intensities = np.copy(intensities[best_grid_index])
  
              # Apply model masks.
              __intensities[~self._model_mask()] = np.nan
  
              chi_sq, dof, model_fluxes = self._chi_sq(theta, data,
                  __intensities=__intensities, __no_precomputed_binning=True)
              del intensities
  
              return (theta, chi_sq, dof, model_fluxes)
  
          # Delete the reference to intensities
          del intensities
          return theta