Addition of a mask to stdlib_experimental_stats function `mean`

[jvdp1]

Current

(Related to #113)

The current spec for stdlib mean is:

mean - mean of array elements

Description

Returns the mean of all the elements of array, or of the elements of array along dimension dim if provided.

Syntax

result = mean(array)

result = mean(array, dim)

Arguments

array: Shall be an array of type integer, or real.

dim: Shall be a scalar of type integer with a value in the range from 1 to n, where n is the rank of array.

Return value

If array is of type real, the result is of the same type as array.
If array is of type integer, the result is of type double precision.

If dim is absent, a scalar with the mean of all elements in array is returned. Otherwise, an array of rank n-1, where n equals the rank of array, and a shape similar to that of array with dimension dim dropped is returned.

Proposal

I would like to propose to add mask into the API of mean as follows (similar to the intrinsic sum):

mean - mean of array elements

Description

Returns the mean of all the elements of array, or of the elements of array along dimension dim if provided, and if the corresponding element in mask is true.

Syntax

result = mean(array [, mask])

result = mean(array, dim [, mask])

Arguments

array: Shall be an array of type integer, or real.

dim: Shall be a scalar of type integer with a value in the range from 1 to n, where n is the rank of array.

mask (optional): Shall be of type logical and either be a scalar or an array of the same shape as array.

Return value

If array is of type real, the result is of the same type as array.
If array is of type integer, the result is of type double precision.

If dim is absent, a scalar with the mean of all elements in array is returned. Otherwise, an array of rank n-1, where n equals the rank of array, and a shape similar to that of array with dimension dim dropped is returned.

If mask is specified, the result is the mean of all elements of array corresponding to true elements of mask. If every element of mask is false , the result is zero.

This definition of mask for mean agrees with the definition of mask of the intrinsic sum in the Fortran Standard draft:

p. 429:

MASK (optional) shall be of type logical and shall be conformable with ARRAY.

and

p. 9:

3.36
conformable
⟨of two data entities⟩ having the same shape, or one being an array and the other being scalar

Implementation

Each function for mean, e.g.,

module function mean_1_sp_sp(x) result(res)
    real(sp), intent(in) :: x(:)
    real(sp) :: res
    res = sum(x) / real(size(x, kind = int64), dp)
end function mean_1_sp_sp

should then become:

module function mean_1_sp_sp(x, mask) result(res)
    real(sp), intent(in) :: x(:)
    real(sp) :: res
    logical, intent(in), optional :: mask

    if(.not.optval(mask, .true.))then
        res = 0._dp
	return
    end if

    res = sum(x) / real(size(x, kind = int64), dp)

end function mean_1_sp_sp

module function mean_1_mask_sp_sp(x, mask) result(res)
    real(sp), intent(in) :: x(:)
    real(sp) :: res
    logical, intent(in) :: mask(:)
    res = sum(x, mask) / real(count(mask, kind = int64), dp)
end function mean_1_mask_sp_sp

I kown there was some discussions about dim and if it should be considered as optional or not in the spec. I guess the same discussion could appear here since mask is not optional in one of the 2 functions. However, it will behave like optional for the users through the interface mean.

[certik]

That looks good to me. It would appear as "optional" due to the interface to the users, so I think that's what matters, and the above is just an implementation detail how to actually implement it.

[certik]

A larger philosophical question is whether we have to implement mask for every array function that operates on arrays, or only on some?

[jvdp1]

A larger philosophical question is whether we have to implement mask for every array function that operates on arrays, or only on some?

Why would you not implement it on everey function? Also, from rank 3 to rank 7 or 15, it would auto-generated.

[certik]

The only downside that I can see is that it will be a lot of functions. But I think we want to support the "mask" argument. The other approach would be that somehow the user would run the mask on the array before passing it to sum or mean (NumPy implements masked arrays). But I don't know how to do that efficiently in Fortran, so the mask argument seems appropriate.

[jvdp1]

I think it would be good to have the same API as sum and other intrinsic functions.
The users will not see all the functions behind the interface, and I think it can be manageable with auto-generation. It could be only (maybe) tricky for the CI, but we can limit the number of dimensions (currently 4). Also, I don't think the CI should stop us to implement a correct/desired code/API.

[certik]

I agree.

[epagone]

Excellent proposal. I only have a problem with the following statement:

If every element of mask is false , the result is zero.

Is this not dangerous? The user might think that zero is a legitimate mean value. Can't we return (a IEEE?) NaN or something along these lines so the user can catch the exception? All this hopefully without hurting performance.

[jvdp1]

Is this not dangerous? The user might think that zero is a legitimate mean value. Can't we return (an IEEE?) NAN or something along these lines so the user can catch the exception? All this hopefully without hurting performance.

Good point. The intrinsic sum returns 0, following the Standard:

Case (ii): The result of SUM (ARRAY, MASK = MASK) has a value equal to a processor-dependent approximation to the sum of the elements of ARRAY corresponding to the true elements of MASK or has the value zero if there are no true elements.

So, for the mean, if mask = .false., it ends up with a division of 0. by 0.. So mean should most likely return NaN. A simple way to return NaN could be:

real(dp)::res
if (.not.optval(mask, .true.) then
    res = res /0._dp
    return
end if

Note that GFortran does not allow expressions like res = 0._dp / 0._dp or res = 1._dp / 0._dp.

[certik]

I would be careful with the nans. I usually turn them off in my production codes. Is there any current intrinsic function in Fortran that can return nan? I think returning zero would work. Stopping the program with an error would also work, but the issue is that then mean could not be a pure function...

[epagone]

I would be careful with the nans. I usually turn them off in my production codes.

I think this resonates in some way with the discussion had about assert. I.e. in debug mode you catch the cases that return NaN and in release mode you can ignore them (at your own risk, trading "safety" with performance).

Still, I think mathematically stands that although the sum of no values is zero, the mean of no values is undefined.

[urbanjost]

Having had to deal with Nan in a lot of programming environments, the most reliable way I found may seem surprising, but is

function nan64(value)
use,intrinsic :: iso_fortran_env, only: real64
implicit none
character(len=*),parameter::ident="@(#)M_units:: nan64(3fp): Returns a NAN (Not a number) of type real64"
character(len=3),save :: STRING='NaN'
real(kind=real64) :: nan64,value
   read(STRING,*)nan64
end function nan64

which is one part of a generic nan(value) where the returned value is the type of value that could be auto-generated. Note the string cannot be a parameter, which seems non-intuitive; but internal reads cannot read from a constant. The rest of the generic definition is quite repetitive but (I think) obvious.

It's a long story, but a lot of "standard" ways of defining a NaN are subject to whether the compiler supports the IEEE standard, whether optimizations change the definition or test (I have an is_nan function too), and so on. So far, reading the string "NaN" has proved portable and is supported by the Fortran standard as far back as I could find (if anyone knows when the earliest support by the standard was, feel free to add it here).

[certik]

@urbanjost interesting. A few questions: what's the purpose of the ident parameter in your nan64 function? Also, why does it accept value as an argument? It does not seem to be used.

[jvdp1]

Still, I think mathematically stands that although the sum of no values is zero, the mean of no values is undefined.

I totally agree with @epagone statement, and I should have thought about it sooner: the mean of no values is undefined, not 0.

Regarding our implementation of mean along a dimension and with a mask being of the same shape as the input x, we have:

module function mean_2_mask_sp_sp(x, dim, mask) result(res)
    real(sp), intent(in) :: x(:,:)
    logical, intent(in) :: mask(:,:)
    integer, intent(in) :: dim
    real(sp) :: res(merge(size(x, 1), size(x, 2), mask = 1 < dim ))

    select case(dim)
     case(1)
       res = sum(x, 1, mask) / real(count(mask, 1), sp)
     case(2)
       res = sum(x, 2, mask) / real(count(mask, 2), sp)
     case default
       call error_stop("ERROR (mean): wrong dimension")
    end select

end function mean_2_mask_sp_sp

This implementation may return an array res with numbers and NaN as means.

If we aim to return 0.0 (or any other value different thanNaN (even for the NaN proposed by @urbanjost) for the following case:

module function mean_2_all_sp_sp(x, mask) result(res)
    real(sp), intent(in) :: x(:,:)
    logical, intent(in), optional :: mask
    real(sp) :: res

    if (.not.optval(mask, .true.)) then
        res = ??? 0._sp / value / @urbanjost NaN????
        return
    end if

    res = sum(x) / real(size(x, kind = int64), sp)

end function mean_2_all_sp_sp

then, for consistency across all mean functions, we will have to catch the possible compiler NaN in mean_2_mask_sp_sp(x, dim, mask) and replace them by ??? value / @urbanjost NaN????, leading to complicate and less efficient code.

Probably, IMHO, the best way to be consistent with mean_2_mask_sp_sp(x, dim, mask) (and to keep this latter function simple and efficient) would be something like:

module function mean_2_all_sp_sp(x, mask) result(res)
    real(sp), intent(in) :: x(:,:)
    logical, intent(in), optional :: mask
    real(sp) :: res

    if (.not.optval(mask, .true.)) then
        res = sum(x, .false.) / count(.false.)
        return
    end if

    res = sum(x) / real(size(x, kind = int64), sp)

end function mean_2_all_sp_sp

With such an implementation, NaN could be still ignore in release mode (while not the case with @urbanjost approach, right?).