docs: def _uneven_utility_transform(systematic_utilities,
                                alt_IDs,
                                rows_to_alts,
                                shape_params,
                                intercept_params,
                                intercept_ref_pos=None,
                                *args, **kwargs):
      """
      Parameters
      ----------
      systematic_utilities : 1D ndarray.
          All elements should be ints, floats, or longs. Should contain the
          systematic utilities of each observation per available alternative.
          Note that this vector is formed by the dot product of the design matrix
          with the vector of utility coefficients.
      alt_IDs : 1D ndarray.
          All elements should be ints. There should be one row per obervation per
          available alternative for the given observation. Elements denote the
          alternative corresponding to the given row of the design matrix.
      rows_to_alts : 2D scipy sparse matrix.
          There should be one row per observation per available alternative and
          one column per possible alternative. This matrix maps the rows of the
          design matrix to the possible alternatives for this dataset. All
          elements should be zeros or ones.
      shape_params : None or 1D ndarray.
          If an array, each element should be an int, float, or long. There
          should be one value per shape parameter of the model being used.
      intercept_params : None or 1D ndarray.
          If an array, each element should be an int, float, or long. If J is the
          total number of possible alternatives for the dataset being modeled,
          there should be J-1 elements in the array.
      intercept_ref_pos : int, or None, optional.
          Specifies the index of the alternative, in the ordered array of unique
          alternatives, that is not having its intercept parameter estimated (in
          order to ensure identifiability). Should only be None if
          `intercept_params` is None.
  
      Returns
      -------
      transformed_utilities : 2D ndarray.
          Should have shape `(systematic_utilities.shape[0], 1)`. The returned
          array contains the transformed utility values for this model. All
          elements will be ints, longs, or floats.
      """
      # Convert the shape parameters back into their 'natural parametrization'
      natural_shapes = np.exp(shape_params)
      natural_shapes[np.isposinf(natural_shapes)] = max_comp_value
      # Figure out what shape values correspond to each row of the
      # systematic utilities
      long_natural_shapes = rows_to_alts.dot(natural_shapes)
  
      # Get the exponentiated neative utilities
      exp_neg_utilities = np.exp(-1 * systematic_utilities)
  
      # Get the log of 1 + exponentiated negative utilities
      log_1_plus_exp_neg_utilitiles = np.log1p(exp_neg_utilities)
  
      # Guard against overflow. Underflow not a problem since we add one to a
      # near zero number and log of one will evaluate to zero
      inf_idx = np.isinf(log_1_plus_exp_neg_utilitiles)
      log_1_plus_exp_neg_utilitiles[inf_idx] = -1 * systematic_utilities[inf_idx]
  
      # Get the exponentiated (negative utilities times the shape parameter)
      exp_neg_shape_utilities = np.exp(-1 *
                                       long_natural_shapes *
                                       systematic_utilities)
  
      # Get the log of 1 + exponentiated (negative utiltiies times the shape)
      log_1_plus_exp_neg_shape_utilities = np.log1p(exp_neg_shape_utilities)
  
      ##########
      # Guard against overflow
      ##########
      # Check for any values which have gone off to positive infinity
      inf_idx = np.isinf(log_1_plus_exp_neg_shape_utilities)
      # Replace those values with an approximation of the true values by ignoring
      # the "1." The idea is that 1 + infinity ~ infinity so the effect of the +1
      # on the log is minimal.
      if np.any(inf_idx):
          log_1_plus_exp_neg_shape_utilities[inf_idx] =\
                -1 * long_natural_shapes[inf_idx] * systematic_utilities[inf_idx]
  
      # Calculate the transformed utility values
      transformed_utilities = (systematic_utilities +
                               log_1_plus_exp_neg_utilitiles -
                               log_1_plus_exp_neg_shape_utilities)
      # Perform a final guard against numbers that are too large to deal with
      transformed_utilities[np.isposinf(transformed_utilities)] = max_comp_value
      transformed_utilities[np.isneginf(transformed_utilities)] = -max_comp_value
      transformed_utilities[np.isneginf(systematic_utilities)] = -max_comp_value
  
      # Account for the outside intercept parameters if there are any.
      if intercept_params is not None and intercept_ref_pos is not None:
          # Get a list of all the indices (or row indices) corresponding to the
          # alternatives whose intercept parameters are being estimated.
          needed_idxs = range(rows_to_alts.shape[1])
          needed_idxs.remove(intercept_ref_pos)
  
          if len(intercept_params.shape) > 1 and intercept_params.shape[1] > 1:
              # Get an array of zeros with shape
              # (num_possible_alternatives, num_parameter_samples)
              all_intercepts = np.zeros((rows_to_alts.shape[1],
                                         intercept_params.shape[1]))
              # For alternatives having their intercept estimated, replace the
              # zeros with the current value of the estimated intercepts
              all_intercepts[needed_idxs, :] = intercept_params
          else:
              # Get an array of zeros with shape (num_possible_alternatives,)
              all_intercepts = np.zeros(rows_to_alts.shape[1])
              # For alternatives having their intercept estimated, replace the
              # zeros with the current value of the estimated intercepts
              all_intercepts[needed_idxs] = intercept_params
  
          # Add the intercept values to f(x, beta, c)
          transformed_utilities += rows_to_alts.dot(all_intercepts)
  
      # Be sure to return a 2D array since other functions will be expecting that
      if len(transformed_utilities.shape) == 1:
          transformed_utilities = transformed_utilities[:, np.newaxis]
  
      return transformed_utilities