Re: How to handle reductions on interval arrays

[Oliver Heimlich]

On 28.06.2017 16:57, Joel Dahne wrote:
> Hi Oliver,
> 
> A new week and some new questions!
> 
> I have almost finished the work on vectorization. The only thing that is
> left is how to handle ctc_union.m, ctc_intersection.m and
> fsolve.m. These functions are a bit more complicated however so I will
> wait a bit with them.

These are the algorithmically most complex functions in the package and
you will need background knowledge in contractor programming.

https://www.youtube.com/watch?v=kxYh2cXHdNQ

You should definitely put them in the end of your work queue.

> Instead I have started to work on folding/reduction functions (Octave
> seems to use the word reduction in the source code at least). I have
> implemented support for N-dimensional arrays in sum.m (and also in
> mpfr_vector_sum_d.cc).

The term “reduction operation” is also used by the IEEE standards on
numeric computing (IEEE 754, IEEE 1788).  So, I find this the most
appropriate term.

> In that change there are also some other changes,
> I have moved all of the vectorization inside the oct-file (yielding
> about a 2 times speed up) and have slightly modified how it handles empty
> arrays to be consistent with how the standard sum-function does it.

Nice, moving code from m-files to oct-files often results in a
performance improvement.

> When I were to work on dot.m I realized that Octaves standard
> implementation is a bit odd. It is inconsistent with both the
> sum-function and Matlab. See bug #51333,
> https://savannah.gnu.org/bugs/?51333, on Savannah for more detailed
> info. I think that the interval package should follow Octave, but in
> this case Octave is weird. Probably we should wait and see what happens
> with my bug report.

Yes, let's see how this is supposed to work in Octave and then mimic the
behavior for intervals.

> I have also worked on prod.m and while doing that I encountered
> something which I think is a bug. But I'm not sure so I thought I would
> ask first. A minimal example is the code
> 
>> x = infsup (1:4);
>> i = logical ([0 1 0 0]);
>> x(i)
> ans = [2]
> 
> where everything is as it should. However when I try to put this code
> inside an interval function, for example testme.m, it does not
> work. With
> 
> function testme (arg)
>   x = infsup (1:4);
>   i = logical ([0 1 0 0]);
>   x(i);
> endfunction
> 
> we get
> 
>> testme (infsup ([]))
> error: x(2): out of bound 1
> error: called from
>     testme at line 7 column 3
> 
> It is important that testme.m is inside the @infsup folder, otherwise it
> seems to work. Do you have any idea of what this is?

Yes, it's a “weird” thing with classes in Octave that you should know
and remember when programming in methods.  Any objects (in the sense of
OOP) in Octave are stored as structure arrays [1] internally.  For the
interval package, every interval object is a 1x1 structure array, so you
may think of a simple structure for each interval object.

The bare interval (@infsup) structure consists of the fields:
  inf = ND array of lower boundary
  sup = ND array of upper boundary

The decorated interval (@infsupdec) structure consists of the fields:
  infsup = parent class (@infsup) object
  dec = ND array of decoration information

For example, the internal representation of a 5x5 interval matrix is a
scalar structure, where field 'inf' is a 5x5 double matrix and field
'sup' is a 5x5 double matrix.

Whenever you do an indexing expression / field access on an object using
the operator syntax, e. g.,

 - x(2)
 - x.inf

The behavior depends on whether you do it inside a method of the class
(inst/@infsup/*.m) or outside a method of the class (inst/*.m or
anywhere else).  Inside a method you operate directly on the internal
representation (on the structure) and since the structure is scalar, you
cannot access x(2).

Outside of the class folder, the expression is resolved by the
interpreter by calling the overridden “subsref” or “subsasgn” methods of
the class.

There is a reason for this: First, you don't want to create an infinite
loop when you override the subsasgn/subsref methods.  Second, you might
want to completely encapsulate and hide the fields of the object, but
your methods require direct access to the internal storage.

For example:
 - x.inf (on the command line, outside of the class methods)
 - calls @infsup/subsref with idx.type = '.'
 - calls @infsup/inf
 - evaluates x.inf (now inside a method of the class)
 - returns the inf field of the object structure

As a solution to your original problem: Inside a method of the class you
can still use subsref/subsasgn to do indexing expressions on the object
(and not on the internal structure).  The substruct function is a
helpful tool for that.

Oliver


[1]
https://www.gnu.org/software/octave/doc/interpreter/Structure-Arrays.html