Provide a more efficient drawing routine for constrained image sets

[lgeretti]

The routine for plotting nonlinear constrained image sets is quite inefficient, which is particularly noticeable when using differential inclusions. We should provide an alternative implementation that sacrifices visualization in favor of a more affordable computation time. We should still allow configuration of the preferred visualization method.

[lgeretti]

Adding a remark to say that the Lorenz system is a good candidate for evaluation of efficient drawing: as it appears, it alternates sections where the step size is small with sections where it is large. An affine draw with just one allowed splitting is really slow, around 12 minutes, while box draw is below 1 second. The graphical results (affine followed by box, for the first two variables):

[lgeretti]

I will work on improving approach 3 according to your suggestions.

I noticed a strange splitting behavior on the harmonic_inclusion example, as can be seen from this figure:

Clearly the bottleneck lies in the thick splitting on the left side. The image was obtained by using the default example file where evolution time is instead 6/4_q. Please note that up to 5/4_q the plotting is rather quick: its performance drops dramatically after that. The noises are 1/16 and 1/32 for x[0] and x[1] respectively.

If on the other hand we use 1/16 for both, we obtain

The first case writes the figure in 22 seconds, while the second case in 5 seconds even if the noise is larger.

Do you have any clues as to the reason of this behavior?

[lgeretti]

I add Pieter comments here:

Our constrained image sets have the form {f(s) | s in D | g(s) in C} where f,g are nonlinear / polynomial and C,D are boxes. D is bounded but may have high dimension. For plotting, we can assume R^2 is the codomain of f.

We currently have three main approaches: Let [f] denotes an interval/box over-approximation to f, and <f|D> an affine over-approximation to f on D.

1. Compute a grid-based over-approximation, by subdividing f into boxes B', and checking whether B' intersects the set using constraint programming.

2. Compute a box-based over-approximation. We can most easily do this by subdividing D into smaller boxes D', and plotting boxed f whenever g intersects C.

3. Compute a polytope-based over-approximation. We do this by again subdividing D into smaller boxes D', and plotting the polyhedra { <f|D'>(s) | s in D' | <g|D'>(s) in C } by finding their boundaries.

   An alternative approach is:

4. Try to find the boundary of the set directly, by following the image of the boundary of the parameter domain, and the image of the critical set on which df/ds is singular.

1) involves nonlinear programming, which is hard. 2) involves subdividing D, which is inefficient in high dimensions. A further disadvantage of 1) and 2) is that boxes give bad pictures. 4) should be quite efficient, but may be hard to implement due to the need to find critical points of f and handle singularities on the boundary. I think 3) is a promising strategy but needs work.

Currently, in 3), we subdivide D to a (user-)given fixed depth. However, it should be possible to adaptively subdivide. In principle, we need only subdivide a variable if f,g are not close enough to linear in that variable. For f, 'close enough' should be defined in terms of the resolution and scaling (change in variable x/y per pixel). If the affine over-approximation introduces an error of at most one pixel, that's definitely close enough. For g it's more complicated, but the allowable error in g can be estimated in terms of the derivative of f.