docs: def read_run(filename, modelfn=None):
      """
      Read chain from a hdf5 saved with `save_run`.
  
      This function will also read the labels, data table, and metadata blobs
      stored in the original sampler. If you want to use the result object with
      `plot_fit` and setting the ``e_range`` parameter, you must provide the
      model function with the `modelfn` argument given that functions cannot be
      serialized in hdf5 files.
  
      Parameters
      ----------
      filename : str
          Filename of the hdf5 containing the chain, lnprobability, and blob
          arrays in the group 'sampler'
  
      modelfn : function, optional
          Model function to be attached to the returned sampler
  
      Returns
      -------
      result : class
          Container object with same properties as an `emcee.EnsembleSampler`
          resulting from a sampling run. This object can be passed onto
          `~naima.plot_fit`, `~naima.plot_chain`, and `~naima.plot_corner` for
          analysis as you would do with a `emcee.EnsembleSampler` instance.
      """
  
      # initialize empty sampler class to return
      result = _result()
      result.modelfn = modelfn
      result.run_info = {}
  
      f = h5py.File(filename, "r")
      # chain and lnprobability
      result.chain = np.array(f["sampler/chain"])
      result.lnprobability = np.array(f["sampler/lnprobability"])
  
      # blobs
      result.blobs = []
      shape = result.chain.shape[:2]
      blobs = []
      blobrank = []
      for i in range(100):
          # first read each of the blobs and convert to Quantities
          try:
              ds = f["sampler/blob{0}".format(i)]
              rank = np.ndim(ds[0])
              blobrank.append(rank)
              if rank <= 1:
                  blobs.append(u.Quantity(ds.value, unit=ds.attrs["unit"]))
              else:
                  blob = []
                  for j in range(np.ndim(ds[0])):
                      blob.append(
                          u.Quantity(
                              ds.value[:, j, :],
                              unit=ds.attrs["unit{0}".format(j)],
                          )
                      )
                  blobs.append(blob)
          except KeyError:
              break
  
      # Now organize in an awful list of lists of arrays
      for step in range(shape[1]):
          steplist = []
          for walker in range(shape[0]):
              n = step * shape[0] + walker
              walkerblob = []
              for j in range(len(blobs)):
                  if blobrank[j] <= 1:
                      walkerblob.append(blobs[j][n])
                  else:
                      blob = []
                      for k in range(blobrank[j]):
                          blob.append(blobs[j][k][n])
                      walkerblob.append(blob)
              steplist.append(walkerblob)
          result.blobs.append(steplist)
  
      # run info
      result.run_info = dict(f["sampler"].attrs)
      result.acceptance_fraction = f["sampler"].attrs["acceptance_fraction"]
      # labels
      result.labels = []
      for i in range(result.chain.shape[2]):
          result.labels.append(f["sampler"].attrs["label{0}".format(i)])
  
      # data
      data = Table(np.array(f["sampler/data"]))
      data.meta.update(f["sampler/data"].attrs)
      for col in data.colnames:
          if f["sampler/data"].attrs[col + "unit"] != "None":
              data[col].unit = f["sampler/data"].attrs[col + "unit"]
      result.data = QTable(data)
  
      return result