correlation.py 1.8 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2
import os
from scipy import integrate, interpolate
3 4
import matplotlib as mpl
mpl.use("agg")
Dion Haefner's avatar
Dion Haefner committed
5 6 7 8 9 10 11 12 13 14 15 16
import matplotlib.pyplot as plt
import numpy as np
import h5py

from dorie.testtools.utilities.decorators import evaluation_script
from dorie.testtools.utilities import statistics
from dorie.utilities.text_to_bool import text_to_bool

@evaluation_script
def evaluate(iniinfo,runtime):
    h5_path = iniinfo["generator.fft.outputPath"] + "/YField.h5"
    with h5py.File(h5_path, 'r') as f:
17
        y_data = np.array(f.get("/YField"), dtype=np.float64)
Dion Haefner's avatar
Dion Haefner committed
18 19 20 21

    print("Operating on file: {}".format(h5_path))

    # Calculate autocorrelation and correlation lengths
Dion Haefner's avatar
Dion Haefner committed
22 23
    ac = statistics.autocorrelation(y_data)
    rank = y_data.ndim
24
    dx = np.array(iniinfo["generator.extensions"].split(), dtype=np.float) / y_data.shape
Dion Haefner's avatar
Dion Haefner committed
25
    assert(len(dx) == rank)
Dion Haefner's avatar
Dion Haefner committed
26
    cl = [integrate.simps(ac,dx=dxi,axis=i).flat[0] for i, dxi in enumerate(dx)][::-1]
Dion Haefner's avatar
Dion Haefner committed
27 28 29 30 31 32 33

    # Compare to given values
    cl_ini = iniinfo["generator.fft.correlationLengths"].split()
    cl_ini = list(map(float,cl_ini))
    if "_correlation.abstol" in iniinfo:
        abstol = float(iniinfo["_correlation.abstol"])
    else:
Dion Haefner's avatar
Dion Haefner committed
34
        abstol = 1E-4
Dion Haefner's avatar
Dion Haefner committed
35 36 37
    if "_correlation.reltol" in iniinfo:
        reltol = float(iniinfo["_correlation.reltol"])
    else:
Dion Haefner's avatar
Dion Haefner committed
38
        reltol = 0.1
Dion Haefner's avatar
Dion Haefner committed
39 40 41 42 43 44 45 46 47 48
    passed = np.allclose(cl,cl_ini,rtol=reltol,atol=abstol)

    # Print output
    print("Measured correlation lengths: " + ", ".join(["{:.2e}".format(c) for c in cl]))
    print("Expected correlation lengths: " + ", ".join(["{:.2e}".format(c) for c in cl_ini]))

    if passed:
        print("- pass\n")
    else:
        print("- fail\n")
Dion Haefner's avatar
Dion Haefner committed
49 50
    if rank == 2:
        plt.imshow(ac,norm=mpl.colors.LogNorm())
Dion Haefner's avatar
Dion Haefner committed
51 52 53
        plt.colorbar()
        plt.savefig("ac-{}.png".format(iniinfo["__name"].split("/")[-1]))
    return bool(passed)