docs: def calc_hessian(beta,
                   design,
                   alt_IDs,
                   rows_to_obs,
                   rows_to_alts,
                   utility_transform,
                   transform_first_deriv_c,
                   transform_first_deriv_v,
                   transform_deriv_alpha,
                   block_matrix_idxs,
                   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.
          There 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 ndarray.
          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 ndarray.
          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.
      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 array
          whose elements contain the derivative of the tranformed utilities with
          respect to the vector of shape parameters. The dimensions of the
          returned vector should be `(design.shape[0], num_alternatives)`.
      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 utilities
          with respect to the vector of shape parameters. The dimensions of the
          returned vector should be `(design.shape[0], num_alternatives - 1)`. If
          `intercept_params == None`, the callable should return None.
      block_matrix_idxs : list of arrays.
          There will be one array per observation. The arrays will note which
          rows correspond to which observations.
      intercept_params : 1D ndarray.
          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 in the dataset.
      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.
          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
      -------
      hess : 2D ndarray.
          It's shape is `(beta.shape[0], 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 = np.ones(design.shape[0])
  
      ##########
      # Get the required matrices
      ##########
      # Differentiate the transformed utilities with respect to the shape params
      # Note that dh_dc will be a 2D scipy sparse matrix
      dh_dc = transform_first_deriv_c(sys_utilities, alt_IDs,
                                      rows_to_alts, shape_params)
      # Differentiate the transformed utilities with respect to the systematic
      # utilities. Note that dh_dv will be a 2D scipy sparse matrix.
      dh_dv = transform_first_deriv_v(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 utility
      # coefficients. Note that dh_db will be a 2D dense numpy matrix
      dh_db = dh_dv.dot(design)
  
      # Differentiate the probabilities with respect to the transformed utilities
      # Note that dp_dh will be a 2D dense numpy array
      block_matrices = create_matrix_blocks(long_probs, block_matrix_idxs)
      dp_dh = block_diag(*block_matrices) * weights[None, :]
  
      ##########
      # Calculate the second and mixed partial derivatives within the hessian
      ##########
      # Calculate the second derivative of the log-likelihood with repect to the
      # utility coefficients. Should have shape (design.shape[0],
      # design.shape[0]). Since dp_dh is a 2D dense numpy array and dh_db is a
      # 2D dense numpy matrix, using the .dot() syntax should work to compute
      # the dot product.
      d2_ll_db2 = -1 * dh_db.T.dot(dp_dh.dot(dh_db))
  
      ##########
      # Form and return the hessian
      ##########
      if shape_params is not None and intercept_params is not None:
          # Calculate the second derivative of the log-likelihood with respect
          # to the shape parameters. Should have shape (shape_params.shape[0],
          # shape_params.shape[0]). Note that since dp_dh is a 2D dense numpy
          # array and dh_dc is a sparse matrix or dense numpy matrix, to
          # compute the dot product we need to use the * operator
          d2_ll_dc2 = -1 * dh_dc.T.dot(dp_dh * dh_dc)
  
          # Calculate the second derivative of the log-likelihood with respect
          # to the intercept parameters. Should have shape (J - 1, J - 1) where
          # J is the total number of observed alternatives for this dataset.
          # Note that since dp_dh is a 2D dense numpy array and dh_d_alpha is a
          # sparse matrix or dense numpy matrix, to compute the dot product
          # we need to use the * operator
          d2_ll_d_alpha2 = -1 * dh_d_alpha.T.dot(dp_dh * dh_d_alpha)
  
          # Calculate the mixed second derivative of the log-likelihood with
          # respect to the intercept and shape parameters. Should have shape
          # (dh_d_alpha.shape[1], dh_dc.shape[1]). Note that since dp_dh is a 2D
          # dense numpy array and dh_dc is a sparse matrix or dense numpy
          # matrix, to compute the dot product we need to use the * operator
          d2_ll_dc_d_alpha = -1 * dh_d_alpha.T.dot(dp_dh * dh_dc)
  
          # Calculate the mixed partial derivative of the log-likelihood with
          # respect to the utility coefficients and then with respect to the
          # shape parameters. Should have shape (design.shape[0],
          # shape_params.shape[0]). Note that since dp_dh is a 2D dense numpy
          # array and dh_dc is a sparse matrix or dense numpy matrix, to
          # compute the dot product we need to use the * operator
          d2_ll_dc_db = -1 * dh_db.T.dot(dp_dh * dh_dc)
  
          # Calculate the mixed partial derivative of the log-likelihood with
          # respect to the utility coefficients and then with respect to the
          # intercept parameters. Should have shape (design.shape[0],
          # intercept_params.shape[0]). Note that since dp_dh is a 2D dense numpy
          # array and dh_d_alpha is a sparse matrix or dense numpy matrix, to
          # compute the dot product we need to use the * operator
          d2_ll_d_alpha_db = -1 * dh_db.T.dot(dp_dh * dh_d_alpha)
  
          # Form the 3 by 3 partitioned hessian of 2nd derivatives
          top_row = np.concatenate((d2_ll_dc2,
                                    d2_ll_dc_d_alpha.T,
                                    d2_ll_dc_db.T), axis=1)
          middle_row = np.concatenate((d2_ll_dc_d_alpha,
                                       d2_ll_d_alpha2,
                                       d2_ll_d_alpha_db.T), axis=1)
          last_row = np.concatenate((d2_ll_dc_db,
                                     d2_ll_d_alpha_db,
                                     d2_ll_db2), axis=1)
          hess = np.concatenate((top_row,
                                 middle_row,
                                 last_row), axis=0)
  
      elif shape_params is not None and intercept_params is None:
          # Calculate the second derivative of the log-likelihood with respect
          # to the shape parameters. Should have shape (shape_params.shape[0],
          # shape_params.shape[0]). Note that since dp_dh is a 2D dense numpy
          # array and dh_dc is a sparse matrix or dense numpy matrix, to
          # compute the dot product we need to use the * operator
          d2_ll_dc2 = -1 * dh_dc.T.dot(dp_dh * dh_dc)
  
          # Calculate the mixed partial derivative of the log-likelihood with
          # respect to the utility coefficients and then with respect to the
          # shape parameters. Should have shape (design.shape[0],
          # shape_params.shape[0]). Note that since dp_dh is a 2D dense numpy
          # array and dh_dc is a sparse matrix or dense numpy matrix, to
          # compute the dot product we need to use the * operator
          d2_ll_dc_db = -1 * dh_db.T.dot(dp_dh * dh_dc)
  
          hess = np.concatenate((np.concatenate((d2_ll_dc2,
                                                 d2_ll_dc_db.T), axis=1),
                                 np.concatenate((d2_ll_dc_db,
                                                 d2_ll_db2), axis=1)), axis=0)
  
      elif shape_params is None and intercept_params is not None:
          # Calculate the second derivative of the log-likelihood with respect
          # to the intercept parameters. Should have shape (J - 1, J - 1) where
          # J is the total number of observed alternatives for this dataset.
          # Note that since dp_dh is a 2D dense numpy array and dh_d_alpha is a
          # sparse matrix or dense numpy matrix, to compute the dot product
          # we need to use the * operator
          d2_ll_d_alpha2 = -1 * dh_d_alpha.T.dot(dp_dh * dh_d_alpha)
  
          # Calculate the mixed partial derivative of the log-likelihood with
          # respect to the utility coefficients and then with respect to the
          # intercept parameters. Should have shape (design.shape[0],
          # intercept_params.shape[0]). Note that since dp_dh is a 2D dense numpy
          # array and dh_d_alpha is a sparse matrix or dense numpy matrix, to
          # compute the dot product we need to use the * operator
          d2_ll_d_alpha_db = -1 * dh_db.T.dot(dp_dh * dh_d_alpha)
  
          hess = np.concatenate((np.concatenate((d2_ll_d_alpha2,
                                                 d2_ll_d_alpha_db.T), axis=1),
                                 np.concatenate((d2_ll_d_alpha_db,
                                                 d2_ll_db2), axis=1)), axis=0)
      else:
          hess = d2_ll_db2
  
      if ridge is not None:
          hess -= 2 * ridge
  
      # Make sure we are returning standard numpy arrays
      if isinstance(hess, np.matrixlib.defmatrix.matrix):
          hess = np.asarray(hess)
  
      return hess