Proposal for clamp

[milancurcic]

Intro

clamp(x, xmin, xmax) is a function that clamps (or, clips) an input scalar or array x so that the result is within given bounds xmin and xmax. For example:

res = clamp(0.9. 0.2, 0.8) ! res is 0.8
res = clamp(0.7. 0.2, 0.8) ! res is 0.7
res = clamp(0.1, 0.2, 0.8) ! res is 0.2

This function is common in standard libraries. I've used it occasionally in Fortran, and I use it often in Python. It provides an easy way to ensure that a variable with a finite (non-NaN, non-Inf) value stays within bounds.

Name

Python has numpy.clip, while C++, Julia, and Rust all have clamp. So clamp is more common. Both clamp and clip are equally meaningful to me. I have slight preference for clip because it's shorter and has a "sharper" sound to it.

• numpy.clip
• C++ clamp
• Julia clamp
• Rust clamp

What inspired me to open this proposal is reading that clamp was stabilized in Rust 1.50.0.

API

elemental real function clamp(x, xmin, xmax) result(res)
  real, intent(in) :: x ! input value to clip
  real, intent(in) :: xmin ! lower bound
  real, intent(in) :: xmax ! higher bound

Alternative names for the bounds could be low and high.

Implementation

There are a few options that I know of, each of which produce different assembly instructions, and performance depending on optimization levels, array sizes, and whether the clamping is in a loop or over a whole array in a single call.

A min/max clamp

This is the simplest implementation I could think of:

elemental real function clamp(x, xmin, xmax) result(res)
  real, intent(in) :: x, xmin, xmax
  res = min(max(x, xmin), xmax)
end function clamp

There are two intrinsic function calls here, but maybe a good compiler will optimize one away. I don't know.

An if/then/else clamp

elemental real function clamp(x, xmin, xmax) result(res)
  real, intent(in) :: x, xmin, xmax
  if (x < xmin) then
    res = xmin
  else if (x > xmax) then
    res = xmax
  else
    res = x
  end if
end function clamp

A branchless clamp

This avoids branching but introduces a floating-point error (~epsilon(x)) for some elements, even if x is not clipped.

elemental real function clamp(x, xmin, xmax) result(res)
  ! A branchless clamp.
  ! Credit: Mark Ransom (https://stackoverflow.com/a/7424900/827297)
  real, intent(in) :: x, xmin, xmax
  res = (x + xmin + abs(x - xmin)) / 2
  res = (res + xmax - abs(res - xmax)) / 2
end function clamp

What module would this go to?

I don't think we have an appropriate existing module yet, but stdlib_numeric or stdlib_math sound meaningful to me.

[ivan-pi]

Do you think there will be much use for passing arrays to the xmin or xmax arguments?

If not, I would propose a fourth implementation option:

A where clamp

The procedures are auto-generated for each rank and type (real, double, ...)

interface clamp
  module procedure clamp_rank0
  module procudure clamp_rank1
  module procedure clamp_rank2
  ...
end interface

function clamp_rank0(x,xmin,xmax)
! ... use min max or ifthen/else version
end function

function clamp_rank1(x,xmin,xmax) result(res)
  real, intent(in) :: x(:)
  real, intent(in) :: xmin, xmax
  real :: res(size(x))
  where (x < xmin)
    res = xmin
  elsewhere (x > xmax)
    res = xmax
  elsewhere
   res = x
  end where
end function

The same logic applies to where/elsewhere logic applies to higher ranks.

[ivan-pi]

Concerning name I also have a personal preference for clip. On the other hand consistency with C++, Julia, and Rust does have some weight. I would also vote in favor of stdlib_numeric vs stdlib_math.

[milancurcic]

Is the usefulness of a where clamp to avoid it being elemental, like in #310 (comment)?

If yes, is it so that clamp could be passed as a procedure pointer? If yes, it makes me think--how do you do it with a generic? I thought that generic can't be used as a procedure pointer because you need to specify the procedure interface in the list of dummy arguments. Or are there other benefits to non-elemental procedures?

[milancurcic]

Do you think there will be much use for passing arrays to the xmin or xmax arguments?

I don't think much but maybe some. I used it a few times in Python for 2-d arrays to clamp only within a radius and to do a masked clamp. There are alternative ways to do that though.

But this could also be provided through specific functions that define xmin and xmax as arrays. So if we chose to do a generic implementation instead of elemental, then allowing xmin and xmax to be arrays would only double the number of specific functions.

[ivan-pi]

Is the usefulness of a where clamp to avoid it being elemental, like in #310 (comment)?

No, in this context I was only playing with the idea of a generic but non-elemental interface to forbid usage such as:

res = clamp(0.1,[0.2,0.3,0.4],0.5)

which if the implementation were min(max(x, xmin), xmax) would return the 3-by-1 vector [0.2,0.3,0.4].

I haven't used the clamp function before to be able to say is it useful or not. Just throwing the concern out for discussion.

[milancurcic]

I see, that didn't occur to me. I can't think of an use for it. I agree that we should somehow forbid it or at least discourage it. I see a few solutions:

• Use the generics and the where clamp to do it at compile time (your suggestion)
• Check that xmin and xmax are scalar at run-time and stop if they are not
• Leave unchecked and document it in user docs

[ivan-pi]

I would prefer to avoid run-time checking. This would make the procedures impure so I don't see any benefit. Having the routine return without warning also make no sense.

This leaves us with the first and third bullet points, meaning either pure procedures for different ranks with a generic interface, or the elemental version (either min/max or the if-then-else versions). Actually it seems like a great experiment to learn about potential performance differences. 👍

Am I right that according to the specified behavior positive Inf values are mapped to xmax, while -Inf is mapped to xmin? NaN values however are preserved?

[milancurcic]

I would prefer to avoid run-time checking. This would make the procedures impure so I don't see any benefit. Having the routine return without warning also make no sense.

Yes, I would too. We can error stop from a pure a procedure (F2018) but it'd still need a check which would hinder performance.

[milancurcic]

My heart prefers the elemental approach (it's one of my favorite features of Fortran) but my head prefers the generic where clamp approach (compile-time enforcing of argument shape). Let's discuss it more broadly on the upcoming call.

[ivan-pi]

My heart prefers the elemental approach (it's one of my favorite features of Fortran) but my head prefers the generic where clamp approach (compile-time enforcing of argument shape). Let's discuss it more broadly on the upcoming call.

I wonder if this is something the standard should address - being able to explicitly enforce that elemental affects only certain variables (perhaps through an attribute such as nonelemental or scalar).

[art-rasa]

Hello all,
I was just thinking if there would be any need for a procedure that would scale the value of an input parameter into an output value defined by min/max parameters. This could be used, for example, scaling an input scalar (or array) into a scaled value describing some physical quantity, such as temperature between defined limits (say -50 ℃ and 100 ℃).

[MarDiehl]

My current implementation with optional upper and lower bounds and 'hint' to the user if lower>upper (inspired by numpy). Probably error stop would be better than NaN, but I'm not clear about performance implications.

real(pReal) pure elemental function math_clip(a, left, right)

  real(pReal), intent(in) :: a
  real(pReal), intent(in), optional :: left, right

  math_clip = a
  if (present(left))  math_clip = max(left,math_clip)
  if (present(right)) math_clip = min(right,math_clip)
  if (present(left) .and. present(right)) &
    math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right

end function math_clip

[milancurcic]

@art-rasa Yes, I think so, I've done this many times in the past (home cooked) and is very common in machine learning. Can you propose it in a separate issue? Ideally with an API, name suggestion, and example implementation.