Site_Factors._interpolate() is not trilinear, and its result depends on the order of the 8 corner points #105

Found while reviewing #74. Pre-existing on main (b764057); #74 does not change the behavior.

Problem

Site_Factors._interpolate() calls scipy.interpolate.griddata(..., method='linear') on the 8 corners of the (Vs30, z1, PGA) cell around the query point. griddata() splits the corners into tetrahedra (a Delaunay triangulation) and interpolates linearly inside one tetrahedron. That is not trilinear interpolation on the cell. _get_results() uses it for every amplification and phase-shift factor:

This may be intended; if so, the docstring does not say so. It only says "High-dimensional interpolation" and refers to the griddata() documentation:

The facts:

• For f = x * y * z at the center of a cube, _interpolate() returns 3000.0; trilinear interpolation gives 3375.0.
• Passing the same 8 corners in a different order changes the result at the center to 4500.0. (_locate_grids() always passes them in the same order, so the library's results are repeatable.)
• On real data, the Fourier amplification of Site_Factors(265, 128, 0.25) differs from trilinear interpolation on the same 8 corners by up to 0.1505 (absolute) and 6.92% (relative).

test_interpolate() cannot tell the two methods apart; that is item 2 of #100.

import itertools
import numpy as np
from scipy.interpolate import RegularGridInterpolator
from PySeismoSoil.class_site_factors import Site_Factors as SF

# 1. f = x * y * z on one cube, value at the center
x, y, z = [1, 2], [10, 20], [100, 200]
corners = list(itertools.product(x, y, z))
f = lambda a, b, c: a * b * c
trilinear = RegularGridInterpolator((x, y, z), f(*np.meshgrid(x, y, z, indexing='ij')))
print(SF._interpolate(corners, [[f(*c)] for c in corners], (1.5, 15, 150)), trilinear((1.5, 15, 150)))
# [3000.] 3375.0

# 2. The same 8 corners in another order give another value
reordered = [corners[i] for i in (0, 1, 2, 3, 5, 6, 4, 7)]
print(SF._interpolate(reordered, [[f(*c)] for c in reordered], (1.5, 15, 150)))  # [4500.]

# 3. Real data: SF(265, 128, 0.25), Fourier amplification, vs. trilinear on the same 8 corners
sf = SF(265, 128, 0.25)
amp = sf.get_amplification().spectrum
grid = [sorted({SF.Vs30_array[c[0]] for c in sf._locate_grids()}),
        sorted({SF.z1_array[c[1]] for c in sf._locate_grids()}),
        sorted({SF.PGA_array[c[2]] for c in sf._locate_grids()})]
values = np.zeros((2, 2, 2, len(amp)))
for (i, vs30), (j, z1), (k, pga) in itertools.product(*(enumerate(g) for g in grid)):
    values[i, j, k] = SF._query('amplif', vs30, z1, pga, data_dir=sf.dir_amplif)[1][: len(amp)]
amp_tri = RegularGridInterpolator(grid, values)((265, 128, 0.25))
print(grid)  # [[250, 300], [75, 150], [0.2, 0.3]]
print(round(np.max(np.abs(amp - amp_tri)), 4), round(np.max(np.abs(amp - amp_tri) / amp_tri), 4))
# 0.1505 0.0692

Fix

1. If the tetrahedral interpolation is intended, say so in the _interpolate() docstring: the result is piecewise linear over a Delaunay split of the 8 corners, it is not trilinear, and it depends on the order of the corners. Then change test_interpolate() as described in Four tests cannot fail for the problem they are named after #100 item 2.
2. If trilinear interpolation is intended, replace the griddata() loop with one RegularGridInterpolator on the 2 x 2 x 2 cell, with all frequencies as the last axis (as in part 3 of the repro). This result does not depend on the order of the corners. With this change patched in, the existing tests in tests/test_class_site_factors.py and tests/test_class_site_effect_adjustment.py still pass. Add a test with f = x * y * z that expects 3375.0 at the cube center for both corner orders above.