From 2827e9aaedcd487e84c3dd2e4dc2dae7ec92638e Mon Sep 17 00:00:00 2001 From: zoziha Date: Sun, 19 Dec 2021 23:03:10 +0800 Subject: [PATCH 1/4] Add procedure `diff`. --- CHANGELOG.md | 2 + doc/specs/stdlib_math.md | 69 ++++++++++++++++++++++++++++ src/CMakeLists.txt | 1 + src/Makefile.manual | 2 + src/stdlib_math.fypp | 21 +++++++++ src/stdlib_math_diff.fypp | 61 ++++++++++++++++++++++++ src/tests/math/test_stdlib_math.fypp | 52 ++++++++++++++++++++- 7 files changed, 207 insertions(+), 1 deletion(-) create mode 100644 src/stdlib_math_diff.fypp diff --git a/CHANGELOG.md b/CHANGELOG.md index 319e80071..119ba1fea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,8 @@ Features available from the latest git source - new module `stdlib_io_npy` [#581](https://github.com/fortran-lang/stdlib/pull/581) - new procedures `save_npy`, `load_npy` +- update module `stdlib_math` + - new procedure `diff` Changes to existing modules diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md index 2ce8bfce6..e621d26c1 100644 --- a/doc/specs/stdlib_math.md +++ b/doc/specs/stdlib_math.md @@ -528,3 +528,72 @@ program demo_math_all_close end program demo_math_all_close ``` + +### `diff` + +#### Description + +Computes differences between adjacent elements of the input array. + +#### Syntax + +```fortran +!> For rank-1 array +Y = [[stdlib_math(module):diff(interface)]](X [, n]) +!> and for rank-2 array +Y = [[stdlib_math(module):diff(interface)]](X [, n, dim]) +``` + +#### Status + +Experimental. + +#### Class + +Pure function. + +#### Arguments + +Note: If the value of `dim` is not equal to `1` or `2`, +`1` will be used by the internal process of `diff`. (Not recommended) + +`X`: Shall be a `real/integer` and `rank-1/rank-2` array. +This argument is `intent(in)`. + +`n`: Shall be an `integer` scalar. +This argument is `intent(in)` and `optional`, which is `1` by default. +It represents to calculate the n-th order difference. + +`dim`: Shall be an `integer` scalar. +This argument is `intent(in)` and `optional`, which is `1` by default. +It represents to calculate the difference along which dimension. + +#### Result value + +Note: That `Y` has one fewer element than `X`. + +Returns a `real/integer` and `rank-1/rank-2` array. + +#### 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] + +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 2db5d2a56..aa3519e60 100644 --- a/src/stdlib_math.fypp +++ b/src/stdlib_math.fypp @@ -15,6 +15,7 @@ module stdlib_math #:endif public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH public :: arange, is_close, all_close + public :: diff integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100 integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50 @@ -326,6 +327,26 @@ module stdlib_math #:endfor #:endfor end interface all_close + + !> Version: experimental + !> + !> Computes differences between adjacent elements of the input 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) result(Y) + ${t1}$, intent(in) :: x(:) + integer, intent(in), optional :: n + ${t1}$, allocatable :: Y(:) + end function diff_1_${k1}$ + pure module function diff_2_${k1}$(X, n, dim) result(Y) + ${t1}$, intent(in) :: X(:, :) + integer, intent(in), optional :: n, dim + ${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..f330ec454 --- /dev/null +++ b/src/stdlib_math_diff.fypp @@ -0,0 +1,61 @@ +#! Originally written in 2016 by Keurfon Luu (keurfonluu@outlook.com) +#! https://github.com/keurfonluu/Forlab/blob/master/src/lib/forlab.f90#L2673 + +#: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 arrays of the ${t1}$ type. + + #:for k1, t1 in RI_KINDS_TYPES + pure module function diff_1_${k1}$(X, n) result(Y) + ${t1}$, intent(in) :: X(:) + integer, intent(in), optional :: n + ${t1}$, allocatable :: Y(:) + integer :: n_, i + + n_ = optval(n, 1) + + Y = X + do i = 1, n_ + Y = Y(2:) - Y(:size(Y) - 1) + end do + + end function diff_1_${k1}$ + + pure module function diff_2_${k1}$(X, n, dim) result(Y) + ${t1}$, intent(in) :: X(:, :) + integer, intent(in), optional :: n, dim + ${t1}$, allocatable :: Y(:, :) + integer :: n_, dim_, i + + n_ = optval(n, 1) + if (present(dim)) then + if (dim == 1 .or. dim == 2) then + dim_ = dim + else + dim_ = 1 + end if + else + dim_ = 1 + end if + + Y = X + if (dim_ == 1) then + do i = 1, n_ + Y = Y(2:, :) - Y(:size(Y, 1) - 1, :) + end do + elseif (dim == 2) then + do i = 1, n_ + Y = Y(:, 2:) - Y(:, :size(Y, 2) - 1) + end do + 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 72c3adf5b..1bf24dc68 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, is_close, all_close + use stdlib_math, only: clip, is_close, all_close, diff use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp implicit none @@ -40,6 +40,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 @@ -294,6 +302,48 @@ 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]) + + call check(error, all_close(diff(x), [${t1}$ :: 5, 10, 15, 20, 25]), & + "diff(x) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + + call check(error, all_close(diff(x, n=2), [${t1}$ :: 5, 5, 5, 5]), & + "diff(x, n=2) 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(x, n=1, dim=2) in test_diff_real_${k1}$ failed") + if (allocated(error)) return + + 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]) + + call check(error, all(diff(x) == [${t1}$ :: 5, 10, 15, 20, 25]), & + "diff(x) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + + call check(error, all(diff(x, n=2) == [${t1}$ :: 5, 5, 5, 5]), & + "diff(x, n=2) 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(A, n=1, dim=2) in test_diff_int_${k1}$ failed") + if (allocated(error)) return + + end subroutine test_diff_int_${k1}$ + #:endfor end module test_stdlib_math From ec00dfeb3ecb19e0aee56071630e001e9df9e2e7 Mon Sep 17 00:00:00 2001 From: zoziha Date: Sun, 2 Jan 2022 09:35:19 +0800 Subject: [PATCH 2/4] Update implementaion of `diff` procedure. --- CHANGELOG.md | 1 + doc/specs/stdlib_math.md | 30 +++++++--- src/stdlib_math.fypp | 14 +++-- src/stdlib_math_diff.fypp | 84 +++++++++++++++++++++++----- src/tests/math/test_stdlib_math.fypp | 59 ++++++++++++++++--- 5 files changed, 152 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9dd28ac91..873cef3e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ Features available from the latest git source - 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 c4b231d2e..6e0417d69 100644 --- a/doc/specs/stdlib_math.md +++ b/doc/specs/stdlib_math.md @@ -654,15 +654,17 @@ end program demo_math_all_close #### Description -Computes differences between adjacent elements of the input array. +Computes differences between adjacent elements of an array. #### Syntax +For rank-1 array ```fortran -!> For rank-1 array -Y = [[stdlib_math(module):diff(interface)]](X [, n]) -!> and for rank-2 array -Y = [[stdlib_math(module):diff(interface)]](X [, n, dim]) +y = [[stdlib_math(module):diff(interface)]](x [, n, prepend, append]) +``` +and for rank-2 array +```fortran +y = [[stdlib_math(module):diff(interface)]](x [, n, dim, prepend, append]) ``` #### Status @@ -675,10 +677,9 @@ Pure function. #### Arguments -Note: If the value of `dim` is not equal to `1` or `2`, -`1` will be used by the internal process of `diff`. (Not recommended) +Note: The `x`, `prepend` and `append` arguments must have same `type`, `kind` and `rank`. -`X`: Shall be a `real/integer` and `rank-1/rank-2` array. +`x`: Shall be a `real/integer` and `rank-1/rank-2` array. This argument is `intent(in)`. `n`: Shall be an `integer` scalar. @@ -689,9 +690,20 @@ It represents to calculate the n-th order difference. This argument is `intent(in)` and `optional`, which is `1` by default. It represents to calculate the difference along which dimension. +`prepend`: Shall be a `real/integer` and `rank-1/rank-2` array, which is no value by default. +This argument is `intent(in)` and `optional`. + +`append`: Shall be a `real/integer` and `rank-1/rank-2` array, which is no value by default. +This argument is `intent(in)` and `optional`. + +Note: +- If the value of `n` is less than or equal to `0`, the return value of `y` is `x`. (Not recommended) +- If the value of `dim` is not equal to `1` or `2`, +`1` will be used by the internal process of `diff`. (Not recommended) + #### Result value -Note: That `Y` has one fewer element than `X`. +Note: That `y` generally has one fewer element than `x`. Returns a `real/integer` and `rank-1/rank-2` array. diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp index 1c4f22ade..3ffda8afd 100644 --- a/src/stdlib_math.fypp +++ b/src/stdlib_math.fypp @@ -362,20 +362,22 @@ module stdlib_math !> Version: experimental !> - !> Computes differences between adjacent elements of the input array. + !> 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) result(Y) + pure module function diff_1_${k1}$(x, n, prepend, append) result(y) ${t1}$, intent(in) :: x(:) integer, intent(in), optional :: n - ${t1}$, allocatable :: Y(:) + ${t1}$, intent(in), optional :: prepend(:), append(:) + ${t1}$, allocatable :: y(:) end function diff_1_${k1}$ - pure module function diff_2_${k1}$(X, n, dim) result(Y) - ${t1}$, intent(in) :: X(:, :) + pure module function diff_2_${k1}$(X, n, dim, prepend, append) result(y) + ${t1}$, intent(in) :: x(:, :) integer, intent(in), optional :: n, dim - ${t1}$, allocatable :: Y(:,:) + ${t1}$, intent(in), optional :: prepend(:, :), append(:, :) + ${t1}$, allocatable :: y(:, :) end function diff_2_${k1}$ #:endfor end interface diff diff --git a/src/stdlib_math_diff.fypp b/src/stdlib_math_diff.fypp index f330ec454..ce2551875 100644 --- a/src/stdlib_math_diff.fypp +++ b/src/stdlib_math_diff.fypp @@ -1,4 +1,4 @@ -#! Originally written in 2016 by Keurfon Luu (keurfonluu@outlook.com) +#! Inspired by original code (MIT license) written in 2016 by Keurfon Luu (keurfonluu@outlook.com) #! https://github.com/keurfonluu/Forlab/blob/master/src/lib/forlab.f90#L2673 #:include "common.fypp" @@ -9,31 +9,65 @@ submodule (stdlib_math) stdlib_math_diff contains - #! `diff` computes differences of arrays of the ${t1}$ type. + #! `diff` computes differences of adjacent elements of an array of the ${t1}$ type. #:for k1, t1 in RI_KINDS_TYPES - pure module function diff_1_${k1}$(X, n) result(Y) - ${t1}$, intent(in) :: X(:) + pure module function diff_1_${k1}$(x, n, prepend, append) result(y) + ${t1}$, intent(in) :: x(:) integer, intent(in), optional :: n - ${t1}$, allocatable :: Y(:) + ${t1}$, intent(in), optional :: prepend(:), append(:) + ${t1}$, allocatable :: y(:) + integer :: size_prepend, size_append, size_x + ${t1}$, allocatable :: 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) + + if (size_x + size_prepend + size_append <= n_) then + allocate(y(0)) + return + end if - Y = X + allocate(work(size_x + size_prepend + size_append)) + 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_ - Y = Y(2:) - Y(:size(Y) - 1) + 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 function diff_1_${k1}$ - pure module function diff_2_${k1}$(X, n, dim) result(Y) - ${t1}$, intent(in) :: X(:, :) + pure module function diff_2_${k1}$(x, n, dim, prepend, append) result(y) + ${t1}$, intent(in) :: x(:, :) integer, intent(in), optional :: n, dim - ${t1}$, allocatable :: Y(:, :) + ${t1}$, intent(in), optional :: prepend(:, :), append(:, :) + ${t1}$, allocatable :: y(:, :) + integer :: size_prepend, size_append, size_x integer :: n_, dim_, i + ${t1}$, allocatable :: work(:, :) 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 @@ -44,15 +78,37 @@ contains dim_ = 1 end if - Y = X + if (present(prepend)) size_prepend = size(prepend, dim_) + if (present(append)) size_append = size(append, dim_) + size_x = size(x, dim_) + + if (size_x + size_prepend + size_append <= n_) then + allocate(y(0, 0)) + return + end if + if (dim_ == 1) then + allocate(work(size_x+size_prepend+size_append, 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_ - Y = Y(2:, :) - Y(:size(Y, 1) - 1, :) + work(1:size(work,1)-i, :) = work(2:size(work)-i+1, :) - work(1:size(x, 1)-i, :) end do - elseif (dim == 2) then + + y = work(1:size(work)-n_, :) + + elseif (dim_ == 2) then + allocate(work(size(x, 1), size_x+size_prepend+size_append)) + 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_ - Y = Y(:, 2:) - Y(:, :size(Y, 2) - 1) + work(:, 1:size(work,2)-i) = work(:, 2:size(work,2)-i+1) - work(:, 1:size(work, 2)-i) end do + + y = work(:, 1:size(work,2)-n_) + end if end function diff_2_${k1}$ diff --git a/src/tests/math/test_stdlib_math.fypp b/src/tests/math/test_stdlib_math.fypp index 6134db56c..2d98de9ca 100644 --- a/src/tests/math/test_stdlib_math.fypp +++ b/src/tests/math/test_stdlib_math.fypp @@ -378,18 +378,43 @@ contains 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(x) in test_diff_real_${k1}$ failed") + "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(x, n=2) in test_diff_real_${k1}$ failed") + "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(x, n=1, dim=2) in test_diff_real_${k1}$ failed") + "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 @@ -399,19 +424,39 @@ contains 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(x) in test_diff_int_${k1}$ failed") + "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(x, n=2) in test_diff_int_${k1}$ failed") + "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(A, n=1, dim=2) 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_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_int_${k1}$ #:endfor From 74c5c111a8099ac269b513a8a7420b2d8194401f Mon Sep 17 00:00:00 2001 From: zoziha Date: Mon, 3 Jan 2022 09:43:08 +0800 Subject: [PATCH 3/4] Improve the efficiency of diff. --- doc/specs/stdlib_math.md | 7 +++- src/stdlib_math_diff.fypp | 60 +++++++++++++++++++--------- src/tests/math/test_stdlib_math.fypp | 18 ++++----- 3 files changed, 55 insertions(+), 30 deletions(-) diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md index 6e0417d69..9baaf6047 100644 --- a/doc/specs/stdlib_math.md +++ b/doc/specs/stdlib_math.md @@ -677,7 +677,7 @@ Pure function. #### Arguments -Note: The `x`, `prepend` and `append` arguments must have same `type`, `kind` and `rank`. +Note: The `x`, `prepend` and `append` arguments must have the same `type`, `kind` and `rank`. `x`: Shall be a `real/integer` and `rank-1/rank-2` array. This argument is `intent(in)`. @@ -688,7 +688,7 @@ It represents to calculate the n-th order difference. `dim`: Shall be an `integer` scalar. This argument is `intent(in)` and `optional`, which is `1` by default. -It represents to calculate the difference along which dimension. +It gives the dimension of the input array along which the difference is calculated, between `1` and `rank(x)`. `prepend`: Shall be a `real/integer` and `rank-1/rank-2` array, which is no value by default. This argument is `intent(in)` and `optional`. @@ -728,5 +728,8 @@ program demo_diff 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/stdlib_math_diff.fypp b/src/stdlib_math_diff.fypp index ce2551875..eb8cb0bc2 100644 --- a/src/stdlib_math_diff.fypp +++ b/src/stdlib_math_diff.fypp @@ -1,5 +1,5 @@ -#! Inspired by original code (MIT license) written in 2016 by Keurfon Luu (keurfonluu@outlook.com) -#! https://github.com/keurfonluu/Forlab/blob/master/src/lib/forlab.f90#L2673 +!> 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 @@ -9,7 +9,7 @@ submodule (stdlib_math) stdlib_math_diff contains - #! `diff` computes differences of adjacent elements of an array of the ${t1}$ type. + !> `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) @@ -17,8 +17,7 @@ contains integer, intent(in), optional :: n ${t1}$, intent(in), optional :: prepend(:), append(:) ${t1}$, allocatable :: y(:) - integer :: size_prepend, size_append, size_x - ${t1}$, allocatable :: work(:) + integer :: size_prepend, size_append, size_x, size_work integer :: n_, i n_ = optval(n, 1) @@ -32,22 +31,31 @@ contains 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_x + size_prepend + size_append <= n_) then + if (size_work <= n_) then allocate(y(0)) return end if - allocate(work(size_x + size_prepend + size_append)) + !> 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) + work(1:size_work-i) = work(2:size_work-i+1) - work(1:size_work-i) end do - - y = work(1:size(work)-n_) + + y = work(1:size_work-n_) + end block end function diff_1_${k1}$ @@ -56,9 +64,8 @@ contains integer, intent(in), optional :: n, dim ${t1}$, intent(in), optional :: prepend(:, :), append(:, :) ${t1}$, allocatable :: y(:, :) - integer :: size_prepend, size_append, size_x + integer :: size_prepend, size_append, size_x, size_work integer :: n_, dim_, i - ${t1}$, allocatable :: work(:, :) n_ = optval(n, 1) if (n_ <= 0) then @@ -81,33 +88,48 @@ contains 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_x + size_prepend + size_append <= n_) then + 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 - allocate(work(size_x+size_prepend+size_append, size(x, 2))) + 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,1)-i, :) = work(2:size(work)-i+1, :) - work(1:size(x, 1)-i, :) + work(1:size_work-i, :) = work(2:size_work-i+1, :) - work(1:size_work-i, :) end do - y = work(1:size(work)-n_, :) + y = work(1:size_work-n_, :) + end block elseif (dim_ == 2) then - allocate(work(size(x, 1), size_x+size_prepend+size_append)) + 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,2)-i) = work(:, 2:size(work,2)-i+1) - work(:, 1:size(work, 2)-i) + work(:, 1:size_work-i) = work(:, 2:size_work-i+1) - work(:, 1:size_work-i) end do - y = work(:, 1:size(work,2)-n_) + y = work(:, 1:size_work-n_) + end block end if diff --git a/src/tests/math/test_stdlib_math.fypp b/src/tests/math/test_stdlib_math.fypp index 2d98de9ca..3ca4c131c 100644 --- a/src/tests/math/test_stdlib_math.fypp +++ b/src/tests/math/test_stdlib_math.fypp @@ -380,7 +380,7 @@ contains ${t1}$ :: A(1, 3) = reshape([${t1}$ :: 1, 3, 5], [1, 3]) ${t1}$ :: B(2) = [${t1}$ :: 1, 2] - #! rank-1 diff + !> 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 @@ -397,7 +397,7 @@ contains "diff(, append=[${t1}$ :: 1]) in test_diff_real_${k1}$ failed") if (allocated(error)) return - #! rank-2 diff + !> 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 @@ -411,7 +411,7 @@ contains &append=reshape([${t1}$ :: 2], [1, 1])) in test_diff_real_${k1}$ failed") if (allocated(error)) return - #! size(B, dim) <= n + !> 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") @@ -426,7 +426,7 @@ contains ${t1}$ :: A(1, 3) = reshape([${t1}$ :: 1, 3, 5], [1, 3]) ${t1}$ :: B(2) = [${t1}$ :: 1, 2] - #! rank-1 diff + !> 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 @@ -444,18 +444,18 @@ contains "diff(, append=[${t1}$ :: 1]) in test_diff_int_${k1}$ failed") if (allocated(error)) return - #! rank-2 diff + !> 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(A, n=1, dim=2) in test_diff_int_${k1}$ failed") + "diff(, n=1, dim=2) 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_real_${k1}$ failed") + !> 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_real_${k1}$ failed") + 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 From cc2344e3f724d4fd7dbe424a0048d8161a68fc37 Mon Sep 17 00:00:00 2001 From: zoziha Date: Sat, 22 Jan 2022 12:37:49 +0800 Subject: [PATCH 4/4] Minor update of diff doc. --- doc/specs/stdlib_math.md | 59 ++++++++++++++-------------- src/tests/math/test_stdlib_math.fypp | 6 +++ 2 files changed, 36 insertions(+), 29 deletions(-) diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md index 9baaf6047..2640d3ed0 100644 --- a/doc/specs/stdlib_math.md +++ b/doc/specs/stdlib_math.md @@ -658,11 +658,11 @@ Computes differences between adjacent elements of an array. #### Syntax -For rank-1 array +For a rank-1 array ```fortran y = [[stdlib_math(module):diff(interface)]](x [, n, prepend, append]) ``` -and for rank-2 array +and for a rank-2 array ```fortran y = [[stdlib_math(module):diff(interface)]](x [, n, dim, prepend, append]) ``` @@ -677,35 +677,36 @@ Pure function. #### Arguments -Note: The `x`, `prepend` and `append` arguments must have the same `type`, `kind` and `rank`. - -`x`: Shall be a `real/integer` and `rank-1/rank-2` array. +`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`: Shall be an `integer` scalar. -This argument is `intent(in)` and `optional`, which is `1` by default. -It represents to calculate the n-th order difference. - -`dim`: Shall be an `integer` scalar. -This argument is `intent(in)` and `optional`, which is `1` by default. -It gives the dimension of the input array along which the difference is calculated, between `1` and `rank(x)`. +`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. -`prepend`: Shall be a `real/integer` and `rank-1/rank-2` array, which is no value by default. -This argument is `intent(in)` and `optional`. +`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. -`append`: Shall be a `real/integer` and `rank-1/rank-2` array, which is no value by default. -This argument is `intent(in)` and `optional`. +`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: -- If the value of `n` is less than or equal to `0`, the return value of `y` is `x`. (Not recommended) -- If the value of `dim` is not equal to `1` or `2`, -`1` will be used by the internal process of `diff`. (Not recommended) +- 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 -Note: That `y` generally has one fewer element than `x`. +#### Result value -Returns a `real/integer` and `rank-1/rank-2` array. +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 @@ -720,16 +721,16 @@ program demo_diff 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] + 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 *, 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] + 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/tests/math/test_stdlib_math.fypp b/src/tests/math/test_stdlib_math.fypp index 3ca4c131c..bf061057b 100644 --- a/src/tests/math/test_stdlib_math.fypp +++ b/src/tests/math/test_stdlib_math.fypp @@ -452,6 +452,12 @@ contains "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