docs: def compute_pmap_stats(pmaps, stats, weights, imtls):
      """
      :param pmaps:
          a list of R probability maps
      :param stats:
          a sequence of S statistic functions
      :param weights:
          a list of ImtWeights
      :param imtls:
          a DictArray of intensity measure types
      :returns:
          a probability map with S internal values
      """
      sids = set()
      p0 = next(iter(pmaps))
      L = p0.shape_y
      for pmap in pmaps:
          sids.update(pmap)
          assert pmap.shape_y == L, (pmap.shape_y, L)
      if len(sids) == 0:
          raise ValueError('All empty probability maps!')
      sids = numpy.array(sorted(sids), numpy.uint32)
      nstats = len(stats)
      curves = numpy.zeros((len(pmaps), len(sids), L), numpy.float64)
      for i, pmap in enumerate(pmaps):
          for j, sid in enumerate(sids):
              if sid in pmap:
                  curves[i, j] = pmap[sid].array[:, 0]
      out = p0.__class__.build(L, nstats, sids)
      for imt in imtls:
          slc = imtls(imt)
          w = [weight[imt] for weight in weights]
          for i, array in enumerate(compute_stats(curves[:, :, slc], stats, w)):
              for j, sid in numpy.ndenumerate(sids):
                  out[sid].array[slc, i] = array[j]
      return out