diff --git a/CHANGELOG.md b/CHANGELOG.md index c5e4b9181..873cef3e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,8 @@ Features available from the latest git source [#488](https://github.com/fortran-lang/stdlib/pull/488) - new procedures `arg`, `argd` and `argpi` [#498](https://github.com/fortran-lang/stdlib/pull/498) + - new procedure `diff` + [#605](https://github.com/fortran-lang/stdlib/pull/605) Changes to existing modules diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md index 665d6ee5f..2640d3ed0 100644 --- a/doc/specs/stdlib_math.md +++ b/doc/specs/stdlib_math.md @@ -649,3 +649,88 @@ program demo_math_all_close end program demo_math_all_close ``` + +### `diff` + +#### Description + +Computes differences between adjacent elements of an array. + +#### Syntax + +For a rank-1 array +```fortran +y = [[stdlib_math(module):diff(interface)]](x [, n, prepend, append]) +``` +and for a rank-2 array +```fortran +y = [[stdlib_math(module):diff(interface)]](x [, n, dim, prepend, append]) +``` + +#### Status + +Experimental. + +#### Class + +Pure function. + +#### Arguments + +`x`: The array to take a difference of. +Shall be a `real/integer` and `rank-1/rank-2` array. +This argument is `intent(in)`. + +`n`: How many times to iteratively calculate the difference. +Shall be an `integer` scalar. +This argument is `intent(in)` and `optional`, and has value of `1` by default. + +`dim`: The dimension of the input array along which to calculate the difference. +Its value must be between `1` and `rank(x)`. +Shall be an `integer` scalar. +This argument is `intent(in)` and `optional` and has a value of `1` by default. + +`prepend`, `append`: Arrays to prepend or append to a along axis prior to performing the difference. +The dimension and shape must match a except along axis. +Shall be a `real/integer` and `rank-1/rank-2` array. +This argument is `intent(in)` and `optional`, which is no value by default. + +Note: +- The `x`, `prepend` and `append` arguments must have the same `type`, `kind` and `rank`. +- If the value of `n` is less than or equal to `0` (which is not recommended), the return value of `diff` is `x`. +- If the value of `dim` is not equal to `1` or `2` (which is not recommended), +`1` will be used by the internal process of `diff`. + + +#### Result value + +Returns the finite difference of the input array. +Shall be a `real/integer` and `rank-1/rank-2` array. +When both `prepend` and `append` are not present, the result `y` has one fewer element than `x` alongside the dimension `dim`. + +#### Example + +```fortran +program demo_diff + + use stdlib_math, only: diff + implicit none + + integer :: i(7) = [1, 1, 2, 3, 5, 8, 13] + real :: x(6) = [0, 5, 15, 30, 50, 75] + integer :: A(3, 3) = reshape([1, 7, 17, 3, 11, 19, 5, 13, 23], [3, 3]) + integer :: Y(3, 2) + + print *, diff(i) ! [0, 1, 1, 2, 3, 5] + print *, diff(x, 2) ! [5.0, 5.0, 5.0, 5.0] + + Y = diff(A, n=1, dim=2) + print *, Y(1, :) ! [2, 2] + print *, Y(2, :) ! [4, 2] + print *, Y(3, :) ! [2, 4] + + print *, diff(i, prepend=[0]) ! [1, 0, 1, 1, 2, 3, 5] + print *, diff(i, append=[21]) ! [0, 1, 1, 2, 3, 5, 8] + +end program demo_diff +``` \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d3a107e54..ed615311e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -42,6 +42,7 @@ set(fppFiles stdlib_math_arange.fypp stdlib_math_is_close.fypp stdlib_math_all_close.fypp + stdlib_math_diff.fypp stdlib_string_type.fypp stdlib_string_type_constructor.fypp stdlib_strings_to_string.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index f8e001377..d61e75cb8 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -39,6 +39,7 @@ SRCFYPP = \ stdlib_math_logspace.fypp \ stdlib_math_is_close.fypp \ stdlib_math_all_close.fypp \ + stdlib_math_diff.fypp \ stdlib_string_type.fypp \ stdlib_string_type_constructor.fypp \ stdlib_strings.fypp \ @@ -205,6 +206,7 @@ stdlib_math_is_close.o: \ stdlib_math_all_close.o: \ stdlib_math.o \ stdlib_math_is_close.o +stdlib_math_diff.o: stdlib_math.o stdlib_stringlist_type.o: stdlib_string_type.o \ stdlib_math.o \ stdlib_optval.o diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp index 82f5961b2..3ffda8afd 100644 --- a/src/stdlib_math.fypp +++ b/src/stdlib_math.fypp @@ -14,7 +14,7 @@ module stdlib_math public :: EULERS_NUMBER_QP #:endif public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH - public :: arange, arg, argd, argpi, is_close, all_close + public :: arange, arg, argd, argpi, is_close, all_close, diff integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100 integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50 @@ -359,6 +359,28 @@ module stdlib_math #:endfor #:endfor end interface all_close + + !> Version: experimental + !> + !> Computes differences between adjacent elements of an array. + !> ([Specification](../page/specs/stdlib_math.html#diff)) + interface diff + #:set RI_KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES + #:for k1, t1 in RI_KINDS_TYPES + pure module function diff_1_${k1}$(x, n, prepend, append) result(y) + ${t1}$, intent(in) :: x(:) + integer, intent(in), optional :: n + ${t1}$, intent(in), optional :: prepend(:), append(:) + ${t1}$, allocatable :: y(:) + end function diff_1_${k1}$ + pure module function diff_2_${k1}$(X, n, dim, prepend, append) result(y) + ${t1}$, intent(in) :: x(:, :) + integer, intent(in), optional :: n, dim + ${t1}$, intent(in), optional :: prepend(:, :), append(:, :) + ${t1}$, allocatable :: y(:, :) + end function diff_2_${k1}$ + #:endfor + end interface diff contains diff --git a/src/stdlib_math_diff.fypp b/src/stdlib_math_diff.fypp new file mode 100644 index 000000000..eb8cb0bc2 --- /dev/null +++ b/src/stdlib_math_diff.fypp @@ -0,0 +1,139 @@ +!> Inspired by original code (MIT license) written in 2016 by Keurfon Luu (keurfonluu@outlook.com) +!> https://github.com/keurfonluu/Forlab + +#:include "common.fypp" +#:set RI_KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES +submodule (stdlib_math) stdlib_math_diff + + implicit none + +contains + + !> `diff` computes differences of adjacent elements of an array. + + #:for k1, t1 in RI_KINDS_TYPES + pure module function diff_1_${k1}$(x, n, prepend, append) result(y) + ${t1}$, intent(in) :: x(:) + integer, intent(in), optional :: n + ${t1}$, intent(in), optional :: prepend(:), append(:) + ${t1}$, allocatable :: y(:) + integer :: size_prepend, size_append, size_x, size_work + integer :: n_, i + + n_ = optval(n, 1) + if (n_ <= 0) then + y = x + return + end if + + size_prepend = 0 + size_append = 0 + if (present(prepend)) size_prepend = size(prepend) + if (present(append)) size_append = size(append) + size_x = size(x) + size_work = size_x + size_prepend + size_append + + if (size_work <= n_) then + allocate(y(0)) + return + end if + + !> Use a quick exit for the common case, to avoid memory allocation. + if (size_prepend == 0 .and. size_append == 0 .and. n_ == 1) then + y = x(2:) - x(1:size_x-1) + return + end if + + block + ${t1}$ :: work(size_work) + if (size_prepend > 0) work(:size_prepend) = prepend + work(size_prepend+1:size_prepend+size_x) = x + if (size_append > 0) work(size_prepend+size_x+1:) = append + + do i = 1, n_ + work(1:size_work-i) = work(2:size_work-i+1) - work(1:size_work-i) + end do + + y = work(1:size_work-n_) + end block + + end function diff_1_${k1}$ + + pure module function diff_2_${k1}$(x, n, dim, prepend, append) result(y) + ${t1}$, intent(in) :: x(:, :) + integer, intent(in), optional :: n, dim + ${t1}$, intent(in), optional :: prepend(:, :), append(:, :) + ${t1}$, allocatable :: y(:, :) + integer :: size_prepend, size_append, size_x, size_work + integer :: n_, dim_, i + + n_ = optval(n, 1) + if (n_ <= 0) then + y = x + return + end if + + size_prepend = 0 + size_append = 0 + if (present(dim)) then + if (dim == 1 .or. dim == 2) then + dim_ = dim + else + dim_ = 1 + end if + else + dim_ = 1 + end if + + if (present(prepend)) size_prepend = size(prepend, dim_) + if (present(append)) size_append = size(append, dim_) + size_x = size(x, dim_) + size_work = size_x + size_prepend + size_append + + if (size_work <= n_) then + allocate(y(0, 0)) + return + end if + + !> Use a quick exit for the common case, to avoid memory allocation. + if (size_prepend == 0 .and. size_append == 0 .and. n_ == 1) then + if (dim_ == 1) then + y = x(2:, :) - x(1:size_x-1, :) + elseif (dim_ == 2) then + y = x(:, 2:) - x(:, 1:size_x-1) + end if + return + end if + + if (dim_ == 1) then + block + ${t1}$ :: work(size_work, size(x, 2)) + if (size_prepend > 0) work(1:size_prepend, :) = prepend + work(size_prepend+1:size_x+size_prepend, :) = x + if (size_append > 0) work(size_x+size_prepend+1:, :) = append + do i = 1, n_ + work(1:size_work-i, :) = work(2:size_work-i+1, :) - work(1:size_work-i, :) + end do + + y = work(1:size_work-n_, :) + end block + + elseif (dim_ == 2) then + block + ${t1}$ :: work(size(x, 1), size_work) + if (size_prepend > 0) work(:, 1:size_prepend) = prepend + work(:, size_prepend+1:size_x+size_prepend) = x + if (size_append > 0) work(:, size_x+size_prepend+1:) = append + do i = 1, n_ + work(:, 1:size_work-i) = work(:, 2:size_work-i+1) - work(:, 1:size_work-i) + end do + + y = work(:, 1:size_work-n_) + end block + + end if + + end function diff_2_${k1}$ + #:endfor + +end submodule stdlib_math_diff \ No newline at end of file diff --git a/src/tests/math/test_stdlib_math.fypp b/src/tests/math/test_stdlib_math.fypp index 4d8bb34a5..bf061057b 100644 --- a/src/tests/math/test_stdlib_math.fypp +++ b/src/tests/math/test_stdlib_math.fypp @@ -4,7 +4,7 @@ module test_stdlib_math use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test - use stdlib_math, only: clip, arg, argd, argpi, arange, is_close, all_close + use stdlib_math, only: clip, arg, argd, argpi, arange, is_close, all_close, diff use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp implicit none @@ -51,6 +51,14 @@ contains , new_unittest("all_close-real-${k1}$", test_all_close_real_${k1}$) & , new_unittest("all_close-cmplx-${k1}$", test_all_close_cmplx_${k1}$) & #:endfor + + !> Tests for `diff` + #:for k1 in REAL_KINDS + , new_unittest("diff-real-${k1}$", test_diff_real_${k1}$) & + #:endfor + #:for k1 in INT_KINDS + , new_unittest("diff-int-${k1}$", test_diff_int_${k1}$) & + #:endfor ] end subroutine collect_stdlib_math @@ -364,6 +372,99 @@ contains end subroutine test_all_close_cmplx_${k1}$ #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + subroutine test_diff_real_${k1}$(error) + type(error_type), allocatable, intent(out) :: error + ${t1}$ :: x(6) = [${t1}$ :: 0, 5, 15, 30, 50, 75] + ${t1}$ :: A(1, 3) = reshape([${t1}$ :: 1, 3, 5], [1, 3]) + ${t1}$ :: B(2) = [${t1}$ :: 1, 2] + + !> rank-1 diff + call check(error, all_close(diff(x), [${t1}$ :: 5, 10, 15, 20, 25]), & + "diff() in test_diff_real_${k1}$ failed") + if (allocated(error)) return + call check(error, all_close(diff(x, n=0), x), & + "diff(, n=0) in test_diff_real_${k1}$ failed") + call check(error, all_close(diff(x, n=2), [${t1}$ :: 5, 5, 5, 5]), & + "diff(, n=2) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + + call check(error, all_close(diff(x, prepend=[${t1}$ :: 1]), [${t1}$ :: -1, 5, 10, 15, 20, 25]), & + "diff(, prepend=[${t1}$ :: 1]) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + call check(error, all_close(diff(x, append=[${t1}$ :: 1]), [${t1}$ :: 5, 10, 15, 20, 25, -74]), & + "diff(, append=[${t1}$ :: 1]) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + + !> rank-2 diff + call check(error, all_close(diff(reshape(A, [3,1]), n=1, dim=1), reshape([${t1}$ :: 2, 2], [2, 1])), & + "diff(, n=1, dim=1) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + call check(error, all_close(diff(A, n=1, dim=2), reshape([${t1}$ :: 2, 2], [1, 2])), & + "diff(, n=1, dim=2) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + + call check(error, all_close(diff(A, n=1, dim=2, prepend=reshape([${t1}$ :: 1], [1, 1]), & + append=reshape([${t1}$ :: 2], [1, 1])), reshape([${t1}$ :: 0, 2, 2, -3], [1, 4])), & + "diff(, n=1, dim=2, prepend=reshape([${t1}$ :: 1], [1, 1]), & + &append=reshape([${t1}$ :: 2], [1, 1])) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + + !> size(B, dim) <= n + call check(error, size(diff(B, 2)), 0, "size(diff(B, 2)) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + call check(error, size(diff(B, 3)), 0, "size(diff(B, 3)) in test_diff_real_${k1}$ failed") + + end subroutine test_diff_real_${k1}$ + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + subroutine test_diff_int_${k1}$(error) + type(error_type), allocatable, intent(out) :: error + ${t1}$ :: x(6) = [${t1}$ :: 0, 5, 15, 30, 50, 75] + ${t1}$ :: A(1, 3) = reshape([${t1}$ :: 1, 3, 5], [1, 3]) + ${t1}$ :: B(2) = [${t1}$ :: 1, 2] + + !> rank-1 diff + call check(error, all(diff(x) == [${t1}$ :: 5, 10, 15, 20, 25]), & + "diff() in test_diff_int_${k1}$ failed") + if (allocated(error)) return + call check(error, all(diff(x, n=0) == x), & + "diff(, n=0) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + call check(error, all(diff(x, n=2) == [${t1}$ :: 5, 5, 5, 5]), & + "diff(, n=2) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + + call check(error, all(diff(x, prepend=[${t1}$ :: 1]) == [${t1}$ :: -1, 5, 10, 15, 20, 25]), & + "diff(, prepend=[${t1}$ :: 1]) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + call check(error, all(diff(x, append=[${t1}$ :: 1]) == [${t1}$ :: 5, 10, 15, 20, 25, -74]), & + "diff(, append=[${t1}$ :: 1]) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + + !> rank-2 diff + call check(error, all(diff(reshape(A, [3,1]), n=1, dim=1) == reshape([${t1}$ :: 2, 2], [2, 1])), & + "diff(, n=1, dim=1) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + call check(error, all(diff(A, n=1, dim=2) == reshape([${t1}$ :: 2, 2], [1, 2])), & + "diff(, n=1, dim=2) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + + call check(error, all(diff(A, n=1, dim=2, prepend=reshape([${t1}$ :: 1], [1, 1]), & + append=reshape([${t1}$ :: 2], [1, 1])) == reshape([${t1}$ :: 0, 2, 2, -3], [1, 4])), & + "diff(, n=1, dim=2, prepend=reshape([${t1}$ :: 1], [1, 1]), & + &append=reshape([${t1}$ :: 2], [1, 1])) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + + !> size(B, dim) <= n + call check(error, size(diff(B, 2)), 0, "size(diff(B, 2)) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + call check(error, size(diff(B, 3)), 0, "size(diff(B, 3)) in test_diff_int_${k1}$ failed") + + end subroutine test_diff_int_${k1}$ + #:endfor end module test_stdlib_math