Source code for dorie.utilities.fuzzy_compare_grid

import numpy as np
from scipy import interpolate

from dorie.utilities.grid import BaseGrid

def _equal_grid_shapes(cells1,cells2):
    shapes_equal = cells1.shape == cells2.shape
    if not shapes_equal:
        print("Grid shapes do not match: {}, {}".format(cells1.shape,cells2.shape))
    return shapes_equal

def _equal_coordinates(cells1,cells2,abstol,reltol):
    if cells1.shape != cells2.shape:
        return False
    coords_equal = np.allclose(cells1,cells2,atol=abstol,rtol=reltol)
    if not coords_equal:
        print("Grids do not match")
        idiff = np.argmax(np.sqrt(np.sum((cells1-cells2)**2, axis=0)))
        print("\tLargest difference: got {} and {}".format(cells1[:,idiff],cells2[:,idiff]))
    return coords_equal

def _equal_variable_keys(grid1,grid2):
    data1, data2 = [g.data() for g in (grid1,grid2)]
    keys_equal = set(data1.keys()) == set(data2.keys())
    if not keys_equal:
        print("Variable keys do not match")
    return keys_equal

def _sort_cells(grid):
    cells = grid.cellCenters()
    ix = np.lexsort(cells.T)
    sorted_cells = cells[ix,...]
    return sorted_cells

def _sort_grid(grid):
    cells = grid.cellCenters()
    ix = np.lexsort(cells.T)
    sorted_cells = cells[ix,...]
    sorted_data = {key: val[ix] for key,val in grid.data().items()}
    return sorted_cells, sorted_data

def _equal_grids(grid1,grid2,abstol,reltol):
    cells1 = _sort_cells(grid1)
    cells2 = _sort_cells(grid2)
    return _equal_grid_shapes(cells1,cells2) and _equal_coordinates(cells1,cells2,abstol,reltol)

def _calc_residuals(grid1,grid2,interpolate_grid):
    cells1 = grid1.cellCenters()
    cells2 = grid2.cellCenters()
    data1 = grid1.data()
    data2 = grid2.data()

    if interpolate_grid:
        residuals = {}
        for key in data1.keys():
            print("Interpolating variable {} ...".format(key))
            interpolator = interpolate.NearestNDInterpolator(cells1,data1[key])
            residuals[key] = interpolator(cells2) - data2[key]
    else:
        _, sorted_data1 = _sort_grid(grid1)
        _, sorted_data2 = _sort_grid(grid2)
        residuals = {var: v1 - v2 for var, v1, v2 in zip(data1.keys(),sorted_data1.values(),sorted_data2.values())}
    return residuals, grid2

def _compare_max_norm(residuals, atol, rtol):
    passed = True
    for var, res in residuals.items():
        if not np.allclose(res, 0., atol=atol, rtol=rtol):
            print(var, res)
            print("Data in variable {} does not match".format(var))
            print("\tLargest difference: {}".format(np.max(np.abs(res))))
            passed = False
    return passed

def _compare_l2_norm(residuals, grid, abstol):
    passed = True
    for var, res in residuals.items():
        l2 = np.sqrt(np.sum(grid.integrate(res**2) / grid.integrate(1)))
        if np.all(l2 > abstol):
            print("Data in variable {} does not match".format(var))
            print("\tL2 error: {}".format(l2))
            passed = False
    return passed


def compare(grid1,grid2,abstol=1E-7,reltol=1E-7,strict=False):
    if not isinstance(grid1,BaseGrid) or not isinstance(grid2,BaseGrid):
        raise TypeError("Inputs must be of type BaseGrid (got: {} and {})"\
                        .format(type(grid1),type(grid2)))

    interpolate_grid = False
    if not _equal_grids(grid1,grid2,abstol,reltol):
        if strict:
            return False
        else:
            interpolate_grid = True

    if not _equal_variable_keys(grid1,grid2):
        return False

    residuals, res_grid = _calc_residuals(grid1,grid2,interpolate_grid)

    if strict:
        return _compare_max_norm(residuals, abstol, reltol)
    else:
        return _compare_l2_norm(residuals, res_grid, abstol)