Quadrature: spec tweaks & Simpson's rule #209

nshaffer · 2020-06-03T06:06:11Z

Implementation of `simps` and `simps_weights`

Interfaces and implementations are for 1d real arrays. The included tests exercise the routines for integrands which are (analytically) exactly integrated by Simpson's rule. Known edge cases (length 0, 1, or 6) arrays are explicitly tested as well. There is a lot of new code, and it is not always simple, but hopefully the tests stand as an empirical demonstration of correctness.

Q: To simplify the implementation, simps and simps_weights have been made recursive. This is only an implementation detail. Am I right to think pure recursive does not limit the routines' usability any more than just pure, e.g., in do concurrent or where constructs?

Changes to `simps` and `simps_weights` optional argument `even`

In the current spec, the optional argument even is an integer that is supposed to let the user specify the index at which to insert a Simpson's 3/8 rule so that even-length arrays can be treated. In implementation, this introduced too many corner cases to be worthwhile. In this implementation, even is still an integer, but its meaning is different. If it is negative, the 3/8 rule is used for the leftmost abscissas. If it is positive, the 3/8 rule is used for the rightmost abscissas. If it is zero or absent, the result is as if the even<0 and even>0 cases were averaged. There are still corner cases and complicated control flow, but it's much simpler than the alternative. The spec has been revised to reflect this new behavior.

Improved `fypp` style

Looking toward future extensions to complex types and multiple ranks, the fypp style has been changed to more closely match that of other heavily preprocessed modules.

`trapz_weights` example "bugfix"

In the spec for trapz_weights, the example had a typo and also used dot_product to evaluate the integral. Looking toward implementation for complex integrands, this seems wrong (introduces a complex conjugate), so the example is changed to use sum instead.

…nto dev-quadrature

nshaffer · 2020-06-03T06:12:04Z

~~It looks like the error tolerance is too strict in some of the tests in the CI. Of course they all pass on my machine ;)~~ (resolved)

jvdp1

Thank you for this implementation. I am not the target user for such functions. However, I had some minor comments. If they are inappropriate, please ignore them.

Note: For an easier review, you could maybe submit the changes of trapz and trapz_weights in another PR (not needed for me though).

src/tests/quadrature/test_simps.f90

src/stdlib_experimental_quadrature.fypp

src/stdlib_experimental_quadrature_simps.fypp

Co-authored-by: Jeremie Vandenplas <[email protected]>

…nto dev-quadrature

nncarlson · 2020-06-03T14:44:44Z

src/stdlib_experimental_quadrature.fypp

+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure module function trapz_x_${k1}$(y, x) result(integral)
+            ${t1}$, dimension(:), intent(in) :: y
+            ${t1}$, dimension(size(y)), intent(in) :: x


Is this a private auxiliary function where you know that the actual argument for x will be contiguous in practice, or is this a public user-callable function? (I haven't looked carefully at the whole component to know.) If it is the latter, I don't think you want to be declaring x as an explicit-shape array, dimension(size(y)). This implies F77-style argument passing, and if the actual argument isn't contiguous it will do an implicit copy-in of that argument to a contiguous temporary and pass it instead. Clearly x and y need to have the same size, but that style of declaration does nothing whatsoever to ensure that that is the case (in fact it hides the true length of the actual argument). This really requires a run-time check or assertion, with x declared as an assumed-shape array.

This routine is public by way of the generic trapz interface, so your comment is germane. I have no complaints about using assumed-shape declaration here, but I need to educate myself more on array argument passing to understand why it matters. It is a little surprising to learn that the explicit-shape declaration introduces a temporary for non-contiguous actual arguments, whereas the assumed-shape declaration does not.

What is described is just the implementation detail of (a) specific compiler(s). Although I didn't have the time to check, I strongly doubt that any of such considerations are present in the standard.

It's quite true that this isn't explicitly stated in the standard, however it is implied. The array ${t1}$, dimension(size(y)), intent(in) :: x is simply contiguous (6.5.4) whereas the corresponding actual argument need not be. The only way for that to work is for the array data to be reorganized into contiguous storage.

Are you sure that the only possible way to match a non-contiguous array with a contiguous one is by implicit copy-in? I don't know the answer.

Thanks for the recommendation! I like the approach of making the in-argument assumed shape, and then later benchmarking against making it contiguous as well.

A related questions is if function results should also be assumed-shape. E.g.,

function trapz_weights_dp(x) result(w) real(dp), dimension(:), intent(in) :: x real(dp), dimension(size(x)) :: w ! explicit-shape result ! real(dp), dimension(:), allocatable :: w ! assumed-shape result ... end function trapz_weights_dp

The discussion in #114 suggests that allocatable function results are undesirable. Do I understand right that an explicit-shape result is still expected to be the better option?

I'm not sure there was any conclusion in #114 as to which performed better. The issue probably being whether the array was allocated on the stack or the heap. I haven't looked to see how you're using this particular function, but personally I've switched to generally avoiding array-valued functions altogether, as I feel this results in a lot of unnecessary allocations/deallocations of memory that passing the result as an argument avoids. I'm not really confident that compiler optimizers are able to eliminate these.

Strictly speaking, there is no such thing as an "assumed-shape" function result. It's just an allocatable.

@klausler, thanks for that excellent advice and explanation. @milancurcic we should put this recommendation somewhere in our tutorial section or FAQ, as this question will keep coming up often.

I have just pushed a commit that makes the x argument of trapz_x_* and simps_x_* assumed-shape with a runtime check on its size. Thanks, @nncarlson, for raising the point, and @epagone and @klausler for the discussion.

jvdp1 · 2020-06-03T14:46:02Z

Thanks for the answers. I have additional minor comments. In the specs, could you:

remove "to be implemented" in the titles?
replace simps in result = simps... by something like result = [[stdlib_experimental_quadrature(module):simps(interface)]](....) for the docs?

nshaffer · 2020-06-04T14:10:14Z

Thanks for the answers. I have additional minor comments. In the specs, could you:

* remove "to be implemented" in the titles?

* replace `simps` in `result = simps...` by something like `result = [[stdlib_experimental_quadrature(module):simps(interface)]](....)` for the docs?

Sure, no problemo. I see that there are now also links for the docs in stdlib_experimental_quadrature.fypp for the trapz and trapz_weights interfaces. I suppose we want the same for simps and simps_weights, but I don't know how to get the URL right. Do you know what to do?

jvdp1 · 2020-06-04T14:22:15Z

I suppose we want the same for simps and simps_weights

It would be good, indeed. You can use a similar link as for trapz:

        !! Integrates sampled values using trapezoidal rule
        !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))

The word description must be adapted (e.g., description_2)

nshaffer · 2020-06-04T14:48:28Z

I suppose we want the same for simps and simps_weights

It would be good, indeed. You can use a similar link as for trapz:
        !! Integrates sampled values using trapezoidal rule
        !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))	
The word description must be adapted (e.g., description_2)

Cool beans, made sense after I installed FORD and built the docs. Should be taken care of now.

jvdp1 · 2020-06-04T18:09:57Z

Cool beans, made sense after I installed FORD and built the docs. Should be taken care of now.

Note that you can load it from here.

certik

This generally looks good to me, I think the API works.

After all feedback and comments are addressed, I am fine with merging.

jvdp1 · 2020-06-07T14:46:41Z

src/stdlib_experimental_quadrature_simps.fypp

+        ${t1}$ :: a, b, c
+
+        n = size(y)
+        if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size."


I suggest you use the stdlib subroutine error_stop to avoid issues with older Fortran standards

Suggested change

if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size."

if (size(x) /= n) call error_stop("simps: Arguments `x` and `y` must be the same size.")

The functions cannot be pure when calling error_stop

I think this is a defect of our error_stop. The only functionality it adds over a bare error stop statement is to optionally put both a message and an error code. Because that message is written to stderr, purity is lost. At least with a bare error stop, we can keep purity with sufficiently recent compilers (e.g., GCC version 8 or newer). I am perfectly happy to pretend pre-f2018 standards do not exist, but I think that does not align with stdlib's current goals.

The alternative is to revert to declaring x to be dimension(size(y)). This has the potential performance penalties described in other comments, but at least the compiler/runtime will do the error checking and reporting for us without running afoul of purity rules.

I don't think reverting to declaring x to be dimension(size(y)) is a good idea, as mentioned in earlier comments.

It has been agreed that Gfortran <=7 were not supported for stdlib. However, recent versions of other compilers may not support the 2018 error stop yet. I am perfectly ok to use the 2018 error stop, but this will restrict stdlib to very recent compilers.

certik · 2020-06-07T17:55:53Z

I would recommend to recommend to use error_stop and open another issue to get rid of it all over stdlib, if desired.

…

On Sun, Jun 7, 2020, at 10:04 AM, Jeremie Vandenplas wrote: ***@***.**** commented on this pull request. In src/stdlib_experimental_quadrature_simps.fypp <#209 (comment)>: > +#:for k1, t1 in REAL_KINDS_TYPES + + pure recursive module function simps_x_${k1}$(y, x, even) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, dimension(:), intent(in) :: x + integer, intent(in), optional :: even + ${t1}$ :: integral + + integer :: i + integer :: n + + ${t1}$ :: h1, h2 + ${t1}$ :: a, b, c + + n = size(y) + if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size." I don't think reverting to declaring `x` to be `dimension(size(y))` is a good idea, as mentioned in earlier comments. It has been agreed that Gfortran <=7 were not supported for stdlib. However, recent versions of other compilers may not support the 2018 `error stop` yet. I am perfectly ok to use the 2018 `error stop`, but this will restrict stdlib to very recent compilers. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#209 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAFAWHJR7OUKKKYDXBNKUDRVO3B3ANCNFSM4NRLGEBA>.

nshaffer · 2020-06-07T23:24:54Z

I would recommend to recommend to use error_stop and open another issue to get rid of it all over stdlib, if desired.
…
On Sun, Jun 7, 2020, at 10:04 AM, Jeremie Vandenplas wrote: @.**** commented on this pull request. In src/stdlib_experimental_quadrature_simps.fypp <#209 (comment)>: > +#:for k1, t1 in REAL_KINDS_TYPES + + pure recursive module function simps_x_${k1}$(y, x, even) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, dimension(:), intent(in) :: x + integer, intent(in), optional :: even + ${t1}$ :: integral + + integer :: i + integer :: n + + ${t1}$ :: h1, h2 + ${t1}$ :: a, b, c + + n = size(y) + if (size(x) /= n) error stop "simps: Arguments x and y must be the same size." I don't think reverting to declaring x to be dimension(size(y)) is a good idea, as mentioned in earlier comments. It has been agreed that Gfortran <=7 were not supported for stdlib. However, recent versions of other compilers may not support the 2018 error stop yet. I am perfectly ok to use the 2018 error stop, but this will restrict stdlib to very recent compilers. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#209 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAAFAWHJR7OUKKKYDXBNKUDRVO3B3ANCNFSM4NRLGEBA.

Ok, done. I've replaced the size check using error_stop with a call to check and dropped the pure attributed from simps_x_* and trapz_x_*.

jvdp1

Thank you. +1 to merge

jvdp1 · 2020-06-15T21:21:14Z

All comments were answered. I'll merge it.
Thank you @nshaffer

certik · 2020-06-15T21:22:17Z

@nshaffer do you want to open a new issue for the "error_stop" thing we discussed above?

jvdp1 · 2020-06-15T21:30:41Z

@nshaffer do you want to open a new issue for the "error_stop" thing we discussed above?

Related to #212

nshaffer added 8 commits February 10, 2020 07:32

attempt at quad API

f73beda

improve fypp usage

eb3b22f

dot_product->sum in trapz_weights example, looking toward complex

b7b273b

change semantics of even in simps* spec

d1db8a7

implement Simpson's rule

cd4ef20

Merge branch 'dev-quadrature' of https://github.com/nshaffer/stdlib i…

b2c6264

…nto dev-quadrature

Merge branch 'master' of https://github.com/fortran-lang/stdlib

85f3df4

update assert -> check

178a897

jvdp1 reviewed Jun 3, 2020

View reviewed changes

nshaffer and others added 4 commits June 3, 2020 08:20

relax tolerance in simps tests

688c81a

Co-authored-by: Jeremie Vandenplas <[email protected]>

remove experimental weighted quad stuff

655693e

revise case (6) comment in simps_x_*

031ecf7

Merge branch 'dev-quadrature' of https://github.com/nshaffer/stdlib i…

15eee77

…nto dev-quadrature

nncarlson reviewed Jun 3, 2020

View reviewed changes

remove "to be implemented" from simps specs

321c9d0

add examples for simps and simps_weights

a57b0cb

add doc links

59abe88

certik approved these changes Jun 5, 2020

View reviewed changes

nshaffer added 4 commits June 6, 2020 07:39

make x arrays assumed-shape w/ runtime size check

f83a935

Merge branch 'master' of https://github.com/fortran-lang/stdlib

7961658

remove unused variables

248d92a

Merge branch 'master' into dev-quadrature

856e9be

jvdp1 reviewed Jun 7, 2020

View reviewed changes

{simps,trapz}_x_*: error stop -> check; drop pure

9b6c3d0

jvdp1 approved these changes Jun 8, 2020

View reviewed changes

jvdp1 merged commit 0cd354a into fortran-lang:master Jun 15, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Quadrature: spec tweaks & Simpson's rule #209

Quadrature: spec tweaks & Simpson's rule #209

nshaffer commented Jun 3, 2020

nshaffer commented Jun 3, 2020 •

edited

Loading

jvdp1 left a comment

nncarlson Jun 3, 2020

nshaffer Jun 3, 2020

epagone Jun 3, 2020

nncarlson Jun 3, 2020 •

edited

Loading

epagone Jun 3, 2020

nshaffer Jun 4, 2020

nncarlson Jun 4, 2020

klausler Jun 4, 2020

certik Jun 4, 2020

nshaffer Jun 6, 2020

jvdp1 commented Jun 3, 2020

nshaffer commented Jun 4, 2020

jvdp1 commented Jun 4, 2020

nshaffer commented Jun 4, 2020

jvdp1 commented Jun 4, 2020

certik left a comment

jvdp1 Jun 7, 2020

nshaffer Jun 7, 2020

jvdp1 Jun 7, 2020

certik commented Jun 7, 2020 via email

nshaffer commented Jun 7, 2020

jvdp1 left a comment

jvdp1 commented Jun 15, 2020

certik commented Jun 15, 2020

jvdp1 commented Jun 15, 2020

	if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size."
	if (size(x) /= n) call error_stop("simps: Arguments `x` and `y` must be the same size.")

Quadrature: spec tweaks & Simpson's rule #209

Quadrature: spec tweaks & Simpson's rule #209

Conversation

nshaffer commented Jun 3, 2020

Implementation of simps and simps_weights

Changes to simps and simps_weights optional argument even

Improved fypp style

trapz_weights example "bugfix"

nshaffer commented Jun 3, 2020 • edited Loading

jvdp1 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nncarlson Jun 3, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jvdp1 commented Jun 3, 2020

nshaffer commented Jun 4, 2020

jvdp1 commented Jun 4, 2020

nshaffer commented Jun 4, 2020

jvdp1 commented Jun 4, 2020

certik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

certik commented Jun 7, 2020 via email

nshaffer commented Jun 7, 2020

jvdp1 left a comment

Choose a reason for hiding this comment

jvdp1 commented Jun 15, 2020

certik commented Jun 15, 2020

jvdp1 commented Jun 15, 2020

Implementation of `simps` and `simps_weights`

Changes to `simps` and `simps_weights` optional argument `even`

Improved `fypp` style

`trapz_weights` example "bugfix"

nshaffer commented Jun 3, 2020 •

edited

Loading

nncarlson Jun 3, 2020 •

edited

Loading