## Monday, 28 August 2017

### Final Report

This is the final report on my work with the interval package during Google Summer of Code 2017. This whole blog have been dedicated to the project and by reading all posts you can follow my work from beginning to end.

The work has been challenging and extremely fun! I have learned a lot about interval arithmetic, the Octave and Octave-forge project, and also how to contribute to open-source in general. I have found the whole Octave community to be very helpful and especially I want to thank my mentor, Oliver Heimlich, and co-mentor, Kai Torben Ohlhus, for helping me during the project.

Here I will give a small introduction to Octave and the interval package for new readers and a summary of how the work has gone and how you can run the code I have actually contributed with.

## Octave and the Interval Package

Octave, or GNU Octave, is a free program for scientific computing. It is very similar to Matlab and its syntax is largely compatible with it. Octave comes with a large set of core functionality but can also be extended with packages from Octave forge. These add new functionality, for example image processing, fuzzy logic or more statistic functions. One of those packages is the interval package which allows you to compute with interval arithmetic in Octave.

## Summary of the Work

The goal with the project was to improve the Octave-forge interval package by implementing support for creating, and working with, N-dimensional arrays of intervals.

The work has gone very well and we have just released version 3.0.0 of the interval package incorporating all the contributions I have made during the project.

The package now has full support for working with N-dimensional arrays in the same way you do with ordinary floating point numbers in Octave. In addition I have also fixed some bugs not directly related to N-dimensional arrays, see for example bug #51783 and #51283.

During the project I have made a total of 108 commits. I have made changes to 332 of the packages 666 files. Some of these changes have only been changes in coding style or (mainly) automated adding of tests. Not counting these I have, manually, made changes to 110 files.

If you want to take a look at all the commits I have contributed with the best way is to download the repository after which you can see all the commits from GSoC with

hg log -u joeldahne@hotmail.com -d "2017-06-01 to 2017-08-29"

Unfortunately I have not found a good way of isolating commits from a specific period and author on sourceforge where the package is hosted. Instead you can find a list of all commits at the end of this blog post.

The NEWS-file from the release of version 3.0.0 is also a pretty good overview of what I have done. While not all of the changes are a result of GSoC quite a lot of them are.

## Running the Code

As mentioned above we have just released version 3.0.0 of the interval package. With the new release it is very easy to test the newly added functionality. If you already have Octave installed the easiest way to install the package is with the command “pkg install -forge interval”. This will install the latest release of the package, at the time of writing this is 3.0.0 but that will of course change in the future. You can also download version 3.0.0 directly from Octave-forge.

If you want you can also download the source code from the official repository and test it with "make run" or install it with "make install". To download the repository, update to version 3.0.0 an run Octave with the package on Linux use the following

hg clone http://hg.code.sf.net/p/octave/interval octave-interval
cd octave-interval
hg update release-3.0.0
make run

## A Package for Taylor Arithmetic

The task took less time than planned so I had time to work upon a project depending on my previous work. I started to work on a project for Taylor arithmetic, you can read my blog post about it. I created a proof of concept implementation as part of my application for Google Summer of Code and I have now started to turn that into a real package. The repository can be found here.

It is still far from complete but my goal is to eventually add it as a package at Octave-Forge. How that goes depends mainly on how much time I have to spend on it the following semesters.

If you want to run the code as it is now you can pull the repository and then run it with "make run", this requires that Octave and version 3.0.0 (or higher) of the interval package is installed.

## List of Commits

Here is a list of all 108 commits I have done to the interval package

https://sourceforge.net/p/octave/interval/ci/68eb1b3281c9

https://sourceforge.net/p/octave/interval/ci/52d6a2565ed2
summary:     @infsupdec/factorial.m: Fix decoration (bug #51783)

https://sourceforge.net/p/octave/interval/ci/cf97d56a8e2a
summary:     mpfr_function_d.cc: Cast int to octave_idx_type

https://sourceforge.net/p/octave/interval/ci/5ed7880c917f
summary:     maint: Fix input to source

https://sourceforge.net/p/octave/interval/ci/637c532ea650
summary:     @infsupdec/dot.m: Fix decoration on empty input

https://sourceforge.net/p/octave/interval/ci/51425fd67692

https://sourceforge.net/p/octave/interval/ci/b4e7f1546e36
summary:     mpfr_function_d.cc, mpfr_vector_dot_d.cc: Fixed bug when broadcasting with one size equal to zero

https://sourceforge.net/p/octave/interval/ci/02db5199932c
summary:     doc: NEWS.texinfo: Info about vectorization for nthroot and pownrev

https://sourceforge.net/p/octave/interval/ci/197d033fc878
summary:     @infsup/pownrev.m, @infsupdec/pownrev.m: Support for vectorization of p

https://sourceforge.net/p/octave/interval/ci/7f8d1945264c
summary:     @infsup/nthroot.m, @infsupdec/nthroot.m, mpfr_function_d.cc: Support for vectorization of n

https://sourceforge.net/p/octave/interval/ci/e841c37383d6
summary:     doc: NEWS.texinfo: Summarized recent changes from GSoC

https://sourceforge.net/p/octave/interval/ci/bc3e9523d42a
summary:     mpfr_vector_dot_d.cc: Fixed bug when broadcasting with one size equal to zero

https://sourceforge.net/p/octave/interval/ci/0160ff2c2134
summary:     doc: examples.texinfo: Updated example for the latest Symbolic package version

https://sourceforge.net/p/octave/interval/ci/f37e1074c3ce
summary:     @infsup/*.m, @infsupdec/*.m: Added missing N-dimensional versions of tests

https://sourceforge.net/p/octave/interval/ci/eb6b9c01bf6a
summary:     @infsup/*.m, @infsupdec/*.m: N-dimensional versions of all ternary tests

summary:     @infsup/*.m, @infsupdec/*.m: N-dimensional versions of all binary tests

https://sourceforge.net/p/octave/interval/ci/9440d748b4aa
summary:     @infsup/*.m, @infsupdec*.m: N-dimensional version of all unary tests

https://sourceforge.net/p/octave/interval/ci/cb5b7b824400
summary:     @infsup/powrev2.m: Fixed bug when called with vector arguments

https://sourceforge.net/p/octave/interval/ci/a3e8fdb85b97
summary:     @infsup/pownrev.m, @infsupdec/pownrev.m: Reworked vectorization test

https://sourceforge.net/p/octave/interval/ci/aae143789b81
summary:     @infsup/pown.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/f05814665570
summary:     @infsup/nthroot.n: Reworked vectorization test

https://sourceforge.net/p/octave/interval/ci/c4c27574e331
summary:     @infsup/nthroot.m, @infsupdec/nthroot.m: Clarified that N must be scalar

https://sourceforge.net/p/octave/interval/ci/3bfd4e7e6600
summary:     @infup/pow.m, @infsupdec/pow.m: Fixed bug when called with vector arguments

https://sourceforge.net/p/octave/interval/ci/e013c55a4e2a
summary:     @infsup/overlap.m, @infsupdec/overlap.m: Fixed formatting of vector test.

summary:     doc: Modified a test so that it now passes

https://sourceforge.net/p/octave/interval/ci/2404686b07fb
summary:     doc: Fixed formatting of example

https://sourceforge.net/p/octave/interval/ci/c078807c6b32
summary:     doc: SKIP an example that always fail in the doc-test

https://sourceforge.net/p/octave/interval/ci/598dea2b8eb8
summary:     doc: Fixed missed ending of example in Getting Started

https://sourceforge.net/p/octave/interval/ci/66ca24edfba4
summary:     Updated coding style for all infsupdec-class functions

https://sourceforge.net/p/octave/interval/ci/4c629b9000fc
summary:     Updated coding style for all infsup-class functions

summary:     Updated coding style for all non-class functions

https://sourceforge.net/p/octave/interval/ci/f911db83df14
summary:     @infsupdec/dot.m: Fixed wrong size of decoration when called with two empty matrices

https://sourceforge.net/p/octave/interval/ci/8e6a9e37ee40
summary:     Small updates to documentation and comments for a lot of function to account for the support of N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/a2ef249d54ae
summary:     doc: A small update to Examples, the interval Newton method can only find zeros inside the initial interval

https://sourceforge.net/p/octave/interval/ci/76431745772f
summary:     doc: Updates to Getting Started, mainly how to create N-dimensional arrays

summary:     doc: Small updates to Preface regarding N-dimensional arrays and fixed one link

summary:     ctc_intersect.m, ctc_union.m: Fixed bugs when used for vectorization and when called with 0 or 1 output arguments

https://sourceforge.net/p/octave/interval/ci/0e01dc19dc75
summary:     @infsup/sumsq.m: Updated to use the new functionality of dot.m

https://sourceforge.net/p/octave/interval/ci/2fba2056ed31
summary:     @infsup/dot.m, @infsupdec/dot.m, mpfr_vector_dot_d.cc: Added support for N-dimensional vectors. Moved all vectorization to the oct-file. Small changes to functionality to mimic how the sum function works.

https://sourceforge.net/p/octave/interval/ci/05fc90112ea9
summary:     ctc_intersect.m, ctc_union.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/2a4d1e9fa43e
summary:     @infsup/fsolve.m: Added support for N-dimensional arrays. Fixed problem with the function in the example. Improved performance when creating the cell used in vectorization.

https://sourceforge.net/p/octave/interval/ci/7a96c346225a
summary:     @infsup/disp.m: Fixed wrong enumeration of submatrices

https://sourceforge.net/p/octave/interval/ci/933117890e45
summary:     Fixed typo in NEWS.texinfo

summary:     @infsup/diag.m: Added description of the previous bug fix in the NEWS file

https://sourceforge.net/p/octave/interval/ci/07ebc81867cd
summary:     @infsup/diag.m: Fixed error when called with more than 1 argument

https://sourceforge.net/p/octave/interval/ci/554c34fb3246
summary:     @infsup/meshgrid.m, @infsupdec/meshgrid.m: Removed these functions, now falls back on standard implementation, also updated index

summary:     @infsup/plot.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/5e0100cdc25f
summary:     @infsup/plot3.m: Small change to allow for N-dimensional arrays as input

https://sourceforge.net/p/octave/interval/ci/695a223ccbb7
summary:     @infsupdec/prod.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/55f570c13f01
summary:     @infsup/prod.m: Added support for N-dimensional arrays. Removed short circuit in simple cases.

https://sourceforge.net/p/octave/interval/ci/445d7e5150aa
summary:     @infsup/sum.m, @infsupdec/sum.m, mpfr_vector_sum_d.cc: Added support for N-dimensional vectors. Moved all vectorization to the oct-file. Small changes to functionality to mimic Octaves standard sum function.

summary:     @infsup/fminsearch.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/19724b3f581e
summary:     mpfr_function_d.cc: Finalized support for N-dimensional arrays with binary functions and added support for it with ternary functions.

https://sourceforge.net/p/octave/interval/ci/023e2788e445

https://sourceforge.net/p/octave/interval/ci/617484113019
summary:     @infsupdec/infsupdec.m: Added full support for creating N-dimensional arrays and added tests

https://sourceforge.net/p/octave/interval/ci/f1c28a3d48b7
summary:     @infsup/subset.m, @infsupdec/subset.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/8e9990e3ed13
summary:     @infsup/strictsubset.m, @infsupdec/strictsubset.m: Fixed coding style and updated documentation

https://sourceforge.net/p/octave/interval/ci/5ca9b26b580d
summary:     @infsup/strictprecedes.m, @infsupdec/strictprecedes.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/1c76721003f1
summary:     @infsup/sdist.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/e044dea7c820
summary:     @infsup/precedes.m, @infsupdec/precedes.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/d48ca9206299
summary:     @infsup/overlap.m, @infsupdec/overlap.m: Fixed coding style and updated documentation

summary:     @infsup/issingleton.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/11dccbcfd97e
summary:     @infsup/ismember.m: Updated documentation

summary:     @infsup/isentire.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/9fc874b71533
summary:     @infsup/isempty.m, @infsupdec/isempty.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/0c5c4eaf263b
summary:     @infsup/iscommoninterval.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/c980da83b634
summary:     @infsup/interior.m, @infsupdec/interior.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/f67d78651aee
summary:     @infsup/idist.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/70a389bd59f5
summary:     @infsup/hdist.m: Fixed coding style and updated documentation.

https://sourceforge.net/p/octave/interval/ci/e34401d2a6d6
summary:     @infsup/sin.m, @infsupdec/sin.m: Added workaround for bug #51283

https://sourceforge.net/p/octave/interval/ci/fea4c6516101
summary:     @infsup/gt.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/ce973432d240
summary:     @infsup/ge.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/a435986d4b0c
summary:     @infsup/lt.m, @infsupdec/lt.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/23fa89e3c461
summary:     @infsup/le.m, @infsupdec/le.m: Updated documentation

https://sourceforge.net/p/octave/interval/ci/730397b9e339
summary:     @infsup/disjoint.m, @infsupdec/disjoint.m: Updated documentation

summary:     mpfr_function_d.cc: Added support for N-dimensional arrays for unary functions. Also temporary support for binary functions.

https://sourceforge.net/p/octave/interval/ci/9661e0256c51
summary:     crlibm_function.cc: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/0ec3d29c0779
summary:     @infsup/infsup.m: Fixed documentation and added missing line continuation

https://sourceforge.net/p/octave/interval/ci/ebb00e763a46
summary:     @infsup/disp.m: Fixed documentation

https://sourceforge.net/p/octave/interval/ci/1ab2749da374
summary:     @infsup/size.m: Fixed documentation

https://sourceforge.net/p/octave/interval/ci/507ca8478b72
summary:     @infsup/size.m: Fixes to the documentation

https://sourceforge.net/p/octave/interval/ci/667da39afece
summary:     nai.m: Small fix to one of the tests

https://sourceforge.net/p/octave/interval/ci/3eb13f91065a
summary:     hull.m: Fixes according to Olivers review

https://sourceforge.net/p/octave/interval/ci/20d784a6605d
summary:     @infsup/display.m: Vectorized loop

https://sourceforge.net/p/octave/interval/ci/fab7aa26410f
summary:     @infsup/disp.m: Fixes according to Olivers review, mainly details in the output

https://sourceforge.net/p/octave/interval/ci/4959566545db
summary:     @infsup/infsup.m: Updated documentation and added test for N-dimensional arrays

summary:     @infsup/infsup.m: Fixed coding style

https://sourceforge.net/p/octave/interval/ci/486a73046d5e
summary:     @infsup/disp.m: Updated documentation and added more tests for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/b6a0435da31f
summary:     exacttointerval.m: Uppdated documentation and added tests for N-dimensional arrays

summary:     @infsup/intervaltotext.m, @infsupdec/intervaltotext.m: Updated documentation and added tests for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/55a3e708aef4
summary:     @infsup/intervaltotext.m: Fixed coding style

https://sourceforge.net/p/octave/interval/ci/536b0c5023ee
summary:     @infsup/subsref.m, @infsupdec/subsref.m: Added tests for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/c9654a4dcd8d
summary:     @infsup/size.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/636709d06194
summary:     @infsup/end.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/e6451e037120
summary:     nai.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/837cccf1d627
summary:     @infsup/resize.m, @infsupdec/resize.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/11cb9006b9ea
summary:     @infsup/reshape.m, @infsupdec/reshape.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/c276af6c42ae

https://sourceforge.net/p/octave/interval/ci/c92c829dc946

https://sourceforge.net/p/octave/interval/ci/0f8f6864123e
summary:     @infsup/meshgrid.m, @infsupdec/meshgrid.m: Added support for outputting 3-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/b7faafef6030
summary:     @infsup/cat.m, @infsupdec/cat.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/d416a17a6d3d
summary:     hull.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/031831f6bdfd
summary:     empty.m, entire.m: Added support for N-dimensional arrays

https://sourceforge.net/p/octave/interval/ci/2816b68b83c4
summary:     @infsup/display.m: Added support for displaying high dimensional arrays

https://sourceforge.net/p/octave/interval/ci/79c8dfa8ac54
summary:     @infsup/disp.m: Added support for displaying high dimensional arrays

https://sourceforge.net/p/octave/interval/ci/cc87924e52ac
summary:     @infsup/disp.m: Fixed coding style

https://sourceforge.net/p/octave/interval/ci/29a6b4ecda2a
summary:     @infsupdec/infsupdec.m: Temporary fix for creating high dimensional arrays

https://sourceforge.net/p/octave/interval/ci/039abcf7623d
summary:     @infsupdec/infsupdec.m: Fixed coding style

## Friday, 18 August 2017

### The Final Polish

We are preparing for releasing version 3.0.0 of the interval package and this last week have mainly been about fixing minor bugs related to the release. I mention two of the more interesting bugs here.

### Compact Format

We (Oliver) recently added support for "format compact" when printing intervals. It turns out that the way to determine if compact format is enabled differs very much between different version of Octave. There are at least three different ways to get the information.

In the older releases (< 4.2.0 I believe) you use "get (0, "FormatSpacing")" but there appear to be a bug for version < 4.0.0 for which this always return "loose".

For the current tip of the development branch you can use "[~, spacing] = format ()" to get the spacing.

Finally in between these two version you use "__compactformat__ ()".

In the end Oliver, probably, found a way to handle this mess and compact format should now be fully supported for intervals. The function to do this is available here https://sourceforge.net/p/octave/interval/ci/default/tree/inst/@infsup/private/__loosespacing__.m.

### Dot-product of Empty Matrices

When updating "dot" to support N-dimensional arrays I also modified it so that it behaves similar to Octaves standard implementation. The difference is in how it handles empty input. Previously we had

> x = infsupdec (ones (0, 2));
> dot (x, x)
ans = 0×2 interval matrix

but with the new version we get

> dot (x, x)
ans = 1×2 interval vector
[0]_com   [0]_com

which is consistent with the standard implementation.

In the function we use "min" to compute the decoration for the result. Normally "min (x)" and "dot (x, x)" returns results of the same size (the dimension along which it is computed is set to 1), but they handle empty input differently. We have

> x = ones (0, 2);
> dot (x, x)
ans =
0   0
> min (x)
ans = [](0x2)

This meant that the decoration would be incorrect since the implementation assumed they always had the same size. Fortunately the solution was very simple. If the dimension along which we are computing the dot-product is zero. the decoration should always be "com". So just adding a check for that was enough.

You could argue that "min (ones (0, 2))" should return "[inf, inf]" similarly to how many of the other reductions, like "sum" or "prod", return their unit for empty input. But this would most likely be very confusing for a lot of people. And it is not compatible with how Matlab does it either.

## Updates on the Taylor Package

I have also had some time to work on the Taylor package this week. The basic utility functions are now completed and I have started to work on functions for actually computing with Taylor expansions. At the moment there are only a limited amount of functions implemented. For example we can calculate the Taylor expansion of order 4 for the functions $\frac{e^x + \log(x)}{1 + x}$ at $x = 5$.

## Create a variable of degree 4 and with value 5
> x = taylor (infsupdec (5), 4)
x = [5]_com + [1]_com X + [0]_com X^2 + [0]_com X^3 + [0]_com X^4

## Calculate the function
> (exp (x) + log (x))./(1 + x)
ans = [25.003, 25.004]_com + [20.601, 20.602]_com X + [8.9308, 8.9309]_com X^2 + [2.6345, 2.6346]_com X^3 + [0.59148, 0.59149]_com X^4

## Friday, 11 August 2017

### Improving the Automatic Tests

Oliver and I have been working on improving the test framework used for the interval package. The package shares a large number of tests with other interval packages through an interval test framework that Oliver created. Here is the repository.

## Creating the Tests

Previously these tests were separated from the rest of the package and you usually ran them with help of the Makefile. Now Oliver has moved them to the m-files and you can run them, together with the other tests for the function, with test @infsup/function in Octave. This makes it much easier to test the functions directly.

In addition to making the tests easier to use we also wanted to extend them to not only test scalar evaluation but also vector evaluation. The test data, input ad expected output, is stored in a cell array and when performing the scalar testing we simply loop over that cell and run the function for each element. The actual code looks like this (in this case for plus)

%!test
%! # Scalar evaluation
%! for testcase = [testcases]'
%!   assert (isequaln (...
%!     plus (testcase.in{1}, testcase.in{2}), ...
%!     testcase.out));
%! endfor

For testing the vector evaluation we simply concatenate the cell array into a vector and give that to the function. Here is what that code looks like

%! # Vector evaluation
%! in1 = vertcat (vertcat (testcases.in){:, 1});
%! in2 = vertcat (vertcat (testcases.in){:, 2});
%! out = vertcat (testcases.out);
%! assert (isequaln (plus (in1, in2), out));

Lastly we also wanted to test evaluation of N-dimensional arrays. This is done by concatenating the data into a vector and then reshape that vector into an N-dimensional array. But what size should we use for the array? Well, we want to have at least three dimensions because otherwise we are not really testing N-dimensional arrays. My solution was to completely factor the length of the vector and use that as size, testsize = factor (length (in1)), and if the length of the vector has two or fewer factors we add a few elements to the end until we get at least three factors. This is the code for that

%!test
%! # N-dimensional array evaluation
%! in1 = vertcat (vertcat (testcases.in){:, 1});
%! in2 = vertcat (vertcat (testcases.in){:, 2});
%! out = vertcat (testcases.out);
%! # Reshape data
%! i = -1;
%! do
%!   i = i + 1;
%!   testsize = factor (numel (in1) + i);
%! until (numel (testsize) > 2)
%! in1 = reshape ([in1; in1(1:i)], testsize);
%! in2 = reshape ([in2; in2(1:i)], testsize);
%! out = reshape ([out; out(1:i)], testsize);
%! assert (isequaln (plus (in1, in2), out));

This works very well, except when the number of test cases is to small. If the number of test is less than four this will fail. But there are only a handful of functions with that few tests so I fixed those independently.

## Running the tests

Okay, so we have created a bunch of new tests for the package. Do we actually find any new bugs with them? Yes!

The function pow.m failed on the vector test. The problem? In one place $\&\&$ was used instead of $\&$. For scalar input I believe these behave the same but they differ for vector input.

Both the function nthroot.m and the function pownrev.m failed the vector test. Neither allowed vectorization of the integer parameter. For nthroot.m this is the same for standard Octave version so it should perhaps not be treated as a bug. The function pownrev.m uses nthroot.m internally so it also had the same limitation. This time I would however treat it as a bug because the function pown.m does allow vectorization of the integer parameter and if that supports it the reverse function should probably also do it. So I implemented support for vectorization of the integer parameter for both nthroot.m and pownrev.m and they now pass the test.

No problems were found with the N-dimensional tests that the vector tests did not find. This is a good indication that the support for N-dimensional arrays is at least partly correct. Always good to know!