Re: Update and questions regarding vecorization

[Oliver Heimlich]

Hi Joel,

On 21.06.2017 16:55, Joel Dahne wrote:
> This week I have begun the work on vectorization and I have some
> questions/comments.
> 
> Where to test vectorization?:
> At the moment the tests in the unary functions (and some others), for
> example sin.m or cos.m, does not test vectorization in any way. They all
> use the same functions for performing the operations,
> crlibm_functions.cc or mpfr_function_d.cc, so that functionality can
> perhaps be tested there. But a lot of the unary functions also do a lot
> of other computations directly in the octave-script. This is true in for
> example sin.m. Should we add tests for this in their own files? We would
> then test with vectors, matrices and ND-arrays as input. I can also
> mention that all of work done in the octave-scripts uses smart indexing
> so they immediately worked for N-dimensional arrays.

There are some vectorized test cases for mpfr_function_d (see
src/mpfr_function_d.cc) and for crlibm_function (see
inst/__check_crlibm__.m).  They load several test cases in column
vectors and call the tested function only once.  However, it only
affects a few functions and does not involve the interval methods.

Of course, you can add some hand-picked test cases at the end of the
individual functions that you'd like to be tested, e. g., in
inst/@infsup/sin.m you may add the cause of bug #51283.

When you want to check all mapping functions (=functions that are
applied element-wise on non-scalar inputs) on a large scale, I would
suggest that we improve ITF1788:  ITF1788 currently generates a long
list of assert statements to test arithmetic correctness of most
functions in the package with many test cases.  The test cases are
defined in the test/*.itl files.  That could be improved as follows:
ITF1788 generates an Octave script, which could then generate test code
for Octave.  The numeric test data could be stored in *.mat files and
could then be used to (1) loop over the test data for scalar test cases
of interval functions like the assert statements do today, to (2) test
the mapping functions on the vectors of all test values at once, and (3)
reshape the test data into more dimensions to also test vectorization in
higher dimensions, and (4) check broadcasting.

Do you want me to look into this during the weekend?  When you are
unfamiliar with ITF1788, it could be distracting to start this side
project yourself now.

> Broadcasting for ternary functions:
> I added support for broadcasting of N-dimensional arrays for binary
> functions in mpfr_function_d.cc. This in practice required a rewrite of
> the whole function, I will write a bit more about this in my blog post
> tomorrow. When I looked at the ternary functions noticed that it does
> not perform full broadcasting at the moment. However the only ternary
> function used is fma.m and that function performs broadcasting already
> in the octave-script. Should I still add broadcasting for ternary
> functions? It might be used in the future. If we choose not to we should
> probably add a comment about it in the code.

Sounds complicated, I'll answer this during the Weekend after looking at
your code for the binary functions.  I don't see broadcasting in
mpfr_function_d as important, since I consider it more of an internal
function.  Do you know of any ternary Octave functions that allow
broadcasting, e. g., does plot3 allow broadcasting?

> rdivide:
> This is not really connected to N-dimensional arrays but I noticed it
> when reading the file. When the numerator is equal to 1 it short
> circuits and performs a more efficient computation. If the numerator is
> an array it only short circuits if all the elements are 1. We should be
> able to short circuit for all elements that are 1 and for the other ones
> we perform the normal computation. Should I implement this? Since it
> does not have to do with my project it might be a better idea to add it
> directly to the interval repository and not to my repository.

Do it if you like.  I have been too lazy to write an indexing expression
to short-circuit the correct subset of the N-D array.  You may either
send a separate patch (patch tracker) or commit it to your repo, from
where I will eventually merge it anyway.

> sin.m:
> See https://savannah.gnu.org/bugs/index.php?51283. We should probably
> add a work around for this. Or perhaps wait and see what Octave does
> with it?

Thanks for the bug report.  I have added a comment.  Please use the work
around to also support current versions of Octave.

> times.m:
> The choice has been made to compute all possible multiplications at all
> times. Normally I agree with this, but broadcasting makes things
> different. If we have a lot of broadcasting it does make sense to check
> the boundaries first I think, one check can then potentially save lots
> of multiplications. For example in
> 
> infsup (rand (1000, 1000)) .* infsup (1, 2);
> 
> we could save 6 million (correct?) computations with one check. This is
> of course an extreme example but you get the point. Since multiplication
> with a scalar is quite a common operation it might be worth it to add a
> short circuit for this case. It does however make the code more complex,
> which is never good.

Did you measure execution times?  When I last checked it with
non-broadcasting vectors, it wasn't worth the effort (except for corner
cases).  With broadcasting this should be no different.

However, this might have changed or I might be wrong.  You may try to
find a faster solution.  The times function is not that complicated and
we have lots of unit tests for regression testing.  If you find a faster
algorithm that'd be good.  However, you should measure execution times
for several input constellations.

> It will try post on my blog tomorrow. On Friday we celebrate Midsummer's
> Eve here in Sweden so then I will not be working.

I'm looking forward to your post and have a happy Midsommar.

Best
Oliver