docs: def create_matrix_block_indices(row_to_obs):
      """
      Parameters
      ----------
      row_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).
  
      Returns
      -------
      output_indices : list of arrays.
          There will be one array per column in `row_to_obs`. The array will note
          which rows correspond to which observations.
      """
      # Initialize the list of index arrays to be returned
      output_indices = []
      # Determine the number of observations in the dataset
      num_obs = row_to_obs.shape[1]
      # Get the indices of the non-zero elements and their values
      row_indices, col_indices, values = scipy.sparse.find(row_to_obs)
      # Iterate over each observation, i.e. each column in row_to_obs, and
      # determine which rows belong to that observation (i.e. the rows with ones
      # in them).
      for col in xrange(num_obs):
          # Store the array of row indices belonging to the current observation
          output_indices.append(row_indices[np.where(col_indices == col)])
  
      return output_indices