Spatial Cross-Correlation in GMFs and Related Issues

[g-weatherill]

A recent PR (#3545) related to issue #3544 introduced spatial cross-correlation into the oq-engine for shakemap applications.

Briefly, spatial cross-correlation is a generalised form of spatial correlation for application when multiple IMTs are required and the correlations between IMTs at the same site and at different sites needs to be taken into account. When implementing this it is no longer possible to build or treat each IMT field separately in the simulation.

The implementation of spatial cross-correlation in #3545 hard-codes one particular approach. However, there are a growing number of models in the literature, some may need to be constructed differently from the approach implemented here, and depending on application it may be necessary to support more of them. Furthermore, implicit within a spatial cross-correlation model is the property that for two different IMTs located at the same site, the correlation should collapse to the known cross-correlation between IMTs. This can be a useful property to consider for many applications as it opens up the eventual possibility of being able to run vector hazard analyses.

To this end it may be beneficial to refactor the correlation model objects to support the more general cases when one wishes to consider spatial correlation, cross-correlation and spatial cross-correlation, whilst still permitting the possibility of running the calculation with only the spatial correlation and not the cross-correlation. This would need some additional changes to the GMFComputer to support this case.

Additionally:

1. For various reasons, many spatial cross-correlation models apply only to the total standard deviation, as they are derived by considering the total residual only and not by explicitly separating our the inter-event residual correlation. OpenQuake's current approach of raising an error when applying spatial correlation to a GMPE that does not have inter- an intra-event standard deviations, rather than just applying the spatial correlation to the total residual, was never actually necessary, as it would have not have been a problem to apply the correlation to the total residual, given that the inter-event residual is constant for a given IMT. This behaviour should be changed anyway and both spatial correlation and spatial cross-correlation should be permitted when the GMPE returns total standard deviation. Depending on the spatial cross-correlation model it should only apply to the total residual in any case.

2. Memory and computation is a big problem for spatial cross-correlation. At the moment the calculation of a spatially correlated field of ground motions requires the application of cholesky factorisation to the full correlation Nsites by Nsites correlation matrix. This holds true for spatial cross-correlation (other methods are available in the literature but each has various cost-benefits) but with the problem that the correlation matrix to which the cholesky factorisation applies is now on the order of (NIMTs * NSites) by (NIMTs * NSites), e.g. for 10,000 sites and 3 IMTs the spatial cross-correlation matrix is 30,000 by 30,000. It is easy to see that this will eventually cause computational problems with RAM.

3. For cholesky factorisation to work the correlation matrix must be positivie-definite. If the cross-correlation model has been developed correctly then this should be implicit; however, depending on the IMTs and the interpolations that are applied (or the sub-set of sites that are selected) it is possible to run into cases where positive definiteness cannot be guaranteed. An elegant way to deal with these cases would be necessary.

[micheles]

Thanks for this comment, Graeme. Do you imply that we should change/extend hazardlib too?

[g-weatherill]

Yes, I think the correlation model object in the hazardlib should be improved to try to address all of the correlation requirements for both hazard and shakemap applications

[mmpagani]

Thanks, Graeme for the message.

I think that the fundamental problem behind your note stands on the fact that in the OQ engine we never worked at abstracting the expected behaviour of a GM correlation model. We can work on defining proper abstract classes for the different correlation models and therefore add more flexibility (and accommodate new models).

Regarding point 2 and 3, we should explore the possibility of parallelising the Cholesky decomposition.

[g-weatherill]

Further thoughts (from attempting implementation):

1. General structure of the BaseCorrelationModel is not necessarily poorly suited to spatial cross-correlation. The main adaptation is the apply_correlation method needs adapting such that it is fed with a list of IMTs rather than a single IMT. The critical distinction here is that the full spatial cross-correlation matrix should be constructed and factorised on the first call, and then on subsequent calls the necessary rows/columns for the sites can then be sliced out. This is working for me with imts as a list and residuals (as input) as samples from the standard Normal distribution of dimension [len(imts) * len(sites), num_events]:

    def apply_correlation(self, sites, imts, residuals):
        """
        Applies the correlation model to the sites
        """
        if not self.cache["corma"]:
            # For first implementation then get the lower
            # matrix
            self.get_lower_triangle_correlation_matrix(sites.complete, imts)
        nsites = len(sites)
        if len(sites.complete) == nsites:
            # No filtering of sites
            corr_residuals = np.matmul(self.cache["corma"], residuals)
        else:
            # Need to locate indices of correlation matrix corresponding
            # to the selected sites
            idx = np.tile(sites.indices, len(imts))
            lb = 0
            ub = len(sites)
            norig_sites = len(sites.complete)
            for i in range(len(imts)):
                idx[lb:ub] += (i * norig_sites) 
                lb += nsites
                ub += nsites
            corr_residuals = np.matmul(self.cache["corma"][np.ix_(idx, idx)], residuals)
       
        # corr_residuals is 2-D array [Nimts * Nsites, Nevents]
        # need to re-arrange back to 3-D shape
        residuals = np.empty([len(imts), nsites, corr_residuals.shape[1]])
        idx = np.arange(nsites)
        for i in range(len(imts)):
            residuals[i] = corr_residuals[idx, :]
            idx += nsites
        return residuals

2. Changes are needed in the GmfComputer and further upstream. If one were to assume all IMTS can be calculated from a desired GSIM then IMTs can be set upon instantiation of the GmfComputer and the loop over the IMTs pushed further down into _compute(...). In the case of no correlation the output would be identical to that in the current form. The problem comes if one were to need to calculated different IMTs from different GSIMS (e.g. Sa(1.0) from GMPE1, PGV from GMPE 2). Leaving aside theoretical questions about this, there is precedent for this application and in some cases this is the only feasible approach. In this case it is better if the user defines not just a GMPE but a dictionary associating GMPEs to IMTs (or vice-versa), e.g.:

<logicTreeBranch branchID="ASB14+T03">
    <uncertaintyModel>{AkkarEtAlRjb2014: [PGA, PGV], TravasarouEtAl2003: [IA]}</uncertaintyModel>
    <uncertaintyWeight>1.0</uncertaintyWeight>
</logicTreeBranch>

In this case the GSIMs and their corresponding IMTs can be kept. Many upstream changes needed in the calculators.

[mmpagani]

Ciao Graeme, thanks for this. I'm fine with point 1. For the second point - unless there are specific and urgent requirements - I suggest proceeding more gradually. i.e. let's have first a working version supporting the simplest case. Later we can extend it to also include more articulated cases.

[g-weatherill]

As with #3833, work on this feature is abandoned for the time being. It may be re-visited in future. Closing the issue for now.