Use active realization list and add difference calculation of GRF posterior minus GRF prior

Author: oddvarlia

src/subscript/field_statistics/field_statistics.py

@@ -1080,6 +1080,7 @@ def get_specifications(
         use_population_stdev = input_dict[key]
 
     use_geogrid_fields = False
+    geogrid_name = None
     zone_names_used = None
     zone_conformity = None
     zone_code_names = None
@@ -1339,6 +1340,11 @@ def calc_stats(
     if not use_geogrid_fields:
         return
 
+    # Get list of active realization (Must be active for all iterations in iter_list)
+    active_real, number_of_skipped = get_active_real(
+        iter_list, ens_path, nreal, geogrid_name
+    )
+
     ensemble_path = ens_path
 
     logger.info(f"Number of realizations: {nreal}")
@@ -1356,8 +1362,8 @@ def calc_stats(
                         (ertbox_size[0], ertbox_size[1], ertbox_size[2], nreal),
                         dtype=np.float32,
                     )
-                    number_of_skipped = 0
-                    for real_number in range(nreal):
+
+                    for real_number in active_real:
                         grid_dimensions, subgrids, property_param = (
                             read_ensemble_realization(
                                 ensemble_path,
@@ -1368,11 +1374,7 @@ def calc_stats(
                                 geogrid_name,
                             )
                         )
-                        if grid_dimensions is None:
-                            txt = f" Skip non-existing realization: {real_number}"
-                            logger.info(txt)
-                            number_of_skipped += 1
-                            continue
+                        assert grid_dimensions is not None
                         ertbox_prop_values = get_values_in_ertbox(
                             grid_dimensions,
                             subgrids,
@@ -1445,8 +1447,8 @@ def calc_stats(
                         (ertbox_size[0], ertbox_size[1], ertbox_size[2], nreal),
                         dtype=np.int32,
                     )
-                    number_of_skipped = 0
-                    for real_number in range(nreal):
+
+                    for real_number in active_real:
                         grid_dimensions, subgrids, property_param = (
                             read_ensemble_realization(
                                 ensemble_path,
@@ -1457,11 +1459,7 @@ def calc_stats(
                                 geogrid_name,
                             )
                         )
-
-                        if grid_dimensions is None:
-                            logger.info(f" Skip realization: {real_number}")
-                            number_of_skipped += 1
-                            continue
+                        assert grid_dimensions is not None
                         ertbox_prop_values = get_values_in_ertbox(
                             grid_dimensions,
                             subgrids,
@@ -1576,7 +1574,8 @@ def calc_temporary_field_stats(
     # Check if any need to continue to calculation
     if not use_temporary_fields:
         return
-
+    # Get list of active realization (Must be active for all iterations in iter_list)
+    active_real, number_of_skipped = get_active_real(iter_list, ens_path, nreal)
     # Import realizations of temporary field parameters
     for param_name in param_list:
         for iteration in iter_list:
@@ -1590,21 +1589,15 @@ def calc_temporary_field_stats(
                 (ertbox_size[0], ertbox_size[1], ertbox_size[2], nreal),
                 dtype=np.float32,
             )
-
-            number_of_skipped = 0
-            for real_number in range(nreal):
+            for real_number in active_real:
                 filepath = (
                     ens_path
                     / Path(
                         "realization-" + str(real_number) + "/iter-" + str(iteration)
                     )
                     / Path(full_param_filename)
                 )
-                if not filepath.exists():
-                    txt = f" Skip non-existing realization: {real_number}"
-                    logger.info(txt)
-                    number_of_skipped += 1
-                    continue
+                assert filepath.exists()
                 property = xtgeo.gridproperty_from_file(filepath, fformat="roff")
                 values = property.values
                 all_values[:, :, :, real_number] = values
@@ -1657,6 +1650,38 @@ def calc_temporary_field_stats(
                 xtgeo_ertbox_stdev.to_file(result_stdev_file_path, fformat="roff")
 
 
+def get_active_real(
+    iter_list: list, ens_path: Path, nreal: int, geogrid_name: str = ""
+):
+    """Get a list of active realizations"""
+    number_of_skipped = 0
+    active_real = []
+    for real_number in range(nreal):
+        real_exist = True
+        for iteration in iter_list:
+            ensemble_path = ens_path / Path(
+                "realization-" + str(real_number) + "/iter-" + str(iteration)
+            )
+            if len(geogrid_name) > 0:
+                grid_path = Path("share/results/grids/" + geogrid_name + ".roff")
+                file_path_grid = ensemble_path / grid_path
+                txt = f" Skip non-existing realization of grid: {real_number}"
+            else:
+                file_path_grid = ensemble_path
+                txt = f" Skip non-existing realization: {real_number}"
+            if not file_path_grid.exists():
+                logger.info(txt)
+                number_of_skipped += 1
+                # No need to check other iterations since active_real should
+                # only be those realizations that exists for all specified
+                # iterations in iter_list
+                real_exist = False
+                continue
+        if real_exist:
+            active_real.append(real_number)
+    return active_real, number_of_skipped
+
+
 def generate_script(
     rms_load_script, ert_config_path, result_path, field_stat_config_file
 ):
@@ -1853,25 +1878,47 @@ def main():
         if init_path and param_names:
             for param_name in param_names:
                 for iteration in iter_list:
-                    new_name = "mean_" + param_name + "_" + str(iteration)
-                    param_file_name = Path(result_path) / Path(new_name + ".roff")
-                    prop_param = xtgeo.gridproperty_from_file(
-                        param_file_name, fformat="roff"
+                    mean_name = "mean_" + param_name + "_" + str(iteration)
+                    std_name = "stdev_" + param_name + "_" + str(iteration)
+
+                    mean_param_file_name = Path(result_path) / Path(mean_name + ".roff")
+                    mean_prop_param = xtgeo.gridproperty_from_file(
+                        mean_param_file_name, fformat="roff"
                     )
-                    print(f"Read: {{new_name}} into {{GRIDNAME}}")
-                    if label:
-                        new_name = new_name + "_" + label
-                    prop_param.to_roxar(PRJ, GRIDNAME, new_name)
+                    print(f"Read: {{mean_name}} into {{GRIDNAME}}")
 
-                    new_name = "stdev_" + param_name + "_" + str(iteration)
-                    param_file_name = Path(result_path) / Path(new_name + ".roff")
-                    prop_param = xtgeo.gridproperty_from_file(
-                        param_file_name, fformat="roff"
+                    std_param_file_name = Path(result_path) / Path(std_name + ".roff")
+                    std_prop_param = xtgeo.gridproperty_from_file(
+                        std_param_file_name, fformat="roff"
                     )
-                    print(f"Read: {{new_name}} into {{GRIDNAME}}")
+                    print(f"Read: {{std_name}} into {{GRIDNAME}}")
+
                     if label:
-                        new_name = new_name + "_" + label
-                    prop_param.to_roxar(PRJ, GRIDNAME, new_name)
+                        new_mean_name = mean_name + "_" + label
+                        new_std_name = std_name + "_" + label
+                        difference_mean_name = "diff_" + new_mean_name
+                        difference_std_name = "diff_" + new_std_name
+
+                    mean_prop_param.to_roxar(PRJ, GRIDNAME, new_mean_name)
+                    std_prop_param.to_roxar(PRJ, GRIDNAME, new_std_name)
+
+                    if iteration == iter_list[0]:
+                        # Init
+                        mean_prop_param_init = mean_prop_param
+                        std_prop_param_init = std_prop_param
+                    elif iteration == iter_list[-1]:
+                        mean_prop_param_upd = mean_prop_param
+                        std_prop_param_upd = std_prop_param
+                        diff_mean_prop_param = mean_prop_param_upd.copy()
+                        diff_std_prop_param = std_prop_param_upd.copy()
+                        diff_mean_prop_param.values = \
+                            mean_prop_param_upd.values - mean_prop_param_init.values
+                        diff_std_prop_param.values = \
+                            std_prop_param_upd.values - std_prop_param_init.values
+                        diff_mean_prop_param.to_roxar(PRJ,
+                            GRIDNAME, difference_mean_name)
+                        diff_std_prop_param.to_roxar(PRJ,
+                            GRIDNAME, difference_std_name)
 
 if __name__ == "__main__":
     main()

tests/test_field_statistics.py

@@ -4,12 +4,14 @@
 from pathlib import Path
 
 import fmu.config.utilities as utils
+import gaussianfft as sim
 import numpy as np
 import pytest
 import xtgeo
 
 from subscript.field_statistics.field_statistics import (
     calc_stats,
+    calc_temporary_field_stats,
     check_disc_param_name_dict,
     check_param_name_dict,
     check_use_zones,
@@ -1101,3 +1103,125 @@ def test_main(tmp_path, config_file, config_dict, print_info=True):
     # For this test not to fail, the CONFIG_DICT and the specified
     # config file in yaml format must define the same setup
     assert compare_with_referencedata(ens_path, result_path, print_check=True)
+
+
+def simulate_ensembles(
+    nreal: int, nx: int, ny: int, nz: int, nreal_lost: int, result_dir: str
+) -> None:
+    sim.seed(123456)
+    variogram = sim.variogram("exponential", 100.0, 50.0, 5.0, 45.0, 0.0)
+    dx = 5.0
+    dy = 5.0
+    dz = 1.0
+    init_ens_field_3d = np.zeros((nx, ny, nz, nreal), dtype=np.float32)
+    final_ens_field_3d = np.zeros((nx, ny, nz, nreal - nreal_lost), dtype=np.float32)
+    for i in range(nreal):
+        gauss_vector = sim.simulate(variogram, nx, dx, ny, dy, nz, dz)
+        field_3d = gauss_vector.reshape((nx, ny, nz), order="F")
+        init_ens_field_3d[:, :, :, i] = field_3d
+        if i < (nreal - nreal_lost):
+            final_ens_field_3d[:, :, :, i] = field_3d
+    write_ensemble(init_ens_field_3d, final_ens_field_3d, result_dir)
+
+
+def write_ensemble(init_ens_field_3d, final_ens_field_3d, result_dir):
+    nx, ny, nz, nreal = init_ens_field_3d.shape
+    _, _, _, nreal_final = final_ens_field_3d.shape
+
+    init_real_defined = np.arange(nreal)
+    final_real_defined = np.arange(nreal_final)
+    for iter in [0, 1]:
+        if iter == 0:
+            for i in init_real_defined:
+                real_path = (
+                    Path(result_dir)
+                    / Path("realization-" + str(i))
+                    / Path("iter-" + str(iter))
+                )
+                real_path /= Path("rms/output/aps")
+                values = init_ens_field_3d[:, :, :, i]
+                name = "GRF"
+                xtgeo_param = xtgeo.GridProperty(
+                    ncol=nx,
+                    nrow=ny,
+                    nlay=nz,
+                    name=name,
+                    roxar_dtype=np.float32,
+                    values=values,
+                )
+                if not real_path.exists():
+                    real_path.mkdir(parents=True, exist_ok=True)
+                real_path /= Path("GRF.roff")
+                xtgeo_param.to_file(real_path, fformat="roff")
+        else:
+            for i in final_real_defined:
+                real_path = (
+                    Path(result_dir)
+                    / Path("realization-" + str(i))
+                    / Path("iter-" + str(iter))
+                )
+                values = final_ens_field_3d[:, :, :, i]
+                name = "GRF"
+                xtgeo_param = xtgeo.GridProperty(
+                    ncol=nx,
+                    nrow=ny,
+                    nlay=nz,
+                    name=name,
+                    roxar_dtype=np.float32,
+                    values=values,
+                )
+                if not real_path.exists():
+                    real_path.mkdir(parents=True, exist_ok=True)
+                real_path /= Path("GRF.roff")
+                xtgeo_param.to_file(real_path, fformat="roff")
+
+
+@pytest.mark.parametrize(
+    "nreal, nreal_lost",
+    [
+        (10, 0),
+        (10, 2),
+        (10, 8),  # Number of realizations must be at least 2
+    ],
+)
+def test_compare_mean_stdev_of_ensembles(tmp_path, nreal, nreal_lost):
+    nx = 5
+    ny = 5
+    nz = 2
+    field_name = "GRF"
+    simulate_ensembles(nreal, nx, ny, nz, nreal_lost, tmp_path)
+    ertbox_size = (nx, ny, nz)
+    input_dict = {
+        "nreal": nreal,
+        "iterations": [0, 1],
+        "use_population_stdev": False,
+        "temporary_ertbox_fields": {
+            "initial_relative_path": "rms/output/aps",
+            "parameter_names": [field_name],
+        },
+    }
+    ens_path = tmp_path
+    result_path = tmp_path
+    ert_config_path = None
+    print("Calculate mean and stdev of ensembles")
+    calc_temporary_field_stats(
+        input_dict, ens_path, result_path, ert_config_path, ertbox_size
+    )
+    compare_ensemble_stats(result_path, field_name)
+
+
+def compare_ensemble_stats(result_path, field_name, tolerance=1e-8):
+    print("Compare mean and stdev of two ensembles")
+    mean_file_name1 = Path(result_path) / Path("mean_" + field_name + "_0.roff")
+    mean_file_name2 = Path(result_path) / Path("mean_" + field_name + "_1.roff")
+    xtgeo_mean1_field = xtgeo.gridproperty_from_file(mean_file_name1, fformat="roff")
+    xtgeo_mean2_field = xtgeo.gridproperty_from_file(mean_file_name2, fformat="roff")
+    diff_values = np.abs(xtgeo_mean1_field.values - xtgeo_mean2_field.values)
+    assert np.all(diff_values < tolerance)
+
+    sdev_file_name1 = Path(result_path) / Path("stdev_" + field_name + "_0.roff")
+    sdev_file_name2 = Path(result_path) / Path("stdev_" + field_name + "_1.roff")
+    xtgeo_sdev1_field = xtgeo.gridproperty_from_file(sdev_file_name1, fformat="roff")
+    xtgeo_sdev2_field = xtgeo.gridproperty_from_file(sdev_file_name2, fformat="roff")
+    diff_values = np.abs(xtgeo_sdev1_field.values - xtgeo_sdev2_field.values)
+    assert np.all(diff_values < tolerance)