feat: add relative equals that defaults to direct equals for small values

bin/funit/pFUnitParser.py

@@ -21,7 +21,8 @@ def __str__(self):
                    'relativelyequal': 2, 'isinfinite': 1, 'isfinite': 1,
                    'isnan': 1, 'ismemberof': 2, 'contains': 2, 'any': 1,
                    'all': 1, 'notall': 1, 'none': 1, 'ispermutationof': 2,
-                   'exceptionraised': 0, 'sameshape': 2, 'that': 2, '_that': 2}
+                   'exceptionraised': 0, 'sameshape': 2, 'that': 2, '_that': 2,
+                   'approx': 2}
 
 def cppSetLineAndFile(line, file):
     if sysconfig.get_platform() == 'mingw':

src/funit/asserts/Assert.F90

@@ -44,6 +44,7 @@ module PF_Assert
    public :: assertLessThan, assertLessThanOrEqual
    public :: assertGreaterThan, assertGreaterThanOrEqual
    public :: assertRelativelyEqual
+   public :: assertApprox
 
    public :: assertIsNan, assertIsFinite
 

src/funit/asserts/AssertUtilities.F90

@@ -21,6 +21,7 @@ module pf_AssertUtilities
    public :: fail_not_greater_than
    public :: fail_not_greater_than_or_equal
    public :: fail_not_relatively_equal
+   public :: fail_not_approx
 
 contains
 
@@ -342,7 +343,33 @@ subroutine fail_not_relatively_equal(expected, actual, difference, unused, index
       call throw(fail_message, location)
 
    end subroutine fail_not_relatively_equal
-
+
+   subroutine fail_not_approx(expected, actual, difference, unused, index, message, location)
+      character(*), intent(in) :: expected
+      character(*), intent(in) :: actual
+      character(*), intent(in) :: difference
+      ! Separator
+      class (KeywordEnforcer), optional, intent(in) :: unused
+      ! Keyword arguments
+      integer, optional, intent(in) :: index(:)
+      character(*), optional, intent(in) :: message
+      type (SourceLocation), optional, intent(in) :: location
+
+      character(len=:), allocatable :: fail_message
+
+      _UNUSED_DUMMY(unused)
+
+      fail_message = base_message('AssertApprox', message, index)
+      fail_message = fail_message // new_line('A')    // '           Expected: <' // expected // '>' 
+      fail_message = fail_message // new_line('A')    // '             Actual: <' // actual // '>'
+      fail_message = fail_message // new_line('A')    // '         Difference: ' // difference
+      if (present(index)) then
+         fail_message = fail_message // new_line('A') // '      at index: ' // toString(index)
+      end if
+
+      call throw(fail_message, location)
+
+   end subroutine fail_not_approx
 
    function base_message(failure_type, user_message, index) result(message)
       character(:), allocatable :: message

src/funit/asserts/Assert_Real.tmpl

@@ -122,6 +122,79 @@
 [(real,  128, 0), (real,  128, rank), (real, 128, 0)]
 @end tkr_parameters
 
+@tkr_parameters with_rel_tol
+! Tolerance is default real
+[(integer,  default, rank), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(real,  default, rank), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(integer,  default, rank), (real,  double, rank), (real, default, 0), (real, default, 0)]
+[(real,  default, rank), (real,  double, rank), (real, default, 0), (real, default, 0)]
+[(real,  double, rank), (real,  double, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, rank), (real,  32, rank), (real, default, 0), (real, default, 0)]
+[(real, default, rank), (real,  32, rank), (real, default, 0), (real, default, 0)]
+[(real,  32, rank), (real,  32, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, rank), (real,  64, rank), (real, default, 0), (real, default, 0)]
+[(real, default, rank), (real,  64, rank), (real, default, 0), (real, default, 0)]
+[(real,  64, rank), (real,  64, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, rank), (real,  80, rank), (real, default, 0), (real, default, 0)]
+[(real, default, rank), (real,  80, rank), (real, default, 0), (real, default, 0)]
+[(real,  80, rank), (real,  80, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, rank), (real,  128, rank), (real, default, 0), (real, default, 0)]
+[(real, default, rank), (real,  128, rank), (real, default, 0), (real, default, 0)]
+[(real,  128, rank), (real,  128, rank), (real, default, 0), (real, default, 0)]
+! Tolerance is same kind as actual
+[(integer,  default, rank), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(real,  default, rank), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(integer,  default, rank), (real,  double, rank), (real, double, 0), (real, double, 0)]
+[(real,  default, rank), (real,  double, rank), (real, double, 0), (real, double, 0)]
+[(real,  double, rank), (real,  double, rank), (real, double, 0), (real, double, 0)]
+[(integer, default, rank), (real,  32, rank), (real, 32, 0), (real, 32, 0)]
+[(real, default, rank), (real,  32, rank), (real, 32, 0), (real, 32, 0)]
+[(real,  32, rank), (real,  32, rank), (real, 32, 0), (real, 32, 0)]
+[(integer, default, rank), (real,  64, rank), (real, 64, 0), (real, 64, 0)]
+[(real, default, rank), (real,  64, rank), (real, 64, 0), (real, 64, 0)]
+[(real,  64, rank), (real,  64, rank), (real, 64, 0), (real, 64, 0)]
+[(integer, default, rank), (real,  80, rank), (real, 80, 0), (real, 80, 0)]
+[(real, default, rank), (real,  80, rank), (real, 80, 0), (real, 80, 0)]
+[(real,  80, rank), (real,  80, rank), (real, 80, 0), (real, 80, 0)]
+[(integer, default, rank), (real,  128, rank), (real, 128, 0), (real, 128, 0)]
+[(real, default, rank), (real,  128, rank), (real, 128, 0), (real, 128, 0)]
+[(real,  128, rank), (real,  128, rank), (real, 128, 0), (real, 128, 0)]
+! And again with conformable rank 0
+[(integer,  default, 0), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(real,  default, 0), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(integer,  default, 0), (real,  double, rank), (real, default, 0), (real, default, 0)]
+[(real,  default, 0), (real,  double, rank), (real, default, 0), (real, default, 0)]
+[(real,  double, 0), (real,  double, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, 0), (real,  32, rank), (real, default, 0), (real, default, 0)]
+[(real, default, 0), (real,  32, rank), (real, default, 0), (real, default, 0)]
+[(real,  32, 0), (real,  32, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, 0), (real,  64, rank), (real, default, 0), (real, default, 0)]
+[(real, default, 0), (real,  64, rank), (real, default, 0), (real, default, 0)]
+[(real,  64, 0), (real,  64, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, 0), (real,  80, rank), (real, default, 0), (real, default, 0)]
+[(real, default, 0), (real,  80, rank), (real, default, 0), (real, default, 0)]
+[(real,  80, 0), (real,  80, rank), (real, default, 0), (real, default, 0)]
+[(integer, default, 0), (real,  128, rank), (real, default, 0), (real, default, 0)]
+[(real, default, 0), (real,  128, rank), (real, default, 0), (real, default, 0)]
+[(real,  128, 0), (real,  128, rank), (real, default, 0), (real, default, 0)]
+[(integer,  default, 0), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(real,  default, 0), (real,  default, rank), (real, default, 0), (real, default, 0)]
+[(integer,  default, 0), (real,  double, rank), (real, double, 0), (real, double, 0)]
+[(real,  default, 0), (real,  double, rank), (real, double, 0), (real, double, 0)]
+[(real,  double, 0), (real,  double, rank), (real, double, 0), (real, double, 0)]
+[(integer, default, 0), (real,  32, rank), (real, 32, 0), (real, 32, 0)]
+[(real, default, 0), (real,  32, rank), (real, 32, 0), (real, 32, 0)]
+[(real,  32, 0), (real,  32, rank), (real, 32, 0), (real, 32, 0)]
+[(integer, default, 0), (real,  64, rank), (real, 64, 0), (real, 64, 0)]
+[(real, default, 0), (real,  64, rank), (real, 64, 0), (real, 64, 0)]
+[(real,  64, 0), (real,  64, rank), (real, 64, 0), (real, 64, 0)]
+[(integer, default, 0), (real,  80, rank), (real, 80, 0), (real, 80, 0)]
+[(real, default, 0), (real,  80, rank), (real, 80, 0), (real, 80, 0)]
+[(real,  80, 0), (real,  80, rank), (real, 80, 0), (real, 80, 0)]
+[(integer, default, 0), (real,  128, rank), (real, 128, 0), (real, 128, 0)]
+[(real, default, 0), (real,  128, rank), (real, 128, 0), (real, 128, 0)]
+[(real,  128, 0), (real,  128, rank), (real, 128, 0), (real, 128, 0)]
+@end tkr_parameters
 
 #include "unused_dummy.fh"
 
@@ -176,6 +249,7 @@ module pf_AssertReal_{rank}d
    @overload(AssertGreaterThanOrEqual, no_tol)
 
    @overload(AssertRelativelyEqual, with_tol)
+   @overload(AssertApprox, with_rel_tol)
    @overload(AssertAssociated, minimal)
 
    @overload(assert_equal, minimal)
@@ -185,6 +259,7 @@ module pf_AssertReal_{rank}d
    @overload(assert_greater_than, minimal)
    @overload(assert_greater_than_or_equal, minimal)
    @overload(assert_relatively_equal, minimal)
+   @overload(assert_approx, minimal)
 
 
    integer, parameter :: MAX_LEN_REAL_AS_STRING = 40
@@ -981,7 +1056,112 @@ contains
    end subroutine {name}
    @end template
 
+   @template(AssertApprox,[expected,actual,rel_tolerance,tolerance])
+   subroutine {name}(expected, actual, rel_tolerance, tolerance, message, location)
+      {expected.type} (kind={expected.kind}), intent(in) :: expected {expected.dims}
+      {actual.type} (kind={actual.kind}), intent(in) :: actual {actual.dims}
+      {rel_tolerance.type} (kind={rel_tolerance.kind}), intent(in) :: rel_tolerance
+      {tolerance.type} (kind={tolerance.kind}), intent(in) :: tolerance
+      character(*), optional, intent(in) :: message
+      type (SourceLocation), optional, intent(in) :: location
+
+      real(kind=kind(actual)) :: rel_t
+      real(kind=kind(actual)) :: t
+      real(kind=kind(actual)), allocatable :: e {actual.dims}
+
+      if (.not. conformable(shape(expected), shape(actual))) then
+         call fail_not_conformable(shape(expected), shape(actual), message=message, location=location)
+         return
+      end if
+
+      ! Trick to get e to have the right shape even if "expected" is a scalar.
+      e = 0*actual + expected
+
+      !if (any([e] == 0)) then
+      !   call fail_generic('Small denominator detected in AssertApprox.',message=message, location=location)
+      !   return
+      !end if
+
+      rel_t = real(rel_tolerance, kind(actual))
+      t = real(tolerance, kind(actual))
+
+      call assert_approx(e, actual, rel_t, t, message=message, location=location)
+
+   end subroutine {name}
+   @end template
 
+   @template(assert_approx,[actual])
+   subroutine {name}(expected, actual, rel_tolerance, tolerance, unused, message, location)
+      {actual.type} (kind={actual.kind}), intent(in) :: expected {actual.dims}
+      {actual.type} (kind={actual.kind}), intent(in) :: actual {actual.dims}
+      {actual.type} (kind={actual.kind}), intent(in) :: rel_tolerance
+      {actual.type} (kind={actual.kind}), intent(in) :: tolerance
+      class (KeywordEnforcer), optional, intent(in) :: unused
+      character(*), optional, intent(in) :: message
+      type (SourceLocation), optional, intent(in) :: location
+
+      real(kind=kind(actual)) :: e, a, d, calc_eps
+
+#if {actual.rank} != 0
+      integer, allocatable :: i(:)
+      logical, allocatable :: within_tolerance {actual.dims}
+#endif
+      character(len=MAX_LEN_REAL_AS_STRING) :: expected_str
+      character(len=MAX_LEN_REAL_AS_STRING) :: actual_str
+      character(len=3*MAX_LEN_REAL_AS_STRING) :: diff_str
+
+      _UNUSED_DUMMY(unused)
+
+#if {actual.rank} == 0
+      ! Scalar
+      ! Handle exact equality separately to avoid subtracting Inf values.
+      ! (Inf - Inf is undefined.)
+      if (actual == expected) return
+
+      if (.not. (abs(actual - expected) <= max(rel_tolerance*abs(expected), tolerance))) then
+         e = expected
+         a = actual
+      else
+         return
+      end if
+
+#else
+      ! Handle exact equality separately to avoid subtracting Inf values.
+      ! (Inf - Inf is undefined.)
+
+      allocate(within_tolerance, mold=(actual==actual))
+      where (actual == expected)
+         within_tolerance = .true.
+      elsewhere (abs(actual-expected) <= max(rel_tolerance*abs(expected), tolerance))
+         within_tolerance = .true.
+      elsewhere
+         within_tolerance = .false.
+      end where
+
+      if (all(within_tolerance)) return
+
+      ! index of first difference is
+      i = maxloc(merge(1,0, .not. abs(actual-expected) <= max(rel_tolerance*abs(expected), tolerance)))
+      e = expected({actual.multi_index})
+      a = actual({actual.multi_index})
+#endif
+      calc_eps = max(rel_tolerance*abs(e), tolerance)
+      d = (a - e)
+
+      ! Wish: allocatable strings were useful as internal files ...
+      write(expected_str,'(g0)') e
+      write(actual_str,'(g0)') a
+      write(diff_str,'("<",g0,"> (greater than calculated tolerance of ",g0,")")') d, calc_eps
+
+#if {actual.rank} == 0
+      call fail_not_approx(trim(expected_str), trim(actual_str), trim(diff_str), &
+           & message=message, location=location)
+#else
+      call fail_not_approx(trim(expected_str), trim(actual_str), trim(diff_str), index=i, &
+           & message=message, location=location)
+#endif
+   end subroutine {name}
+   @end template 
 
    @instantiate(AssertEqual, with_tol)
    @instantiate(AssertNotEqual, with_tol)
@@ -990,6 +1170,7 @@ contains
    @instantiate(AssertGreaterThan, with_tol)
    @instantiate(AssertGreaterThanOrEqual, with_tol)
    @instantiate(AssertRelativelyEqual, with_tol)
+   @instantiate(AssertApprox, with_rel_tol)
 
    @instantiate(AssertEqual, no_tol)
    @instantiate(AssertNotEqual, no_tol)
@@ -1007,6 +1188,7 @@ contains
    @instantiate(assert_greater_than, minimal)
    @instantiate(assert_greater_than_or_equal, minimal)
    @instantiate(assert_relatively_equal, minimal)
+   @instantiate(assert_approx, minimal)
 
 end module pf_AssertReal_{rank}d
 

tests/funit-core/CMakeLists.txt

@@ -28,6 +28,7 @@ set(pf_tests
   Test_AssertEqual_Real.pf
   Test_AssertNotEqual_Real.pf
   Test_AssertRelativelyEqual_Real.pf
+  Test_AssertApprox_Real.pf
   Test_AssertLessThan_Real.pf
   Test_AssertLessThanOrEqual_Real.pf
   Test_AssertGreaterThan_Real.pf