docs: def apply_correlation(self, sites, imt, residuals, stddev_intra=0):
          """
          Apply correlation to randomly sampled residuals.
  
          :param sites:
              :class:`~openquake.hazardlib.site.SiteCollection` residuals were
              sampled for.
          :param imt:
              Intensity measure type object, see :mod:`openquake.hazardlib.imt`.
          :param residuals:
              2d numpy array of sampled residuals, where first dimension
              represents sites (the length as ``sites`` parameter) and
              second one represents different realizations (samples).
          :param stddev_intra:
              Intra-event standard deviation array. Note that different sites do
              not necessarily have the same intra-event standard deviation.
          :returns:
              Array of the same structure and semantics as ``residuals``
              but with correlations applied.
  
          NB: the correlation matrix is cached. It is computed only once
          per IMT for the complete site collection and then the portion
          corresponding to the sites is multiplied by the residuals.
          """
          # intra-event residual for a single relization is a product
          # of lower-triangle decomposed correlation matrix and vector
          # of N random numbers (where N is equal to number of sites).
          # we need to do that multiplication once per realization
          # with the same matrix and different vectors.
          try:
              corma = self.cache[imt]
          except KeyError:
              corma = self.get_lower_triangle_correlation_matrix(
                  sites.complete, imt)
              self.cache[imt] = corma
          if len(sites.complete) == len(sites):
              return numpy.dot(corma, residuals)
          # it is important to allocate little memory, this is why I am
          # accumulating below; if S is the length of the complete sites
          # the correlation matrix has shape (S, S) and the residuals (N, s),
          # where s is the number of samples
          return numpy.sum(corma[sites.sids, sid] * res
                           for sid, res in zip(sites.sids, residuals))