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
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.