docs: def create_matrix_blocks(long_probs, matrix_block_indices):
      """
      Parameters
      ----------
      long_probs : 1D ndarray.
          There should be one row per observation per available alternative. The
          elements of the array will indicate the probability of the alternative
          being the outcome associated with the corresponding observation.
      matrix_block_indices : list of arrays.
          There will be one array per observation. The arrays will note which
          rows correspond to which observations.
  
      Returns
      -------
      output_matrices : list of matrices.
          Each matrix will contain the derivative of P_i with respect to H_i, and
          there will be one matrix for each observations i. P_i is the array of
          probabilities of each observation being associated with its available
          alternatives. H_i is the array of transformed index values for each
          alternative that is available to observation i.
      """
      # Initialize the list of matrices that is to be returned.
      output_matrices = []
      # Iterate over each observation, i.e. over each list of rows that
      # corresponds to each observation.
      for indices in matrix_block_indices:
          # Isolate P_i, the vector of probabilities of each alternative that
          # is associated with the current observation
          current_probs = long_probs[indices]
          # Get the outer product of the current probabilities
          # probability_outer_product = np.outer(current_probs, current_probs)
          probability_outer_product = robust_outer_product(current_probs,
                                                           current_probs)
  
          # Create the desired dP_i/dh_i matrix
          dP_i_dh_i = np.diag(current_probs) - probability_outer_product
          # Ensure that the diagonal is positive and non-zero, since it must be.
          diag_idx = np.diag_indices_from(dP_i_dh_i)
          diag_values = dP_i_dh_i[diag_idx].copy()
          diag_values[diag_values == 0] = min_comp_value
          dP_i_dh_i[diag_idx] = diag_values
          # Guard against underflow on the off-diagonals
          underflow_idxs = np.where(dP_i_dh_i == 0)
          for i in xrange(underflow_idxs[0].size):
              row_idx, col_idx = underflow_idxs[0][i], underflow_idxs[1][i]
              if row_idx != col_idx:
                  # Since this type of underflow essentially comes from
                  # multiplying two very small numbers, just set the overall
                  # result to a small number
                  dP_i_dh_i[row_idx,
                            col_idx] = -1 * min_comp_value
          # Store the desired dP_i/dh_i matrix
          output_matrices.append(dP_i_dh_i)
  
      return output_matrices