docs: def _uneven_transform_deriv_shape(systematic_utilities,
                                    alt_IDs,
                                    rows_to_alts,
                                    shape_params,
                                    output_array=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.
      output_array : 2D scipy sparse array.
          The array should have shape `(systematic_utilities.shape[0],
          shape_params.shape[0])`. It's data is to be replaced with the correct
          derivatives of the transformation vector with respect to the vector of
          shape parameters. This argument is NOT optional.
  
      Returns
      -------
      output_array : 2D scipy sparse array.
          The shape of the returned array is `(systematic_utilities.shape[0],
          shape_params.shape[0])`. The returned array specifies the derivative of
          the transformed utilities with respect to the shape parameters. All
          elements are ints, floats, or longs.
      """
      # 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_shapes = rows_to_alts.dot(natural_shapes)
  
      # Get the exponentiated (utilities times the shape parameter)
      exp_shape_utilities = np.exp(long_shapes * systematic_utilities)
  
      # Calculate the derivative of h_ij with respect to shape_j.
      derivs = (systematic_utilities / (1.0 + exp_shape_utilities))
      # Guard against overflow. Only for cases of systematic_utilities becomming
      # huge. It is unlikely this safeguard will be needed.
      derivs[np.isposinf(systematic_utilities)] = 0
      # Guard against underflow from v --> -inf.
      huge_index = np.isneginf(systematic_utilities)
      derivs[huge_index] = -max_comp_value
  
      # Return the matrix of dh_dshapes. Note the matrix should be of dimension
      # (systematic_utilities.shape[0], shape_params.shape[0])
      # Note that the "* long_shapes" accounts for the fact that the derivative
      # of the natural shape parameters with resepect to the actual shape
      # parameters being estimated is simply
      # exp(actual shape parameters) = natural shape parameters. The
      # multiplication comes from the chain rule.
      output_array.data = derivs * long_shapes
      return output_array