Discussion of `diff` procedure

[zoziha]

Motivation

Here are some questions about the diff procedure:

• Is it worth introducing diff procedure?
• API of diff procecure?
  Function? Subroutine?
• Return value: allocatbable array or automatic array?

#605 is about the submission of diff, it seems that the efficiency of this implementation needs to be improved. Its API needs to be discussed and reached a consensus.

@gareth-nx has some suggestions in #605. (see here )

• Efficiency
• API
• User-oriented

Prior Art

• diff in Matlab

Y = diff(X)
Y = diff(X,n)
Y = diff(X,n,dim)

• diff in numpy

Y = numpy.diff(a, n=1, axis=-1, prepend=<no value>, append=<no value>)

Note:
prepend, append: array_like, optional: New in numpy version 1.16.0.

Additional Information

No response

[gareth-nx]

The optional arguments prepend and append in the numpy routine give a flexible treatment of the end-points issue raised in #605.

The matlab version doesn't offer that - it provides output that is of reduced size compared to the input. Numpy defaults to similar behaviour (if prepend and append are not given). R also includes a diff routine that produces an output array of reduced size.

IMO something like the numpy version could be quite flexible -- using optional arguments with defaults like the above.

If this were done with a subroutine interface then the problem of output allocation on every call could be avoided. OTOH there is some convenience trade-off, vs a function interface.

[LKedward]

• Is it worth introducing diff procedure?

Yes, I use this a lot in MATLAB.

• API of diff procecure?

IMO, it makes sense to stick closely to the API of prior art with MATLAB and numpy because this is what people will be used to.
I didn't know numpy provided the prepend and append options, they look useful.

Function? Subroutine?

How has this question (function vs subroutine) been approached in the rest of stdlib for math operations?
(The argument has been made elsewhere that the additional overheads of array function calls is simply a question of compiler implementation, which could be solved in the compiler)

[gareth-nx]

How has this question (function vs subroutine) been approached in the rest of stdlib for math operations?

One option would be to expose both function (diff) and subroutine (difs) interfaces? Internally the former could take care of allocation and then call the latter. difs could assume that the user-provided intent(out) array has the right shape.

[jvdp1]

Thank you for this proposition.

• Is it worth introducing diff procedure?

Yes, I think so. E.g., it could be used in a function mode in the module stdlib_stats.

• API of diff procecure?
  I suggest to stay as close as possible to the ones implemented in other languages, like matlab, octave or Julia.

Supporting an optional mask might be an advantage.

Function? Subroutine?

• Return value: allocatbable array or automatic array?

Usually such procedures are implemented as functions with automatic arrays in stdlib.

[zoziha]

Thank you very much for all your feedback~

If we decide to introduce the prepend and append arguments (array_like) in the diff function, it seems that we may need to consider introducing a similar append function to stitch the arrays (prepend, diff and append arrays) ?? :

• numpy.append (which is used in the implement of numpy.diff, see here )
• matlab. horzcat / vertcat / cat

[gareth-nx]

Usually such procedures are implemented as functions with automatic arrays in stdlib.

If we are to support prepend and append as optional arguments then I don't think this will work, as the size of the output arrays will be size(input) + size(prepend) + size(append). I assume we can't use that expression to declare the automatic output array size if prepend and append are optional (?).

If we decide to introduce the prepend and append arguments (array_like) in the diff function, it seems that we may need to consider introducing a similar append function to stitch the arrays (prepend, diff and append arrays) ??

I assumed that, rather than perform array stitching, we could compute the differences associated with prepend and append in a separate part of the code within e.g. if(present(append)) then ... and similarly for prepend. The approach you suggest would be more elegant, but might be slower due to the array stitching.

Both issues above would be avoided if we don't support append and prepend. This is one argument for not supporting those options.

OTOH, diff is simple enough for most users to implement themselves efficiently inline (x(2:N) - x(1:N-1)). Given that, IMO it is preferable for the stdlib implementation to be "possible to use efficiently" (even if not as efficient as user-inlined code in every case) and provide a few conveniences to users. On those grounds I like the idea of append, prepend, and a subroutine interface (or some other interface that can avoid the allocation).

[zoziha]

I tried to build on the previous diff commits and made a poor implementation only for our discussion (see stdlib_math_diff.fypp), adding support for the prepend and append arguments, but I believe this does increase the complexity of the internal implementation of diff.

https://github.com/fortran-lang/stdlib/blob/8561c45325589b1e49a90b5e7b0e804a81786054/src/stdlib_math_diff.fypp#L92~L99 (Outdated)
https://github.com/fortran-lang/stdlib/blob/9ebec0194f163600632eccaaff9c0ff6aafb8a9e/src/stdlib_math_diff.fypp#L78~L82

I still hope to implement functions like numpy.append in stdlib to facilitate stitching arrays.
This will also improve the existing diff implementation, which is not only elegant, but also clean and easy to maintain.

[gareth-nx]

Thanks @zoziha for the discussion implementation.

I wonder if we should only support rank 1 and rank 2 arrays (as in the draft implementation)?

Both matlab and numpy support arrays of any rank. OTOH, R only supports rank1 arrays.

Is there much value to the rank > 1 case (as compared to users doing the required calculation in a loop)? Rather than supporting rank=1,2 but not other ranks, it might be conceptually simpler to just support rank=1, and suggest looping for other ranks (??). Or support all ranks?

[gareth-nx]

FWIW I've been exploring some alternative implementations (very preliminary, only for real, rank=1 arrays). Edited to fix the test code, new timings below.

The gist here contains:

• a build script
• a program run_test.f90 (that also uses the draft stdlib version of diff), and
• a module with an alternative function (diff) and subroutine (difs) implementations, where the latter is inflexible, and the former is like the stdlib rank=1 version.

If you copy them all to the same directory, modify build.sh to include your stdlib environmental variables, then it can be compiled by running source build.sh.

I used gfortran -O3 to compile it (which is the same that I am using for stdlib).

When we use prepend, append I find the timings are often a bit faster (1.5-5.3) for the alternative versions compared to the current stdlib draft, with greater speedups for small arrays. When we don't use prepend, append, the time difference shows the same general trend but with less difference (and in the case with a very large array, the timings of all implementations are fairly close).

Test case with N=    10000000 ; reps=          10 ; size_append_prepend=           1
 Time ratio (stdlib_function / local_subroutine) :   2.43964362    
 Time ratio (stdlib_function / local_function  ) :   1.67561471    
 Test case with N=       10000 ; reps=       10000 ; size_append_prepend=           1
 Time ratio (stdlib_function / local_subroutine) :   4.08620548    
 Time ratio (stdlib_function / local_function  ) :   3.27745557    
 Test case with N=         100 ; reps=     1000000 ; size_append_prepend=           1
 Time ratio (stdlib_function / local_subroutine) :   5.29873419    
 Time ratio (stdlib_function / local_function  ) :   3.70539761    
 Test case with N=    10000000 ; reps=          10 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   1.11118972    
 Time ratio (stdlib_function / local_function  ) :  0.756159544    
 Test case with N=       10000 ; reps=       10000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   1.96070278    
 Time ratio (stdlib_function / local_function  ) :   1.54513323    
 Test case with N=         100 ; reps=     1000000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   2.54248381    
 Time ratio (stdlib_function / local_function  ) :   1.45795929    

Very much an experiment, but maybe something is of use.

[zoziha]

The operating efficiency of subroutines is higher than that of functions. I think array functions are more suitable for creating a new array and data, and subroutines are more suitable for overwriting and refreshing array data.

From this perspective, creating a subroutine difs is efficient and demanding?

[gareth-nx]

I re-ran the timing test mentioned above, using the recently pushed version of the pull request.

In the case where we use append/prepend:

• For mid sized arrays (size of 100 - 10000) it still has a noticeable time penalty compared to the function version in the gist, although is faster than the previous stdlib pull-request version.
• For large arrays (10-million entries) it is slightly faster than the function version in the gist, and slightly slower than the inflexible subroutine version.

I wonder if this difference is related to the work array being allocatable in the stdlib version, vs automatic in the gist version? The gist version achieves that by having a worker subroutine with no optional arguments. It would be nice if we could make the stdlib version at least as efficient (or, understand why we can't).

If we find there are persistent speed advantages to the subroutine interface, I would lean toward exposing a subroutine interface (even if it's less convenient) as well as a function interface.

 Test case with N=    10000000 ; reps=          10 ; size_append_prepend=           1
 Time ratio (stdlib_function / local_subroutine) :   1.09848988    
 Time ratio (stdlib_function / local_function  ) :  0.752216458    
 Test case with N=       10000 ; reps=       10000 ; size_append_prepend=           1
 Time ratio (stdlib_function / local_subroutine) :   2.02271032    
 Time ratio (stdlib_function / local_function  ) :   1.64654362    
 Test case with N=         100 ; reps=     1000000 ; size_append_prepend=           1
 Time ratio (stdlib_function / local_subroutine) :   1.98440468    
 Time ratio (stdlib_function / local_function  ) :   1.46916664    
 Test case with N=    10000000 ; reps=          10 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   1.07091498    
 Time ratio (stdlib_function / local_function  ) :  0.734303594    
 Test case with N=       10000 ; reps=       10000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   2.03846216    
 Time ratio (stdlib_function / local_function  ) :   1.62488389    
 Test case with N=         100 ; reps=     1000000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   2.11474347    
 Time ratio (stdlib_function / local_function  ) :   1.43239534    

[gareth-nx]

It's also interesting to compare the performance vs writing(x(2:) - x(1:size(x)-1) in-line , in the common case where prepend, append are not provided, and n=1.

It turns out the inline version is about 6x faster (or in the best cases, 2.5-3x vs the subroutine). Here is the output where inline_without_prepend_or_append refers to the special case.

 Test case with N=    10000000 ; reps=          10 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   1.08094656    
 Time ratio (stdlib_function / local_function  ) :  0.750260413    
 Time ratio (stdlib_function / inline_without_append_prepend_or_n  ) :   6.21810484    
 Test case with N=       10000 ; reps=       10000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   1.96804166    
 Time ratio (stdlib_function / local_function  ) :   1.24381065    
 Time ratio (stdlib_function / inline_without_append_prepend_or_n  ) :   5.96953821    
 Test case with N=         100 ; reps=     1000000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   2.40701604    
 Time ratio (stdlib_function / local_function  ) :   1.19632530    
 Time ratio (stdlib_function / inline_without_append_prepend_or_n  ) :   6.17020941    

However much of the performance can be regained by giving special treatment of that case in the function and subroutine implementation. Using the versions in this gist I get:

 Test case with N=    10000000 ; reps=          10 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   6.24738979    
 Time ratio (stdlib_function / local_function  ) :   1.99576485    
 Time ratio (stdlib_function / inline_without_append_prepend_or_n  ) :   6.25094366    
 Test case with N=       10000 ; reps=       10000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   6.19736195    
 Time ratio (stdlib_function / local_function  ) :   3.57656074    
 Time ratio (stdlib_function / inline_without_append_prepend_or_n  ) :   6.12597227    
 Test case with N=         100 ; reps=     1000000 ; size_append_prepend=           0
 Time ratio (stdlib_function / local_subroutine) :   3.87429070    
 Time ratio (stdlib_function / local_function  ) :   3.07483053    
 Time ratio (stdlib_function / inline_without_append_prepend_or_n  ) :   5.73714638    

Summary:

• There seems to be an argument for implementing special-cases efficiently in the stdlib routines.
• I imagine this approach could be generalised to treat the case where n=1 even when we do have prepend, append (since in that case, one could avoid the work allocation easily enough).