How to handle reductions on interval arrays

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

How to handle reductions on interval arrays

Urathai
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.

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). 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.

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.

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?

Best,
Joel Dahne

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to handle reductions on interval arrays

Urathai

Oliver Heimlich writes:

> 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.

If I recall we took the same MOOC about interval analysis some year
ago. But yes, I will have to go back and repeat the theory behind SIVIA.

>> 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.

Then I will wait for that!

> 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.

Okay, that makes sense. In the end I moved away from using such
indexing. I will write a little about that on my blog tomorrow.

Best,
Joel

Loading...