Interval Picard continues

Last week Jan continued the implementation of an interval Picard-based ODE IVP solver using the framework for safe numerics developed with Veronica. The design follows the program (described in previous posts) of keeping things generic in the numeric base type and function representation to facilitate the comparison and evaluation of various alternatives. He will continue this work over the next weeks.

Floating sequences

Continuing their work on generic differential equation solving in Scala Veronica and Jan have made sequences of Floating values an instance of Floating. This makes it possible to use algorithms for solving linear ODEs to be applied to higher order ODEs by translating them to systems of linear ODEs. When the original solver algorithm is written polymorphically over the real number type the extension comes at no extra cost as the Floating instance to be used is obtained from the code for the field and initial condition of the system of ODEs. Coming work aims at completing the supported set of Interval Floating methods and adding Floating instances for Floating-valued functions to the framework.

Reducing noise in failing test cases

During the past week Adam has looked at test case minimization, also known as shrinking, and how to apply it to test cases generated from external specifications. Shrinking aims to make debugging easier by removing noise, that is information which is not relevant to the error which the a failing test case provokes. Given a failing test case and the specification that was used to generate it, Adam investigated ways to obtain shrunk testcases which are still valid with respect to the specification. This included a simple approach, where the specification was only used to validate shrunk test case variants, and a more sophisticated method where the specification was used to guide the shrinking. In both cases the number of cases to consider is very large, as an optimal solution really requires that all possible shrunk variants are enumerated. In the latter case, however, the task is reduced since variants which are invalid with repsect to the specification can be omitted. 

Sampling the strings of a regular language

This week Adam focused on the generation of strings which match a specific regular expression. Several ways of approaching this has problem have been described in the literature. For example, two methods are given in the functional pearl Enumerating the strings of regular languages by M. Douglas McIlroy, but focus has been on the enumeration of all such strings, an approach that does not scale well. A practical alternative is to give up on the constraint that each string should be generated with equal probability. An approach taken by the tool Xeger is to generate the strings by transforming the regular expression into a non-deterministic finite automaton, and then to generate the string simulating the automaton by giving each option in a non-deterministic choice equal probability. This corresponds to a probabilistic automaton whose stochastic matrix consists of uniform distributions.

Generating generators

This week Adam started looking at QuickCheck-style property-based testing. He has been focusing on understanding how to create generators "from scratch" as well as how to use combinators to create new generators from existing ones.

One problem which Adam faced was how to create a generator for a nested data structure. Perhaps unsurprisingly, the solution was to use recursion in the comprehension which constructs the generator.

Parametric IVP solver

Building on lazy analytic functions Veronica and Jan have now completed an implementation of a generic solver for initial value problems. The ideas for the solver are borrowed from the work on validated numerics by Davis and Tucker. The solver computes an interval enclosure of the solution and takes the following parameters:
a)  the field as a function of type (A,A) => A for any type A that can be used for endpoints of intervals and as a representation of the real numbers,
b) one such type
c) the initial value of type A
d) the degree for the Taylor expansion used to approximate the solution at each step
e) the step size
Our implementations are in Scala where we could make use parametric polymorphism, overloading and a lazy datastructure in the procedure for automatic differentiation. As a result the taylor coefficients for the sought function needed at the points where the solution is calculated and at intervals are expressed as a mutually recursive definition using the field and the initial value in the lazy data structure. We have performed some experiments with both constructive reals (a Java implementation by Boehm) and BigDecimals (a Java library). 

The following graph summarizes one experiment with the IVP given by x'(t) = -x(t)^3 with x(0)=1up to t = 4. In the graph below we show the execution time and the width of the enclosing interval at time = 4 for  constructive reals and big decimals with precision 10. 

Future work includes a) completing the library for using big decimals as endpoints, b) completing the library for intervals and c) using these libraries in an implementation of a Picard solver. We will also start writing about the program, the experiments and how to use the solver in the context of the parametric version of Acumen.

Interval endpoints should be rounded

Continuing their work on parametrised safe numerics in Scala, Veronica and Jan now have an implementation of intervals parameterized on the type for the endpoints. This type is required  to have numeric operations in 2 versions: for rounding up  and for rounding down. These operations can then be used to define interval extensions of numeric operations that are safe, in the sense that they enclose the result.

Consider the example of addition of intervals [a,b] + [c,d] = [a+c,b+d] when a,b,c and d belong to a type A in which the sum of two values is not necessarily in A. To make the resulting operation safe we need to use a lower bound plusL(a,c) on a+c and an upper bound plusU(b,d) on b+d as the endpoints of the result. Making this a special case of a generic scheme we define plus[Interval[A]]([a,b],[c,d]) as [plusL[A](a,c),plus[A]U(b,d)].

Veronica and Jan have tested the scheme with  BigDecimal as an  instance of the the type of endpoints. The size of the mantissa can be provided as an argument to the instance.  Coming work will focus on a) adding to the set of operations supported by the BigDecimal instance, b) completing the numeric operations for intervals (trigonometric functions) and c) implementing a validated IVP solver using Taylor expansions and intervals of constructive reals.