Fix error in triangular distribution

equinor · Jan 16, 2020 · 3585997 · 3585997
1 parent ddfd865
commit 3585997
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 3 deletions.
diff --git a/lib/enkf/tests/trans_func.cpp b/lib/enkf/tests/trans_func.cpp
@@ -40,6 +40,20 @@ void test_triangular() {
   stringlist_free(args);
 }
 
+void test_triangular_assymetric() {
+  stringlist_type * args = stringlist_alloc_new();
+  stringlist_append_copy(args , "TRIANGULAR");
+  stringlist_append_copy(args, "0");
+  stringlist_append_copy(args,"1.0");
+  stringlist_append_copy(args, "4.0");
+
+  trans_func_type * trans_func = trans_func_alloc(args);
+  test_assert_double_equal( trans_func_eval(trans_func, -1.0), 0.7966310411513150456286);
+  test_assert_double_equal( trans_func_eval(trans_func, 1.1), 2.72407181575270778882286);
+  trans_func_free( trans_func );
+  stringlist_free(args);
+}
+
 void test_create() {
   {
     stringlist_type * args = stringlist_alloc_new();
@@ -79,5 +93,6 @@ void test_create() {
 int main(int argc , char ** argv) {
   test_create();
   test_triangular();
+  test_triangular_assymetric();
 }
 
diff --git a/lib/enkf/trans_func.cpp b/lib/enkf/trans_func.cpp
@@ -182,14 +182,15 @@ static double trans_triangular(double x, const double_vector_type * arg) {
   double xmode = double_vector_iget(arg, 1);
   double xmax = double_vector_iget(arg,2);
 
-  double inv_norm = (xmax - xmin) * (xmode- xmin);
+  double inv_norm_left = (xmax - xmin) * (xmode - xmin);
+  double inv_norm_right = (xmax - xmin) * (xmax - xmode);
   double ymode = (xmode - xmin) / (xmax - xmin);
   double y = 0.5*(1 + erf(x/sqrt(2.0)));           /* 0 - 1 */
 
   if (y < ymode)
-    return xmin + sqrt(y * inv_norm);
+    return xmin + sqrt(y * inv_norm_left);
   else
-    return xmax - sqrt((1 - y)*inv_norm);
+    return xmax - sqrt((1 - y)*inv_norm_right);
 }