code: def stack_nll(shape, components, ylims, weights=None):
          
          n_bins = shape[0]
          n_vals = shape[1]
  
          if weights is None:
              weights = np.ones((len(components)))
  
          norm_vals = np.zeros(shape)
          nll_vals = np.zeros(shape)
          nll_offsets = np.zeros((n_bins))
          for i in range(n_bins):
              log_min = np.log10(ylims[0])
              log_max = np.log10(ylims[1])
              norm_vals[i, 1:] = np.logspace(log_min, log_max, n_vals - 1)
              for c, w in zip(components, weights):
                  nll_vals[i] += w * c[i].interp(norm_vals[i]) - c.nll_offsets[i]
  
              # Reset the offsets
              nll_obj = LnLFn(norm_vals[i], nll_vals[i])
              ll_offset = nll_obj.fn_mle()
              nll_vals[i] -= ll_offset
              nll_offsets[i] = -ll_offset
              pass
  
          return norm_vals, nll_vals, nll_offsets
docs: Combine the log-likelihoods from a number of components.
  
          Parameters
          ----------
          shape    :  tuple
             The shape of the return array
  
          components : `~fermipy.castro.CastroData_Base`
             The components to be stacked
  
          weights : array-like
  
          Returns
          -------
          norm_vals : `numpy.ndarray`
             N X M array of Normalization values
  
          nll_vals  : `numpy.ndarray`
             N X M array of log-likelihood values
  
          nll_offsets :  `numpy.ndarray`
             N array of maximum log-likelihood values in each bin