docs: def calc_fisher_info_matrix(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.
          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.
      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 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 utility tranformation vector with respect to the
           vector of systematic utilities. The dimensions of the returned vector
           should be `(design.shape[0],  design.shape[0])`.
      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
      -------
      fisher_matrix : 2D numpy array.
          It will be a square matrix, with one row and one column for each shape,
          intercept, and index coefficient. Contains the BHHH approximation to
          the Fisher Information matrix of the log likelihood.
      """
      # Calculate the systematic utility for each alternative for each individual
      sys_utilities = design.dot(beta)
  
      # Calculate the probability that the individual associated with a given row
      # chooses the alternative associated with a given row.
      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])
      weights_per_obs =\
          np.max(rows_to_obs.toarray() * weights[:, None], axis=0)
  
      ##########
      # Get the required matrices
      ##########
      # Differentiate the transformed utilities with respect to the shape params
      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
      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. This should be a dense numpy array.
      dh_db = np.asarray(dh_dv.dot(design))
      # Differentiate the log likelihood w/ respect to the transformed utilities
      d_ll_dh = (choice_vector - long_probs)[np.newaxis, :]
  
      ##########
      # Create the matrix where each row represents the gradient of a particular
      # observations log-likelihood with respect to the shape parameters and
      # beta, depending on whether there are shape parameters being estimated
      ##########
      if shape_params is not None and intercept_params is not None:
          if isinstance(dh_dc, np.matrixlib.defmatrix.matrix):
              # Note that the '.A' transforms the matrix into a numpy ndarray
              gradient_vec = d_ll_dh.T * np.concatenate((dh_dc.A,
                                                         dh_d_alpha.toarray(),
                                                         dh_db), axis=1)
          else:
              gradient_vec = d_ll_dh.T * np.concatenate((dh_dc.toarray(),
                                                         dh_d_alpha.toarray(),
                                                         dh_db), axis=1)
      elif shape_params is not None and intercept_params is None:
          if isinstance(dh_dc, np.matrixlib.defmatrix.matrix):
              # Note that the '.A' transforms the matrix into a numpy ndarray
              gradient_vec = d_ll_dh.T * np.concatenate((dh_dc.A, dh_db), axis=1)
          else:
              gradient_vec = d_ll_dh.T * np.concatenate((dh_dc.toarray(),
                                                         dh_db), axis=1)
      elif shape_params is None and intercept_params is not None:
          # Note '.to_array()' is used because dh_d_alpha will be a sparse
          # matrix, but we want a 2D numpy array.
          gradient_vec = d_ll_dh.T * np.concatenate((dh_d_alpha.toarray(),
                                                     dh_db), axis=1)
      else:
          gradient_vec = d_ll_dh.T * dh_db
  
      # Make sure that we calculate the gradient PER OBSERVATION
      # and then take the outer products of those gradients.
      # Note that this is different than taking the outer products of the
      # gradient of the log-likelihood per available alternative per observation
      gradient_vec = rows_to_obs.T.dot(gradient_vec)
  
      # Compute and return the outer product of each row of the gradient
      # with itself. Then sum these individual matrices together. The line below
      # does the same computation just with less memory and time.
      fisher_matrix =\
          gradient_vec.T.dot(np.multiply(weights_per_obs[:, None], gradient_vec))
  
      if ridge is not None:
          # The rational behind adding 2 * ridge is that the fisher information
          # matrix should approximate the hessian and in the hessian we add
          # 2 * ridge at the end. I don't know if this is the correct way to
          # calculate the Fisher Information in ridge regression models.
          fisher_matrix -= 2 * ridge
  
      return fisher_matrix