code: def apply_correlation(self, sites, imt, residuals, stddev_intra=0):
          
          # 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))
docs: 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.