Should we not aim at PURE procedures in stdlib?

1

I was just trying to use pdf_norm from the stdlib until my compiler kindly gave me the bad news that this function is impure… As such, it can’t be used in pure procedures and do concurrent loops, which, IMO, significantly limits its practical application. Probably, the same occurs in many other stdlib modules/procedures.
As far as I can see it, the impure attribute comes from the choice of calling the impure procedure error_stop()… Hence, my question:

  • Is the “benefit” of calling the impure procedure error_stop() – versus what could be achieved with a pure error stop message – really worth the cost of downgrading the whole procedure to impure?
    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
    !
    ! Normal distribution probability density function
    !
        ${t1}$, intent(in) :: x, loc, scale
        ${t1}$ :: res
        ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$))
        if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal"       &
            //"distribution scale parameter must be non-zero")
        res = exp(- 0.5_${k1}$ * ((x - loc) / scale) * (x - loc) / scale) /    &
              (sqrt_2_Pi * scale)
    end function pdf_norm_${t1[0]}$${k1}$
    #:endfor
4 Likes
2

I don’t think so. I don’t recall why it was implemented like this, but I suspect it had to do with supporting an older compiler version or OS that didn’t provide the F2018 pure error stop statement. Maybe it’s time to re-evaluate this.

3 Likes
3

This type of statements is common if the dummy arguments need to be restricted. The question is if there is a mechanism to handle the situation so that the target subroutine is pure while the restriction is still in force.

It is interesting to know error stop statement in F2018 is pure. Not sure it is implemented in any compiler.

4

The latest versions of both ifort and gfortran allows the use of error stop in pure procedures.

2 Likes
5

To answer the question in the title, generally yes. However, writing robust and safe-to-use mathematical routines can be hard, and the people behind stdlib are currently just volunteers.

If this were a compiler built-in routine, you would expect it perform a domain check for real literal values, e.g.

! print_sqrt.f90
print *, sqrt(-1.0)
end

$ gfortran print_sqrt.f90 
print_sqrt.f90:1:14:

    1 | print *, sqrt(-1.0)
      |              1
Error: Argument of SQRT at (1) has a negative value

But at runtime you can’t be sure what you get:

! print_cmd_sqrt.f90
character(len=64) :: rstr
real :: r

call get_command_argument(1,rstr)

read(rstr,*) r
print *, "sqrt(",r,") = ", sqrt(r)
end

$ gfortran print_sqrt.f90 -o print_cmd_sqrt
$ ./print_cmd_sqrt 1.0
 sqrt(   1.00000000     ) =    1.00000000    
$ ./print_cmd_sqrt 0.5
 sqrt(  0.500000000     ) =   0.707106769    
$ ./print_cmd_sqrt 0.0
 sqrt(   0.00000000     ) =    0.00000000    
$ ./print_cmd_sqrt -0.0
 sqrt(  -0.00000000     ) =   -0.00000000    
$ ./print_cmd_sqrt -1.0
 sqrt(  -1.00000000     ) =               NaN

What happens with the normal distribution function as the standard deviation tends to zero? Intuitively, it becomes a Dirac delta, so the result should be x-loc (maybe double check my reasoning). However if scale = 0.0 we get division by zero, which will be not-a-number. What happens in case both the numerator and denominator are zero? What about infinity? What about negative zero? What do you do in case of underflow or overflow?

Working out the correct behavior for all of the possible edge cases and writing checks for all of them is a lot of painful work… Not to mention the test harness. Now some user might come along and say, he/she is sure he/she will never have scale == 0.0 and the extra checks are slowing him/her down. Suddenly we need a second “unsafe” version of the function (“fast math”), with a new specification and set of tests.

What started as a question of purity, is actually a much tougher issue of working out the acceptable behavior for each possible floating point value. If you just want to move on with your work, copy the file into your project, delete the scale check, and switch impure to pure.

6

Unless there have been other conversations, the only point raised by OP appears to be only about the PURE attribute of stdlib procedures. The considerations with floating-point behavior with binary arithmetic and robustness under various extreme scenarios need not affect the PURE nature of the methods unless certain IEEE support is sought in which case Fortran 202X support may be needed that extends further their use in PURE procedures.

Given the updates in processors toward the Fortran 2018 standard revision that allows ERROR STOP in PURE procedures, perhaps OP too can extend some help with stdlib and initiate the work that leads to a PR for the stdlib subprograms as applicable?

2 Likes
7

I have the impression we are somehow mixing oranges and potatoes.

  • What is the “correct” value of the normal distribution for sigma=0? That was far from being my question… We can nevertheless ask a math pro! I would probably dare to say that pdf_norm(x, mean, sigma=0) equals +Inf for x->mean and 0 elsewhere. Just for info, Octave gives NaN everywhere (but that is really not the subject).
  • Do robust procedures require error management? Yes, for sure, they do!
  • Does error management require a call to an impure custom error procedure? IMO, not really… we can achieve the same result with error stop message, which almost certainly was introduced in the standard to avoid this kind of dilemma. :slight_smile:
  • Can I make my own pdf_norm procedure? Yes, but I don’t think we want to go down this road… it defies the whole purpose of having a standard library for high-performance Fortran. (And, actually, I already have my own function – I was trying to switch to stdlib to have something better).

Above all, my point is that if stdlib’s goal is to get wide practical acceptance by the community, we can’t afford to downgrade a procedure to impure based on the small added benefit of having a custom error management function.

2 Likes
8

I have only experience with the latest gfortran, where it works like a charm:

 pure subroutine rktvd(fu, u, t, tout, dt, order, itask, istate)
...
    character(:), allocatable :: errmsg

    ! Check input conditions
    if (is_done(t, tout, dt)) return

    if (istate == 1) then
       if (order < 1 .or. order > 3) then
          errmsg = "Invalid input 'order' in 'rktvd'. Valid range: 1 <= k <= 3."
          error stop errmsg
       end if
       if (itask < 1 .or. itask > 2) then
          errmsg = "Invalid input 'itask' in 'rktvd'. Valid set: {1, 2}."
          error stop errmsg
       end if
       istate = 2
    else if (istate < 1 .or. istate > 2) then
       errmsg = "Invalid value 'istate' in 'rktvd'. Valid set: {1, 2}."
       error stop errmsg
    end if
...
1 Like
9

If I can help, I am glad to do so.

10

My impression was the procedure errors out for the value 0.0 because this is easier than dealing with the extreme scenario.

In favor of branchless code execution, you could also get rid of the scale check (and make the procedure pure), and make it the responsibility of the caller to provide valid input. :man_shrugging: I think it’s an important debate what do you want stdlib to prioritize. But of course you can argue I’m digressing at this point… sqrt doesn’t stop your program for an invalid input. Why should the Gaussian function do so?

11

Seems to work well with GFortran back to 5.5.0, so the GFortran version for this particular feature seems not an issue.

❯ gfortran-5 --version
GNU Fortran (GCC) 5.5.0
Copyright (C) 2015 Free Software Foundation, Inc.

GNU Fortran comes with NO WARRANTY, to the extent permitted by law.
You may redistribute copies of GNU Fortran
under the terms of the GNU General Public License.
For more information about these matters, see the file named COPYING

❯ cat test.f90
call test
contains
  pure subroutine test
    character(:), allocatable :: msg
    msg = "Hello, world!"
    error stop msg
  end subroutine test
end
❯ gfortran-5 -Wl,-rpath,/usr/lib/gcc/x86_64-pc-linux-gnu/5.5.0 test.f90
❯ ./a.out
ERROR STOP Hello, world!

The original reason for having call error_stop instead of error stop was the support of exit codes, which are not that well supported. In this case we don’t actually make use of an exit code and therefore we can get the same functionality with a plain error stop.

Getting rid of the check altogether would be an alternative option, in this case we should try to produce at least predictable output in case invalid arguments are provided.

4 Likes
12

This error handling situation is always tricky. Do you want to stop on error, or return an error code, or return some in-band value(s) that require the calling program to test for errors. There are situations where a programmer might want different choices at different times.

Perhaps the solution is to have two different libraries, one that stops on errors or returns error codes and one that returns some in-band value and continues on (requiring the calling program to test).

In the case of a gaussian function evaluation, the main problem seems to be overflow, which can occur for both x==0.0 and x.ne.0.0. Maybe it is a good idea to separate those two cases, but why? Overflow is overflow. Maybe the width==0.0 case is a little special, but not much, it seems to be the same problem. So for in-band errors on an IEEE machine, maybe +INF is the right return value in those overflow cases. On other machines, then HUGE() maybe is the next best thing. Those values could be returned from a pure procedure.

1 Like
13

I would vote for keeping the DRY principle in mind and avoid multiple implementations. I also don’t think we need to reinvent the wheel; we can get inspiration for what other codes do. This is how octave reacts:

>> pkg load statistics

>> s = [1, 2, 0, 3];
>> normpdf(0,1,s)
ans =

   0.2420   0.1760      NaN   0.1258

>>

IMO, this is the expected (or less surprising) behavior: it computes where it can, without stopping the whole process…

3 Likes
14

Ok, I decided to give it a try (HugoMVale/stdlib: Fortran Standard Library (github.com)):

  • converted to pure
  • commented out scale check

It seems to be doing what it is supposed to:

program test_pdf_norm
   use stdlib_stats_distribution_normal, only: pdf_normal
   implicit none

   print *, pdf_normal(0., 1., 1.)
   print *, pdf_normal(0., 1., 2.)
   print *, pdf_normal(0., 1., 0.)

end program

PS ~> fpm run --example
test_pdf_norm.f90                      done.
test_pdf_norm.exe                      done.
[100%] Project compiled successfully.
  0.241970733    
  0.176032662    
              NaN
PS ~>
1 Like
15

pdf_norm looks like may involve generating random numbers. Basically anything involves using random number generator, is mostly likely have some k th step and k+1 th step dependence, and not pure, and cannot be put inside things like do concurrent.

17

pdf_norm is purely deterministic. It can be made PURE, and I’ve done so. Your confusion probably results from the fact that the (same) module stdlib_stats_distribution_normal also contains rvs_normal (Normal Distribution Random Variates).

18

I’m still not sure that NaN is the correct return condition for scale==0.0 and x==loc. Wouldn’t +INF be more consistent? Also, what should happen in the scale==0.0 case when x.ne.loc? Should that be simply 0.0, with no exception or error flag?

19

In the limit of scale->0, we obtain a dirac distribution delta(x-loc). So, the exact mathematical answer should indeed be +inf at x=loc, and 0 elsewhere. But I can hardly believe that a user would intentionally set scale=0. And that is not all – we are actually missing half of the real domain! What if a user sets a negative scale?
That is why, IMO, the best compromise is to return NaN if scale<=0. This clearly tells the user that something is not ok, without stopping the whole calculation.
I will probably propose a PR and then those interested in the topic can continue to exchange there.

stdlib/stdlib_stats_distribution_normal.fypp at stats_dist_normal_pure · HugoMVale/stdlib (github.com)

   #:for k1, t1 in REAL_KINDS_TYPES
      elemental function pdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
         !
         ! Normal distribution probability density function
         !
         ${t1}$, intent(in) :: x, loc, scale
         ${t1}$ :: res
         ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$*acos(-1.0_${k1}$))
         if (scale <= 0._${k1}$) then
            res = ieee_value(1._${k1}$, ieee_quiet_nan)
         else
            res = exp(-0.5_${k1}$*((x - loc)/scale)*(x - loc)/scale)/ &
                  (sqrt_2_Pi*scale)
         end if
      end function pdf_norm_${t1[0]}$${k1}$
   #:endfor
2 Likes