Merge pull request #534 from xjules/resample_w_extrapolation

Support for extrapolation in the resample function.

lib/ecl/ecl_sum.cpp

@@ -779,37 +779,20 @@ const char * ecl_sum_get_general_var_unit( const ecl_sum_type * ecl_sum , const
 /*****************************************************************/
 
 
-ecl_sum_type * ecl_sum_alloc_resample(const ecl_sum_type * ecl_sum, const char * ecl_case, const time_t_vector_type * _times) {
-  time_t start_time = ecl_sum_get_data_start(ecl_sum);
-  time_t end_time = ecl_sum_get_end_time(ecl_sum);
-  time_t input_t0 = time_t_vector_get_first( _times );
-  time_t_vector_type * times = time_t_vector_alloc(0,0);
-
+ecl_sum_type * ecl_sum_alloc_resample(const ecl_sum_type * ecl_sum, const char * ecl_case, const time_t_vector_type * times, bool lower_extrapolation, bool upper_extrapolation) {
   /*
-    This is a temporary hack - the time argument used for resampling are
-    internally clamped to the [start, end] interval of the simulation; however
-    when outputting the original time values from the input argument _times is
-    used.
-
-    The real fix for this is to update the complete time interpolation code to
-    return something sensible when the time arguments are outside of simulation
-    time range.
+    If lower and  / or upper extrapolation is set to true it makes sure that resampling returns the first / last value of the simulation
+    or in the case of derivate / rate then it gets zero. if these are set to false, we jus throw exception
   */
-  for (int i=0; i < time_t_vector_size(_times); i++) {
-    time_t t = time_t_vector_iget(_times, i);
-
-    if (t < start_time)
-      t = start_time;
-
-    if (t > end_time)
-      t = end_time;
 
-    time_t_vector_append(times, t);
-  }
+  time_t start_time = ecl_sum_get_data_start(ecl_sum);
+  time_t end_time = ecl_sum_get_end_time(ecl_sum);
+  time_t input_start = time_t_vector_get_first( times );
+  time_t input_end = time_t_vector_get_last( times );
 
-  if ( time_t_vector_get_first(times) < start_time )
+  if ( !lower_extrapolation && input_start < start_time )
     return NULL;
-  if ( time_t_vector_get_last(times) > ecl_sum_get_end_time(ecl_sum) )
+  if ( !upper_extrapolation && input_end > end_time)
     return NULL;
   if ( !time_t_vector_is_sorted(times, false) )
     return NULL;
@@ -821,35 +804,52 @@ ecl_sum_type * ecl_sum_alloc_resample(const ecl_sum_type * ecl_sum, const char *
   if ( util_string_equal(smspec_node_get_unit(node), "DAYS" ) )
     time_in_days = true;
 
-  ecl_sum_type * ecl_sum_resampled = ecl_sum_alloc_writer( ecl_case , ecl_sum->fmt_case , ecl_sum->unified , ecl_sum->key_join_string , start_time , time_in_days , grid_dims[0] , grid_dims[1] , grid_dims[2] );
-
-
+  //create elc_sum_resampled with TIME node only
+  ecl_sum_type * ecl_sum_resampled = ecl_sum_alloc_writer( ecl_case , ecl_sum->fmt_case , ecl_sum->unified , ecl_sum->key_join_string , input_start , time_in_days , grid_dims[0] , grid_dims[1] , grid_dims[2] );
 
+  //add remaining nodes
   for (int i = 0; i < ecl_smspec_num_nodes(ecl_sum->smspec); i++) {
     const smspec_node_type * node = ecl_smspec_iget_node_w_node_index(ecl_sum->smspec, i);
-    if (util_string_equal(smspec_node_get_gen_key1(node), "TIME"))
-      continue;
-
+    if (util_string_equal(smspec_node_get_keyword(node), "TIME"))
+        continue;
     ecl_sum_add_smspec_node(  ecl_sum_resampled, node );
   }
 
   /*
     The SMSPEC header structure has been completely initialized, it is time to
     start filling it up with data.
+
   */
   ecl_sum_vector_type * ecl_sum_vector = ecl_sum_vector_alloc(ecl_sum, true);
   double_vector_type * data = double_vector_alloc( ecl_sum_vector_get_size(ecl_sum_vector) , 0);
 
-
   for (int report_step = 0; report_step < time_t_vector_size(times); report_step++) {
-    time_t t       = time_t_vector_iget(times, report_step);
-    time_t input_t = time_t_vector_iget(_times, report_step);
-
-    /* Look up interpolated data in the original case. */
-    ecl_sum_get_interp_vector( ecl_sum, t, ecl_sum_vector, data);
+    time_t input_t = time_t_vector_iget(times, report_step);
+    if (input_t < start_time) {
+      //clamping to the first value for t < start_time or if it is a rate than derivative is 0
+      for (int i=1; i < ecl_smspec_num_nodes(ecl_sum->smspec); i++) {
+        double value = 0;
+        const smspec_node_type * node = ecl_smspec_iget_node_w_node_index(ecl_sum->smspec, i);
+        if (!smspec_node_is_rate(node))
+            value = ecl_sum_iget_first_value(ecl_sum, smspec_node_get_params_index(node));
+        double_vector_iset(data, i-1, value);
+      }
+    } else if (input_t > end_time) {
+      //clamping to the last value for t > end_time or if it is a rate than derivative is 0
+      for (int i=1; i < ecl_smspec_num_nodes(ecl_sum->smspec); i++) {
+        double value = 0;
+        const smspec_node_type * node = ecl_smspec_iget_node_w_node_index(ecl_sum->smspec, i);
+        if (!smspec_node_is_rate(node))
+          value = ecl_sum_iget_last_value(ecl_sum, smspec_node_get_params_index(node));
+        double_vector_iset(data, i-1, value);
+      }
+    } else {
+      /* Look up interpolated data in the original case. */
+      ecl_sum_get_interp_vector( ecl_sum, input_t, ecl_sum_vector, data);
+    }
 
     /* Add timestep corresponding to the interpolated data in the resampled case. */
-    ecl_sum_tstep_type * tstep = ecl_sum_add_tstep( ecl_sum_resampled , report_step , input_t - input_t0);
+    ecl_sum_tstep_type * tstep = ecl_sum_add_tstep( ecl_sum_resampled , report_step , input_t - input_start);
     for (int data_index = 0; data_index < ecl_sum_vector_get_size(ecl_sum_vector); data_index++) {
       double value = double_vector_iget(data,data_index);
       int params_index = data_index + 1;  // The +1 shift is because the first element in the tstep is time value.
@@ -858,7 +858,6 @@ ecl_sum_type * ecl_sum_alloc_resample(const ecl_sum_type * ecl_sum, const char *
   }
   double_vector_free( data );
   ecl_sum_vector_free( ecl_sum_vector );
-  time_t_vector_free(times);
   return ecl_sum_resampled;
 }
 

lib/ecl/tests/ecl_sum_alloc_resampled_test.cpp

@@ -37,7 +37,7 @@ void test_correct_time_vector() {
   time_t_vector_append(t, util_make_date_utc( 4,1,2010 ));
   time_t_vector_append(t, util_make_date_utc( 6,1,2010 ));
   time_t_vector_append(t, util_make_date_utc( 8,1,2010 ));
-  ecl_sum_type * ecl_sum_resampled = ecl_sum_alloc_resample(ecl_sum, "kk", t);
+  ecl_sum_type * ecl_sum_resampled = ecl_sum_alloc_resample(ecl_sum, "kk", t, false, false);
   test_assert_int_equal(  ecl_sum_get_report_time(ecl_sum_resampled, 2)  , util_make_date_utc( 6,1,2010 ));
 
   const ecl_smspec_type * smspec_resampled = ecl_sum_get_smspec(ecl_sum_resampled);
@@ -57,23 +57,56 @@ void test_correct_time_vector() {
   ecl_sum_free(ecl_sum);
 }
 
+
+void test_resample_extrapolate_rate() {
+
+    ecl_sum_type * ecl_sum = test_alloc_ecl_sum();
+
+    time_t_vector_type * t = time_t_vector_alloc( 0 , 0 );
+    time_t_vector_append(t, util_make_date_utc( 1,1,2009 ));
+    time_t_vector_append(t, util_make_date_utc( 4,1,2010 ));
+    time_t_vector_append(t, util_make_date_utc( 12,1,2010 ));
+
+    ecl_sum_type * ecl_sum_resampled = ecl_sum_alloc_resample(ecl_sum, "kk", t, true, true);
+
+
+    const ecl_smspec_type * smspec_resampled = ecl_sum_get_smspec(ecl_sum_resampled);
+    const smspec_node_type * node1 = ecl_smspec_iget_node_w_params_index(smspec_resampled, 1);
+    const smspec_node_type * node3 = ecl_smspec_iget_node_w_params_index(smspec_resampled, 3);
+
+
+    //testing extrapolation for rate wrt. 3 dates: lower, inside and upper
+    test_assert_double_equal(0, ecl_sum_get_from_sim_time( ecl_sum_resampled, util_make_date_utc( 1, 1, 2009), node3));
+    test_assert_double_equal(10.000, ecl_sum_get_from_sim_time( ecl_sum_resampled, util_make_date_utc( 4, 1, 2010), node3));
+    test_assert_double_equal(0, ecl_sum_get_from_sim_time( ecl_sum_resampled, util_make_date_utc( 12, 1, 2010), node3));
+
+    //testing extrapolation for variable wrt. 3 dates: lower, inside and upper
+    test_assert_double_equal(0, ecl_sum_get_from_sim_time( ecl_sum_resampled, util_make_date_utc( 1, 1, 2009 ), node1) );
+    test_assert_double_equal(2.000, ecl_sum_get_from_sim_time( ecl_sum_resampled, util_make_date_utc( 4, 1,2010 ), node1) );
+    test_assert_double_equal(6.000, ecl_sum_get_from_sim_time( ecl_sum_resampled, util_make_date_utc( 12, 1,2010 ), node1) );
+
+
+    ecl_sum_free(ecl_sum_resampled);
+    time_t_vector_free(t);
+    ecl_sum_free(ecl_sum);
+}
+
 void test_not_sorted() {
   ecl_sum_type * ecl_sum = test_alloc_ecl_sum();
   time_t_vector_type * t = time_t_vector_alloc( 0 , 0 );
   time_t_vector_append(t, util_make_date_utc( 1,1,2010 ));
   time_t_vector_append(t, util_make_date_utc( 3,1,2010 ));
   time_t_vector_append(t, util_make_date_utc( 2,1,2010 ));
-  test_assert_NULL( ecl_sum_alloc_resample( ecl_sum, "kk", t) );
+  test_assert_NULL( ecl_sum_alloc_resample( ecl_sum, "kk", t, false, false) );
   time_t_vector_free(t);
   ecl_sum_free(ecl_sum);
 }
 
 
 int main() {
-  fprintf(stderr,"The ecl_sum resample code is currently broken - this should be fizxed\n");
-  exit(0);
-  // test_correct_time_vector();
-  // test_not_sorted();
-  // return 0;
+  test_correct_time_vector();
+  test_resample_extrapolate_rate();
+  test_not_sorted();
+  return 0;
 }
 

lib/include/ert/ecl/ecl_sum.hpp

@@ -93,7 +93,7 @@ typedef struct ecl_sum_struct       ecl_sum_type;
   ecl_sum_type   * ecl_sum_fread_alloc_case(const char *  , const char * key_join_string);
   ecl_sum_type   * ecl_sum_fread_alloc_case__(const char * input_file , const char * key_join_string , bool include_restart);
   ecl_sum_type   * ecl_sum_fread_alloc_case2__(const char *  , const char * key_join_string , bool include_restart, bool lazy_load, int file_options);
-  ecl_sum_type   * ecl_sum_alloc_resample(const ecl_sum_type * ecl_sum, const char * ecl_case, const time_t_vector_type * times);
+  ecl_sum_type   * ecl_sum_alloc_resample(const ecl_sum_type * ecl_sum, const char * ecl_case, const time_t_vector_type * times, bool lower_extrapolation, bool upper_extrapolation);
   bool             ecl_sum_case_exists( const char * input_file );
 
   /* Accessor functions : */

python/ecl/summary/ecl_sum.py

@@ -92,7 +92,7 @@ class EclSum(BaseCClass):
     _fread_alloc                   = EclPrototype("void*     ecl_sum_fread_alloc(char*, stringlist, char*, bool)", bind=False)
     _create_restart_writer         = EclPrototype("ecl_sum_obj  ecl_sum_alloc_restart_writer2(char*, char*, int, bool, bool, char*, time_t, bool, int, int, int)", bind = False)
     _create_writer                 = EclPrototype("ecl_sum_obj  ecl_sum_alloc_writer(char*, bool, bool, char*, time_t, bool, int, int, int)", bind = False)
-    _resample                      = EclPrototype("ecl_sum_obj  ecl_sum_alloc_resample( ecl_sum, char*, time_t_vector)")
+    _resample                      = EclPrototype("ecl_sum_obj  ecl_sum_alloc_resample( ecl_sum, char*, time_t_vector, bool, bool)")
     _iiget                         = EclPrototype("double   ecl_sum_iget(ecl_sum, int, int)")
     _free                          = EclPrototype("void     ecl_sum_free(ecl_sum)")
     _data_length                   = EclPrototype("int      ecl_sum_get_data_length(ecl_sum)")
@@ -1537,8 +1537,8 @@ def export_csv(self, filename, keys=None, date_format="%Y-%m-%d", sep=";"):
 
 
 
-    def resample(self, new_case_name, time_points):
-        new_case = self._resample(new_case_name, time_points)
+    def resample(self, new_case_name, time_points, lower_extrapolation=False, upper_extrapolation=False):
+        new_case = self._resample(new_case_name, time_points, lower_extrapolation, upper_extrapolation)
         if new_case is None:
             raise ValueError("Failed to create new resampled case:{}".format(new_case_name))
 

python/tests/ecl_tests/test_sum.py

@@ -620,3 +620,49 @@ def test_directory_conflict(self):
             case.fwrite()
             os.mkdir("UNITS")
             case2 = EclSum("./UNITS")
+
+
+    def test_resample_extrapolate(self):
+        """
+        Test resampling of summary with extrapolate option of lower and upper boundaries enabled
+        Note: When performing resampling on cp_simple3 test case, it fails to duplicate node 251 so using mocked ecl_sum instead
+        path = os.path.join(self.TESTDATA_ROOT, "local/ECLIPSE/cp_simple3/SIMPLE_SUMMARY3")
+        ecl_sum = EclSum( path, lazy_load=True )
+        """
+        from ecl.util.util import TimeVector, CTime
+
+        time_points = TimeVector()
+        ecl_sum = create_case(data_start=datetime.date(2010, 1, 1))
+        start_time = ecl_sum.get_data_start_time() - datetime.timedelta(seconds=86400)
+        end_time = ecl_sum.get_end_time() + datetime.timedelta(seconds=86400)
+        delta = end_time - start_time
+
+        N = 25
+        time_points.initRange( CTime(start_time),
+                               CTime(end_time),
+                               CTime(int(delta.total_seconds()/(N - 1))))
+        time_points.append(CTime(end_time))
+        resampled = ecl_sum.resample( "OUTPUT_CASE", time_points, lower_extrapolation=True, upper_extrapolation=True )
+
+        for key in ecl_sum.keys():
+            self.assertIn( key, resampled )
+
+        self.assertEqual(ecl_sum.get_data_start_time() -  datetime.timedelta(seconds=86400), resampled.get_data_start_time())
+
+        key_not_rate = "FOPT"
+        for time_index,t in enumerate(time_points):
+            if t < ecl_sum.get_data_start_time():
+                self.assertFloatEqual(resampled.iget( key_not_rate, time_index), ecl_sum._get_first_value(key_not_rate))
+            elif t >  ecl_sum.get_end_time():
+                self.assertFloatEqual(resampled.iget( key_not_rate, time_index), ecl_sum.get_last_value( key_not_rate))
+            else:
+                self.assertFloatEqual(resampled.iget( key_not_rate, time_index), ecl_sum.get_interp_direct( key_not_rate, t))
+
+        key_rate = "FOPR"
+        for time_index,t in enumerate(time_points):
+            if t < ecl_sum.get_data_start_time():
+                self.assertFloatEqual(resampled.iget( key_rate, time_index), 0)
+            elif t >  ecl_sum.get_end_time():
+                self.assertFloatEqual(resampled.iget( key_rate, time_index), 0)
+            else:
+                self.assertFloatEqual(resampled.iget( key_rate, time_index), ecl_sum.get_interp_direct( key_rate, t))

python/tests/ecl_tests/test_sum_statoil.py

@@ -525,7 +525,6 @@ def test_resample(self):
             for time_index,t in enumerate(time_points):
                 self.assertFloatEqual(resampled.iget( key, time_index), self.ecl_sum.get_interp_direct( key, t))
 
-
     def test_summary_units(self):
         self.assertEqual(self.ecl_sum.unit_system, EclUnitTypeEnum.ECL_METRIC_UNITS)