docs: def calc_gradient(beta,
                    design,
                    alt_IDs,
                    rows_to_obs,
                    rows_to_alts,
                    choice_vector,
                    utility_transform,
                    transform_first_deriv_c,
                    transform_first_deriv_v,
                    transform_deriv_alpha,
                    intercept_params,
                    shape_params,
                    ridge,
                    weights):
      """
      Parameters
      ----------
      beta : 1D ndarray.
          All elements should by ints, floats, or longs. Should have 1 element
          for each utility coefficient being estimated (i.e. num_features).
      design : 2D ndarray.
          Tjere should be one row per observation per available alternative.
          There should be one column per utility coefficient being estimated. All
          elements should be ints, floats, or longs.
      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_obs : 2D scipy sparse array.
          There should be one row per observation per available alternative and
          one column per observation. This matrix maps the rows of the design
          matrix to the unique observations (on the columns).
      rows_to_alts : 2D scipy sparse array
          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.
      choice_vector : 1D ndarray.
          All elements should be either ones or zeros. There should be one row
          per observation per available alternative for the given observation.
          Elements denote the alternative which is chosen by the given
          observation with a 1 and a zero otherwise.
      utility_transform : callable.
          Must accept a 1D array of systematic utility values, a 1D array of
          alternative IDs, and miscellaneous args and kwargs. Should return a 1D
          array whose elements contain the appropriately transformed systematic
          utility values, based on the current model being evaluated.
      transform_first_deriv_c : callable.
          Must accept a 1D array of systematic utility values, a 1D array of
          alternative IDs, the `rows_to_alts` array, (shape parameters if there
          are any) and miscellaneous args and kwargs. Should return a 2D matrix
          or sparse array whose elements contain the derivative of the tranformed
          utility vector with respect to the vector of shape parameters. The
          dimensions of the returned vector should be
          `(design.shape[0], num_alternatives)`. If there are no shape parameters
          then the callable should return None.
      transform_first_deriv_v : callable.
          Must accept a 1D array of systematic utility values, a 1D array of
          alternative IDs, (shape parameters if there are any) and miscellaneous
          args and kwargs. Should return a 2D array whose elements contain the
          derivative of the tranformed utility vector with respect to the vector
          of systematic utilities. The dimensions of the returned vector should
          be `(design.shape[0], design.shape[0])`.
      transform_deriv_alpha : callable.
          Must accept a 1D array of systematic utility values, a 1D array of
          alternative IDs, the `rows_to_alts` array, (intercept parameters if
          there are any) and miscellaneous args and kwargs. Should return a 2D
          array whose elements contain the derivative of the tranformed utility
          vector with respect to the vector of shape parameters. The dimensions
          of the returned vector should be
          `(design.shape[0], num_alternatives - 1)`. If there are no intercept
          parameters, the callable should return None.
      intercept_params : 1D numpy array or None.
          If an array, each element should be an int, float, or long. For
          identifiability, there should be J- 1 elements where J is the total
          number of observed alternatives for this dataset. Default == None.
      shape_params : 1D ndarray or None.
         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.
         Default == None.
      ridge : int, float, long, or None.
          Determines whether or not ridge regression is performed. If an int,
          float or long is passed, then that scalar determines the ridge penalty
          for the optimization. Default = None.
      weights : 1D ndarray or None.
          Allows for the calculation of weighted log-likelihoods. The weights can
          represent various things. In stratified samples, the weights may be
          the proportion of the observations in a given strata for a sample in
          relation to the proportion of observations in that strata in the
          population. In latent class models, the weights may be the probability
          of being a particular class.
  
      Returns
      -------
      gradient : 1D ndarray.
         It's shape is (beta.shape[0], ). It is the second derivative of the log-
         likelihood with respect to beta.
      """
      # Calculate the systematic utility for each alternative for each individual
      sys_utilities = design.dot(beta)
  
      # Calculate the probability of each individual choosing each available
      # alternative for that individual.
      long_probs = calc_probabilities(beta,
                                      design,
                                      alt_IDs,
                                      rows_to_obs,
                                      rows_to_alts,
                                      utility_transform,
                                      intercept_params=intercept_params,
                                      shape_params=shape_params,
                                      return_long_probs=True)
  
      # Calculate the weights for the sample
      if weights is None:
          weights = 1
  
      ##########
      # Get the required matrices
      ##########
      # Differentiate the transformed utilities with respect to the shape params
      # Note that dh_dc should be a sparse array
      dh_dc = transform_first_deriv_c(sys_utilities, alt_IDs,
                                      rows_to_alts, shape_params)
      # Differentiate the transformed utilities by the intercept params
      # Note that dh_d_alpha should be a sparse array
      dh_d_alpha = transform_deriv_alpha(sys_utilities, alt_IDs,
                                         rows_to_alts, intercept_params)
      # Differentiate the transformed utilities with respect to the systematic
      # utilities. Note that dh_dv should be a sparse matrix
      dh_dv = transform_first_deriv_v(sys_utilities, alt_IDs,
                                      rows_to_alts, shape_params)
      # Differentiate the transformed utilities with respect to the utility
      # coefficients. Note that dh_db should be a dense **matrix**, not a dense
      # 2D array. This is because the dot product of a 2D scipy sparse array and
      # a 2D dense numpy array yields a 2D dense numpy matrix
      dh_db = dh_dv.dot(design)
      # Differentiate the log likelihood w/ respect to the transformed utilities
      # Note that d_ll_dh will be a dense 2D numpy array.
      d_ll_dh = np.multiply(weights, choice_vector - long_probs)[np.newaxis, :]
  
      # Calculate the gradient of the log-likelihood with respect to the betas
      d_ll_d_beta = d_ll_dh.dot(dh_db)
  
      ##########
      # Form and return the gradient
      ##########
      if shape_params is not None and intercept_params is not None:
          # Note that we use d_ll_dh * dh_dc and d_ll_dh * dh_d_alpha because
          # that is how one computes the dot product between a dense 2D numpy
          # array and a 2D sparse matrix. This is due to numpy ndarrays and
          # scipy sparse matrices not playing nicely together. However, numpy
          # ndarrays and numpy matrices can be dot producted together,
          # hence d_ll_dh.dot(dh_db).
  
          # Note that the 'np.asarray' is because dll_dh * dh_dc will be a row
          # matrix, but we want a 1D numpy array.
          gradient = np.concatenate((np.asarray(d_ll_dh * hstack((dh_dc,
                                                                  dh_d_alpha),
                                                                 format='csr')),
                                     d_ll_d_beta), axis=1).ravel()
          params = np.concatenate((shape_params, intercept_params, beta),
                                  axis=0)
  
      elif shape_params is not None and intercept_params is None:
          # Note that we use d_ll_dh * dh_dc because that is how one computes
          # the dot product between a dense 2D numpy array and a 2D sparse matrix
          # This is due to numpy ndarrays and scipy sparse matrices not playing
          # nicely together. However, numpy ndarrays and numpy matrices can be
          # dot producted together, hence d_ll_dh.dot(dh_db).
  
          # Note that the 'np.asarray' is because dll_dh * dh_dc will be a row
          # matrix, but we want a 1D numpy array.
          gradient = np.concatenate((np.asarray(d_ll_dh * dh_dc), d_ll_d_beta),
                                    axis=1).ravel()
          params = np.concatenate((shape_params, beta), axis=0)
  
      elif shape_params is None and intercept_params is not None:
          # Note that we use d_ll_dh * dh_d_alpha because that's how one computes
          # the dot product between a dense 2D numpy array and a 2D sparse matrix
          # This is due to numpy ndarrays and scipy sparse matrices not playing
          # nicely together. However, numpy ndarrays and numpy matrices can be
          # dot producted together, hence d_ll_dh.dot(dh_db).
  
          # Note 'np.asarray' is used because dll_dh * dh_d_alpha will be a row
          # matrix, but we want a 1D numpy array.
          gradient = np.concatenate((np.asarray(d_ll_dh * dh_d_alpha),
                                     d_ll_d_beta), axis=1).ravel()
          params = np.concatenate((intercept_params, beta), axis=0)
  
      else:
          gradient = d_ll_d_beta.ravel()
          params = beta
  
      if ridge is not None:
          gradient -= 2 * ridge * params
  
      return gradient