docs: def apply_correlation(self, sites, imt, residuals, stddev_intra):
          """
          Apply correlation to randomly sampled residuals.
  
          See Parent function
          """
          # stddev_intra is repeated if it is only 1 value for all the residuals
          if stddev_intra.shape[0] == 1:
              stddev_intra = numpy.matlib.repmat(
                  stddev_intra, len(sites.complete), 1)
          # Reshape 'stddev_intra' if needed
          stddev_intra = stddev_intra.squeeze()
          if not stddev_intra.shape:
              stddev_intra = stddev_intra[None]
  
          if self.uncertainty_multiplier == 0:   # No uncertainty
  
              # residuals were sampled from a normal distribution with
              # stddev_intra standard deviation. 'residuals_norm' are residuals
              # normalized, sampled from a standard normal distribution.
              # For this, every row of 'residuals' (every site) is divided by its
              # corresponding standard deviation element.
              residuals_norm = residuals / stddev_intra[sites.sids, None]
  
              # Lower diagonal of the Cholesky decomposition from/to cache
              try:
                  cormaLow = self.cache[imt]
              except KeyError:
                  # Note that instead of computing the whole correlation matrix
                  # corresponding to sites.complete, here we compute only the
                  # correlation matrix corresponding to sites.
                  cormaLow = numpy.linalg.cholesky(
                         numpy.diag(stddev_intra[sites.sids]) *
                         self._get_correlation_matrix(sites, imt) *
                         numpy.diag(stddev_intra[sites.sids]))
                  self.cache[imt] = cormaLow
  
              # Apply correlation
              return numpy.dot(cormaLow, residuals_norm)
  
          else:   # Variability (uncertainty) is included
              nsim = len(residuals[1])
              nsites = len(residuals)
  
              # Re-sample all the residuals
              residuals_correlated = residuals * 0
              for isim in range(0, nsim):
                  corma = self._get_correlation_matrix(sites, imt)
                  cov = (numpy.diag(stddev_intra[sites.sids]) * corma *
                         numpy.diag(stddev_intra[sites.sids]))
                  residuals_correlated[0:, isim] = (
                      numpy.random.multivariate_normal(
                          numpy.zeros(nsites), cov, 1))
  
              return residuals_correlated