Repository: equinor/libres (language: python)

Issue:
Row scale
**Issue**
Part of: #1070


**Approach**
This PR is an initial refactoring to start support the localization described in #1070. The two first commits are 100% pure shuffling of code.


Discussion / hints:
Can one of the admins verify this patch?
Is it defensible to add new C code? We're trying to move away from the legacy code and use Python as much as possible.

In new code, we should:
1. Avoid `util_abort`. Use `std::optional` or error codes. Errors must be recoverable.
2. Avoid C getters. They account for a lot of headache when reading code. Prefer public structs. (This is a recent change of direction and can be discussed how to implement in practice)
3. Have consistent code style. The spacing between functions is inconsistent as well as the space between parens in function calls. You have both `foo( 1, 2 )` and `foo(1, 2)`. There's no standard at this time but at least choose one or the other for this PR.
4. Don't put double-underscore functions in public headers. These functions aren't public and shouldn't be available publicly.
5. Use `black` for new Python code
6. Use `pytest` for new or refactored Python unit tests
> Is it defensible to add new C code? We're trying to move away from the legacy code and use Python as much as possible.

I see that point. I will not argue the case either way - but for the task I have been assigned I really only see as a viable short term solution to proceed with C; if you prefer to let refactoring to Python proceed further before this problem is addressed I can back off for now.

> Avoid util_abort. Use std::optional or error codes. Errors must be recoverable.

I have used `std::exception` instead of `util_abort()`. I can look into using `std::optional<>` instead - but that certainly shifts the responsability for error checking. I don't think there is error checking code which can handle error returns from the very inner loops I am targeting here - but I will certainly keep it in mind.

> Avoid C getters. They account for a lot of headache when reading code. Prefer public structs.

I will do as told - but this is certainly a big change; for good or for bad most of libres is built with private access the way I have started on also with the new `row_scaling_type`.

> Have consistent code style. 

Sure - I will try harder!

> Don't put double-underscore functions in public headers. These functions aren't public and shouldn't be available publicly.

Well - there are patterns and there are patterns; one of the patterns in `libecl` and `libres` is that functions which take a `(void * arg)` argument and subsequently cast and call the real function have a trailing  `__` - i.e. like this:

```C++

void xxx_type_someting( xxx_type * arg ) {
    ....
}

void xxx_type_something__(void * arg) {
       xxx_type * xxx = xxx_type_safe_cast( arg );
       xxx_type_something( xxx );
}
```
The `xxx_type_something__()` function is to typically used to register a suitable `free( )` function in the `vector` and `hash` containers. I am not arguing the virtues of this particular pattern, but right now the quite some `__` functions need to be public to stick with that pattern.

> Use black [and Pytest] for new Python code 

Sure!


I will probably not make any code updates before January - but in general I will eventually try my best to address your comments.



> I see that point. I will not argue the case either way - but for the task I have been assigned I really only see as a viable short term solution to proceed with C; if you prefer to let refactoring to Python proceed further before this problem is addressed I can back off for now.

Yep, I agree. :+1: 

> I have used `std::exception` instead of `util_abort()`. I can look into using `std::optional<>` instead - but that certainly shifts the responsability for error checking. I don't think there is error checking code which can handle error returns from the very inner loops I am targeting here - but I will certainly keep it in mind.

Imo, even just returning a `bool` is fine. The `util_abort`s are caused programmer errors so it's debatable whether it makes sense to return anything. I think a `std::logic_error` or `std::argument_error` as you've used other places is fine.

> 
> > Avoid C getters. They account for a lot of headache when reading code. Prefer public structs.
> 
> I will do as told - but this is certainly a big change; for good or for bad most of libres is built with private access the way I have started on also with the new `row_scaling_type`.

This one's a bit tricky and I'll need to investigate how to properly interface the C code with the Python code without getters. I'll come up with an alternative.

> > Don't put double-underscore functions in public headers. These functions aren't public and shouldn't be available publicly.
> 
> Well - there are patterns and there are patterns; one of the patterns in `libecl` and `libres` is that functions which take a `(void * arg)` argument and subsequently cast and call the real function have a trailing `__` - i.e. like this:

Oh okay. I read `__` as private, but I guess it's fine as long as there's a use for this function.

> I will probably not make any code updates before January - but in general I will eventually try my best to address your comments.

Thanks so much for you effort!


Thank you for the review comments - I will address them although it might take some time. I would like to "make a statement" to the general question of C vs C++, homebrew containers vs stl and so on:

1. I have aimed to write code in the existing C style with homemade containers and so on[1]  - not in any way because I consider it superior, but because that is the style of the existing code.
2. If you have wishes in this regard - e.g. use stl for christ sake or something similar I am all ears, but I am a bit hesitant to devise my own solutions in this regard as I consider it highly improbable that I will hit your sweet spot.

[1]: At least in the API.
OK - I I now consider this ready for review. Qualitative remarks:

1. I have implemented the new class `row_scaling` as a C++ class and added a C layer for Python interop. For this particular case this was quite simple.
2. I have tried to address all of @dotfloat's comments - with the exception of public struct - which is maybe not particularlt relevant any longer.

When this is merged the next step is to actually start using this code, that will probably induce additional changes.
Looks good from my side in terms of the C++ stuff. :+1: 
Test this please!
Overall looks good. Some minor comments. Please tag me and @autt here if there is any questions/comments :)
@ManInFez & @autt : I have addressed all the comments except the one about automatic resizing; changes are force pushed.
Test this please!
@joakim-hove Very nice! 

Fix patch:
diff --git a/lib/enkf/enkf_main.cpp b/lib/enkf/enkf_main.cpp
index 6385f6b99..330a373a7 100644
--- a/lib/enkf/enkf_main.cpp
+++ b/lib/enkf/enkf_main.cpp
@@ -517,8 +517,7 @@ typedef struct {
 
 
 static void serialize_node( enkf_fs_type * fs ,
-                            const ensemble_config_type * ensemble_config,
-                            const char * key ,
+                            const enkf_config_node_type * config_node,
                             int iens ,
                             int report_step ,
                             int row_offset ,
@@ -526,7 +525,6 @@ static void serialize_node( enkf_fs_type * fs ,
                             const active_list_type * active_list,
                             matrix_type * A) {
 
-  const enkf_config_node_type * config_node = ensemble_config_get_node( ensemble_config , key );
   enkf_node_type * node = enkf_node_alloc( config_node );
   node_id_type node_id = {.report_step = report_step, .iens = iens  };
   enkf_node_serialize( node , fs , node_id , active_list , A , row_offset , column);
@@ -536,14 +534,13 @@ static void serialize_node( enkf_fs_type * fs ,
 
 static void * serialize_nodes_mt( void * arg ) {
   serialize_info_type * info = (serialize_info_type *) arg;
-  int iens;
-  for (iens = info->iens1; iens < info->iens2; iens++) {
+  const auto * node_info = info->node_info;
+  const enkf_config_node_type * config_node = ensemble_config_get_node( info->ensemble_config , node_info->key );
+  for (int iens = info->iens1; iens < info->iens2; iens++) {
     int column = int_vector_iget( info->iens_active_index , iens);
     if (column >= 0) {
-      const auto * node_info = info->node_info;
       serialize_node( info->src_fs ,
-                      info->ensemble_config,
-                      node_info->key ,
+                      config_node,
                       iens ,
                       info->report_step ,
                       node_info->row_offset ,
@@ -627,15 +624,13 @@ static int enkf_main_serialize_dataset( const ensemble_config_type * ens_config
 
 static void deserialize_node( enkf_fs_type * target_fs,
                               enkf_fs_type * src_fs,
-                              const ensemble_config_type * ensemble_config,
-                              const char * key ,
+                              const enkf_config_node_type * config_node,
                               int iens,
                               int target_step ,
                               int row_offset ,
                               int column,
                               const active_list_type * active_list,
                               matrix_type * A) {
-  const enkf_config_node_type * config_node = ensemble_config_get_node( ensemble_config , key );
 
   node_id_type node_id = {.report_step = target_step, .iens = iens  };
   enkf_node_type * node = enkf_node_alloc( config_node );
@@ -653,13 +648,12 @@ static void deserialize_node( enkf_fs_type * target_fs,
 
 static void * deserialize_nodes_mt( void * arg ) {
   serialize_info_type * info = (serialize_info_type *) arg;
-  int iens;
-  for (iens = info->iens1; iens < info->iens2; iens++) {
+  const auto * node_info = info->node_info;
+  const enkf_config_node_type * config_node = ensemble_config_get_node(info->ensemble_config, node_info->key);
+  for (int iens = info->iens1; iens < info->iens2; iens++) {
     int column = int_vector_iget( info->iens_active_index , iens );
-    if (column >= 0) {
-      const auto * node_info = info->node_info;
-      deserialize_node( info->target_fs , info->src_fs, info->ensemble_config , node_info->key , iens , info->target_step , node_info->row_offset , column, node_info->active_list , info->A );
-    }
+    if (column >= 0)
+      deserialize_node( info->target_fs , info->src_fs, config_node, iens , info->target_step , node_info->row_offset , column, node_info->active_list , info->A );
   }
   return NULL;
 }
@@ -1145,12 +1139,17 @@ static void enkf_main_analysis_update( enkf_main_type * enkf_main ,
             for (int ikw=0; ikw < scaled_keys.size(); ikw++) {
               const auto& key = scaled_keys[ikw];
               const active_list_type * active_list = local_dataset_get_node_active_list(dataset, key.c_str());
+              const auto matrix_rows = matrix_get_rows(A);
+              const auto* config_node = ensemble_config_get_node(serialize_info->ensemble_config, key.c_str());
+              const auto node_size = enkf_config_node_get_data_size(config_node, serialize_info->report_step);
+              if (matrix_get_rows(A) < node_size)
+                  matrix_resize( A , node_size , active_ens_size , false);
+
               for (int iens = 0; iens < bool_vector_size(ens_mask); iens++) {
                 int column = int_vector_iget( serialize_info->iens_active_index , iens);
                 if (column >= 0) {
                   serialize_node(serialize_info->src_fs,
-                                 serialize_info->ensemble_config,
-                                 key.c_str(),
+                                 config_node,
                                  iens,
                                  serialize_info->report_step,
                                  0,
@@ -1172,8 +1171,7 @@ static void enkf_main_analysis_update( enkf_main_type * enkf_main ,
                 if (column >= 0) {
                   deserialize_node(serialize_info->target_fs,
                                    serialize_info->src_fs,
-                                   serialize_info->ensemble_config,
-                                   key.c_str(),
+                                   ensemble_config_get_node(serialize_info->ensemble_config, key.c_str()),
                                    iens,
                                    serialize_info->target_step,
                                    0,


Test patch:
diff --git a/tests/res/enkf/test_row_scaling_case.py b/tests/res/enkf/test_row_scaling_case.py
index 0fc4d843c..1cd59d503 100644
--- a/tests/res/enkf/test_row_scaling_case.py
+++ b/tests/res/enkf/test_row_scaling_case.py
@@ -18,18 +18,22 @@
 import os
 import math
 import numpy as np
+import tempfile
+import shutil
 
 from tests import ResTest
+from tests.utils import tmpdir
 
 from ecl.util.util import BoolVector
+from ecl.grid import EclGridGenerator
 
 from res.test import ErtTestContext
-from res.enkf import RowScaling
+from res.enkf import RowScaling, ResConfig
 from res.enkf import EnkfNode, NodeId
 from res.enkf import ESUpdate
 from res.enkf import FieldConfig, ErtRunContext
 from res.enkf.enums import RealizationStateEnum
-
+from res.enkf import EnKFMain
 
 # This function will initialize the data in the case before the actual row
 # scaling test can be performed. The function will initialize the PORO field
@@ -51,6 +55,7 @@
 def init_data(main):
     fsm = main.getEnkfFsManager()
     init_fs = fsm.getFileSystem("init")
+    grid = main.eclConfig().getGrid()
 
     # Model: bhp = poro * 1000
     poro_mean = 0.15
@@ -68,7 +73,7 @@ def init_data(main):
         with open("fields/poro{}.grdecl".format(i), "w") as f:
             poro = random.gauss(poro_mean, poro_std)
             f.write("PORO")
-            for i in range(200):
+            for i in range(grid.get_num_active()):
                 if i % 10 == 0:
                     f.write("\n")
 
@@ -388,8 +393,7 @@ def test_row_scaling_using_assign_vector(self):
     # where the same field is updated in both steps. Because the
     # obs_data_allocE() function uses random state it is difficult to get
     # identical results from one ministep updating everything and two ministeps
-    # updating different parts of the field. The final assert is therefor
-    # unfortunately only approximative.
+    # updating different parts of the field.
     def test_2ministep(self):
         with ErtTestContext("row_scaling", self.config_file) as tc:
             main = tc.getErt()
@@ -456,5 +460,56 @@ def test_2ministep(self):
                 for iv, v1, v2 in zip(init_field, field1, field2):
                     assert iv != v1
 
-                    diff = abs(v1 - v2)
-                    assert diff < 0.05
+    # In the enkf_main_update() routine the A matrix is allocated with a
+    # default size; and should grow when exposed to a larger node. At a stage
+    # there was a bug in code for combination of row_scaling and
+    # matrix_resize(). The purpose of this test is to ensure that we create a
+    # sufficiently large node to invoke the rescaling.
+    @tmpdir()
+    def test_large_case(self):
+        with open("config", "w") as fp:
+            fp.write(
+                """NUM_REALIZATIONS 10
+GRID             CASE.EGRID
+FIELD            PORO    PARAMETER    poro.grdecl INIT_FILES:fields/poro%d.grdecl
+SUMMARY          WBHP
+OBS_CONFIG       observations.txt
+TIME_MAP timemap.txt
+"""
+            )
+
+        for f in ["timemap.txt", "observations.txt"]:
+            src_file = self.createTestPath(os.path.join("local/row_scaling", f))
+            shutil.copy(src_file, "./")
+        # The grid size must be greater than 250000 (the default matrix size in
+        # enkf_main_update())
+        grid = EclGridGenerator.create_rectangular((70, 70, 70), (1, 1, 1))
+        grid.save_EGRID("CASE.EGRID")
+        res_config = ResConfig(user_config_file="config")
+        main = EnKFMain(res_config)
+        init_fs = init_data(main)
+
+        # Configure the local updates
+        local_config = main.getLocalConfig()
+        local_config.clear()
+        local_data = local_config.createDataset("LOCAL")
+        local_data.addNode("PORO")
+        obs = local_config.createObsdata("OBSSET_LOCAL")
+        obs.addNode("WBHP0")
+        ministep = local_config.createMinistep("MINISTEP_LOCAL")
+        ministep.attachDataset(local_data)
+        ministep.attachObsset(obs)
+        updatestep = local_config.getUpdatestep()
+        updatestep.attachMinistep(ministep)
+
+        # Apply the row scaling
+        row_scaling = local_data.row_scaling("PORO")
+        ens_config = main.ensembleConfig()
+        poro_config = ens_config["PORO"]
+        field_config = poro_config.getFieldModelConfig()
+        grid = main.eclConfig().getGrid()
+        row_scaling.assign(field_config.get_data_size(), ScalingTest(grid))
+        es_update = ESUpdate(main)
+        update_fs = main.getEnkfFsManager().getFileSystem("target2")
+        run_context = ErtRunContext.ensemble_smoother_update(init_fs, update_fs)
+        es_update.smootherUpdate(run_context)