Rigorous Function Calculi
Theory, Practice and Problems
Pieter Collins
Department of Advanced Computing Sciences
Maastricht University

pieter.collins@maastrichtuniversity.nl
Computing with Infinite Data (CID)
Kyoto, 13 Febrary 2023
This project has received funding from the European Union’s Horizon 2020 research andinnovation programme under the Marie Skłodowska-Curie grant agreement No 731143.
Kyoto, 13 February 2023

Rigorous Function Calculi – 1 / 48

Introduction• Outline
Motivation
TheoryPractice
Problems
Conclusions

Introduction

Kyoto, 13 February 2023

Rigorous Function Calculi – 2 / 48

Outline
Introduction• Outline
Motivation
TheoryPractice
Problems
Conclusions

Kyoto, 13 February 2023

• Motivation — Dynamic Systems• Theory — Computable Analysis• Practice — Rigorous Numerics• Problems — To work on together!
• Concluding Remarks

Rigorous Function Calculi – 3 / 48

Introduction
Motivation• Function calculus• Applications• Hybrid systems• Function calculus in CID• Other tools• A quick look
TheoryPractice

Motivation — Dynamic Systems

Problems
Conclusions

Kyoto, 13 February 2023

Rigorous Function Calculi – 4 / 48

What is function calculus?
A function calculus is a computational toolbox for working with functions on
continuous spaces.
Many different function spaces: continuous, smooth, measurable, Sobolev,
piecewise, etc.
We would like to be able to work with functions in a natural, rigorous, and efficientway!
Rigour is especially important in mathematical proofs, verification of safety-criticalsystems, and long chains of reasoning.

Kyoto, 13 February 2023

Rigorous Function Calculi – 5 / 48

Why function calculus?
Many problems in applied mathematics are formulated in terms of functions:• Trajectories and flow tubes of ordinary differential equations.
• Trajectory sets of trajectories for differential inclusions.
• Probability densities of stochastic systems.

• State spaces and solutions of partial differential equations.
• Feedback for control systems.

• Crossing times for hybrid systems.
• Parametrised families of solutions.
• Reachable and safe sets.

Kyoto, 13 February 2023

Rigorous Function Calculi – 6 / 48

Verification of hybrid systems
A problem requiring many different function calculus operations is verification of
hybrid systems. A hybrid system is a dynamic system in which
the continuous state in location qi follows a differential equation ẋ = fi (x),until some guard condition gij (x) > 0 becomes true,
at which time the state jumps to x′ = rij (x) in location qj .

The safety verification problem is to determine whether, given a set X0 of initialstates, and a set S of safe states, the evolution starting in X0 remains in S .
Kyoto, 13 February 2023

Rigorous Function Calculi – 7 / 48

Function calculus for hybrid systems
For application to verification of hybrid dynamic systems, our function calculus
should therefore at the least support:• Solution of ordinary differential equations φ̇(x, t) = f (φ(x, t)).

• Solution of single-equation implicit function problems g(φ(x, τ (x))) = 0.

• Composition of functions r ◦ φ ◦ (id x , τ )?

• Constrained feasibility ψ(D) ∩ C = ∅.
• Over approximation of the range f (D).
These operations must be implemented rigorously and efficiently.

Kyoto, 13 February 2023

Rigorous Function Calculi – 8 / 48

Function Calculus in CID
Frameworks for function calculi in EU project “Computing with Infinite Data”:• A RIADNE (Collins, Geretti, Villa et al.) for verification of hybrid systems (C++).

http://www.ariadne-cps.org/• AERN tool (Konec̆ný et al.) for effective real computation (Haskell).
http://michalkonecny.github.io/aern/ site/• iRRAM package (Müller, Brauße) for real number arithmetic (C++).
http://irram.uni-trier.de/• ERC language (Ziegler, Park et al.) for exact real computation.

Kyoto, 13 February 2023

Rigorous Function Calculi – 9 / 48

Other Function Calculus Approaches
Other tools for rigorous numerics:• Cosy Infinity (Berz & Makino)• CAPD Library (Mrozek, Zgliczynski, Wilczak et al.)• Flow* (Chen, Abraham et al.)• VNode (Nedialkov et al.)• CORA (Althoff et al.)

• JuliaReach (Benet et al.)• Ellipsoidal calculus (Kurzhanski & Valyi)• Set-valued viability (Cardaliaguet, Quincampoix & Saint-Pierre)• Taylor-ellipsoid models (Houska, Villanueva & Chachuat)• AWA (Lohner)• DynIBEX (Chabert, Jaulin; Chapoutot,Alexandre dit Sandretto et al.)

Kyoto, 13 February 2023

Rigorous Function Calculi – 10 / 48

A quick look at A RIADNE (in C++)////F i l e c o m p u t e a r e a l . cpp
Build using :c l a n g ++c o m p u t e a r e a l . cpp − l a r i a d n e−ocomput# include < ariadne / ariadne . hpp >using namespace Ariadne ;# define PRINT ( expr ) { std :: cout < <# expr < < " : " < <( expr ) < < " \ n " ; }int main () {auto r = 6* atan (1/ sqrt (3 _q ));////Define a
The ’ q ’r e a l number .
c o n v e r t s t o an

AriadneRational

PRINT ( r );PRINT ( r . compute ( Accuracy (123 _bits )));//

Computewitha m axi m um

error

of

1/2ˆ123

PRINT ( r . compute ( Effort (123)));//

Compute e . g .

using

123

bits

of

precision .

PRINT ( r . compute ( Effort (123)). get ( precision (75)));//Compute ,andreturnwith

lessprecision}
Kyoto, 13 February 2023

Rigorous Function Calculi – 11 / 48

A quick look at A RIADNE (in Python)#Filec o m p u t e a r e a l . py

from ariadne import *if __name__ == ’ __main__ ’:r = 6* atan (1/ sqrt (3))##

D e f i n e a r e a l number .
sqrt ( . . . ) converts toanAriadneReal

print ( r )print ( r . compute ( Accuracy ( bips =123)))#

Computewitha m axi m um

error

of

1/2ˆ123

print ( r . compute ( Effort (123)))#

Compute e . g .

using

123

bits

of

precision .

print ( r . compute ( Effort (123)). get ( precision (75)))#
Kyoto, 13 February 2023

Compute ,andreturnwith

lessprecision
Rigorous Function Calculi – 12 / 48

Introduction
Motivation
Theory• Effective information• Validated information• Concrete objects• Computable analysis• Computable functions• Continuous functions• Discontinuous functions• Measurable functions

Theory — Computable Analysis

Practice
Problems
Conclusions

Kyoto, 13 February 2023

Rigorous Function Calculi – 13 / 48

Effective and symbolic information
Objects from uncountable spaces need an infinite amount of data to represent
exactly.
e.g. Real numbers can be specified by their decimal expansion,
such as π = 3.14159 · · · .
There may be more natural descriptions of an object, but they can all be encoded
as a sequence over some alphabet Σ.

Kyoto, 13 February 2023

Rigorous Function Calculi – 14 / 48

Effective and symbolic information
Objects from uncountable spaces need an infinite amount of data to represent
exactly.
e.g. Real numbers can be specified by their decimal expansion,
such as π = 3.14159 · · · .
There may be more natural descriptions of an object, but they can all be encoded
as a sequence over some alphabet Σ.
It’s slightly problematic to specify an infinite amount of information in practice...

Kyoto, 13 February 2023

Rigorous Function Calculi – 14 / 48

Effective and symbolic information
Objects from uncountable spaces need an infinite amount of data to represent
exactly.
e.g. Real numbers can be specified by their decimal expansion,
such as π = 3.14159 · · · .
There may be more natural descriptions of an object, but they can all be encoded
as a sequence over some alphabet Σ.
It’s slightly problematic to specify an infinite amount of information in practice...
An object is computable if it is possible to compute a complete description from afinite amount of information.
e.g. π = 4 limn→∞

Kyoto, 13 February 2023

Pn

(−1)k
k=0 2k+1 = 4 atan(1).

Rigorous Function Calculi – 14 / 48

Validated and approximate information
Since we generally don’t want to wait forever for our computations to terminate, we
work with finite-precision approximations.
Classical numerical packages work with floating-point numbers and do not control
the errors.

• e.g. In double-precision, π ∼= 3.141592653589793.

A common approximation is π ≈ 22/7 ∼= 3.142857142857143.

Results of predicates on approximate objects could be wrong, so should not be
used!.
Instead, we provide an error bound, or lower and upper bounds.
e.g. Given π = 3.14159 · · · , we know π ∈ [3.14159 : 3.14160].
Such information is validated, and guaranteed to be correct.
Obtain rigorous results by working with validated objects throught.
(Approximate objects can still be useful for preconditioning rigorous algorithms.)
Kyoto, 13 February 2023

Rigorous Function Calculi – 15 / 48

Algebraic and concrete objects
Algbraic objects from countable spaces like Z, Q, and the dyadic numbers Q2 canbe described with a finite amout of data, support exact operations, and can bedecidably compared and tested for equality.
Working with these numbers directly is inefficient, as memory allocation isexpensive, and representations get large.
4 1• e.g. f (x) = 23
x(1
−
x)
,
f
( 2 ) = 29529021591697849/30814043149172736.
6

It is fastest to work with builtin objects of a fixed size, like double.
“Multiple-sized” objects, such as floating-point numbers provided by the MPFR
library can be allocated efficiently and have reasonable performance.
However, fixed-size types are finite and cannot support exact arithmetic.
Instead, arithmetic is rounded, either upwards, downwards or to the nearest
representable value.
Using outward-rounded arithmetic preserves bounds for quantities.
e.g. If x ∈ [x : x] and y ∈ [y : y] then x − y ∈ [x −↓ y : x −↑ y].
Kyoto, 13 February 2023

Rigorous Function Calculi – 16 / 48

Computable analysis
Encode objects from continuous spaces by streams of data.

• A representation of a space X using alphabet Σ is a partial surjective functionδ : Σω →7 X satisfying additional admissibility properties.

• A δ -name of x ∈ X is a sequence p such that δ(p) = x.

Representations δ1 , δ2 of X are equivalent if for any x ∈ X a δ2 -name of x canbe computed from any valid δ1 -name, and vice-versa.
A type X = (X, [δ]) is a space with an equivalence class of representations.
A representation δ of X induces the quotient topology on X , so any type also is a
topological space.
An operation f : X → Y is computable if there is a Turing-computable functionµ : Σω →7 Σω such that given any valid name p ∈ dom(δX ), we haveδY (µ(p)) = f (δX (p)).

Kyoto, 13 February 2023

Rigorous Function Calculi – 17 / 48

Computable functions
One of the fundamental results of computable analysis is that:
Any computable function f : X → Y is continuous.

Note that “computable” means that f can be effectively evaluated. It may still bepossible to give a complete description of f .

• e.g. Let h(x) = 0 for x < 0 and h(x) = 1 for x > 1. Then h is
uncomputable! What is the value of h(−1 + 3 × 0.3333 · · · )?

There is a standard way of encoding continuous functions f : X → Y such that:• The corresponding represenation γδX ,δY : Σω →7 C(X; Y ) is admissible,• f has a computable name if, and only if, it is a computable function, and• The evaluation map ε(f, x) : C(X; Y ) × X → Y taking (f, x) 7→ f (x) is
computable.

Hence given types X, Y, there is a natural type of continuous functions X → Y,
also denoted C(X; Y) or YX .
Kyoto, 13 February 2023

Rigorous Function Calculi – 18 / 48

Continuous function types
There are many important subtypes of (continuous) functions.
Functions f : X → Y are defined by:• Evaluation ε(f, x) = f (x).

Combining evaluations yields:• Composition [f ◦ g](x) = f (g(x)).

Continuous real functions on compact domains D have a natural uniform norm

• Norm ||f || ∞,D = supx∈D |f (x)|.
By the Weierstrass Approximation Theorem, can uniformly approximate acontinuous function on a bounded domain by a polynomial• Polynomial f˜D,ǫ = {f | kf − f˜k∞,D 6 ǫ.
Differentiation is uncomputable! Need to specify derivative information separately:• Differential f (6n) (x) = hf (x), f ′ (x), f ′′ (x), . . . , f (n) (x)i.
Often functions are defined by a formula involving elementary operations.

• Formula f (x) generated by 1, x, +, −, ×, ÷, |·|, exp, log, sin, cos, atan.

Kyoto, 13 February 2023

Rigorous Function Calculi – 19 / 48

Discontinuous function types
Multivalued functions take values in (e.g. compact) sets, defined by:• Evaluation F : X → K(Y).
Compact-valued functions are naturally upper-semicontinuous.
Piecewise-continuous functions f |Ai = fi based on topological partition
S{Ai | i ∈ I} for closed Ai with interior Ui = Ai \ j6=i Aj .

• Semicontinuous compact-set evaluation: F (x) = {fi (x) | x ∈ Ai } ∈ K(Y).

• Singleton evaluation in the interior of a continuity domain: f (x) = fi (x) if x ∈ Ui .

Integrable function spaces Lp defined as the effective (Cauchy) completion ofpolynomial/continuous/piecewise-continuous functions under the p-norms.

• Norms kf kp,D =

R

D



1/p
f (x)p dx.

Sobolev-spaces W k,p of functions whose k -th derivative is p-integrable defined interms of integrable functions and weak derivative or antiderivatives.
Measurable functions from (X, µ) to Y defined by:• Preimage f −1 : O(Y ) → M< (X, µ)Kyoto, 13 February 2023

Rigorous Function Calculi – 20 / 48

Measurable function typeFix a measure µ on X , defined as a valuation on open sets µ : O(X) → R+<.
Define the µ-lower-measurable sets as limits of sequences (Un ) satisfying∀m > n, µ(Um ∩ Un ) > µ(Un ) − 2−n .
The limit U∞ satisfies µ(Un \ U∞ ) 6 2−nIntersection, countable union and Cartesian product are computable operations.
The type of measurable functions from (X, µ) to (Y, τ ) is the type of maps
f −1 : O(Y ) → M< (X, µ) such that

V1 ⊂ V2 =⇒ f −1 (V1 ) ⊂ f −1 (V2 )Denote µ-measurable functions from X to Y by f : X µ
Clearly, if f : Y → Z is continuous, and g : X µcomposition f ◦ g : X µ Z can be computed.

Y.

Y is µ-measurable, then the

If f : X µ Y and g : X µ Z the product f × g : X µ Y × Z is computable
if X, Y are countably-based spaces, since (f × g)−1 (U × V ) = f −1 (U ) ∩ g −1 (V ).
Kyoto, 13 February 2023

Rigorous Function Calculi – 21 / 48

Introduction
Motivation
TheoryPractice• Numbers and logic• Function types• Function models• Algebraic equations• Differential equations• Function sets
• Hybrid systems

Practice — Rigorous Numerics

Problems
Conclusions

Kyoto, 13 February 2023

Rigorous Function Calculi – 22 / 48

Numbers and Logic
A RIADNE supports Integer, Dyadic and Rational classes:
w = Dyadic (5 ,3 u ) # w i s 5 / 2 ˆ 3
The Real number class can be defined by a symbolic formula:r = Real (6* atan (1/ sqrt (3)))
1
Bounds on a real number can be extracted with a given precision:x = Bounds [ FloatMP ]( r . get ( precision (128)))
Comparing real numbers is undecidable, so must return a Kleenean object:
k = Kleenean ( sin ( r ) >0)This can be checked using a given amount of computational “effort”:
v = k . check ( Effort (2))
The result can be converted to a boolean if reaquired:b = possibly ( v )Kyoto, 13 February 2023

Rigorous Function Calculi – 23 / 48

Function types
A RIADNE currently supports functions on Euclidean space, distinguishing Scalarand Vector arguments and results. Function types are tagged by the information
provided, which can be Effective, Validated or Approximate.
a = Dyadic ( " 1.875 " ); b = Decimal ( " 0.3 " )x = EffectiveV e ct o rM u l ti v ar i a te F un c ti o n . identity (2)h = EffectiveV e ct o rM u l ti v ar i a te F un c ti o n ([ a - x [0]* x [0] - b * x [1] , x [0]])Functions can be evaluated on (vectors of) concrete numbers:v = Vector [ FloatDP ]([ " 0.5 " ,1] , dp )h(v)evaluate (h , v )Automatic differentiation is used to compute the derivatives at a point:h . jacobian ( v )jacobian (h , v )h . differential (v ,3)

Kyoto, 13 February 2023

Rigorous Function Calculi – 24 / 48

Function models
Most concrete function approximations are only valid over bounded domains.
A powerful rigorous calculus for continuous functions Rn → R is based aroundthe Makino & Berz’s Taylor models: approximations kf − p ◦ s−1 k∞ 6 e. wheres is a scaling function, p is a polynomial and e is an error bound.
f (x)
s−1
e
p(z)
a

b

x

−1

0

+1 z

Standard functions operations, including evaluation, arithmetic, composition, andantidifferentiation, are available for Taylor models.
The efficiency of Taylor models relies on sweeping terms of p with small
coefficients into the error e.
Taylor models can be computed using a Taylor series with remainder term:p(z) =

Pn−1 1

Kyoto, 13 February 2023

k=0 k!

f (k) (0) z k ;e = rad1(n) ]([−1, +1]) .
[fn!
Rigorous Function Calculi – 25 / 48

Function models
Compute a polynomial approximation of h over a finite domain.
dom = BoxDomainType ([[0 ,1] ,[ dy_ (0.5) , dy_ (1.5)]])th = Validate d V e ct o r M ul t i v a ri a t e T ay l o r F un c t i on M o d e lD P (dom ,h , ThresholdSweeper ( dp ,1 e -4));Pre-compose this by the original h to obtain the second iterate.
thh = compose (h , th );Evaluate at the orignal x.
evaluate ( thh , v );h ( h ( v ));Compute the norm and range:
thh . range ()norm ( thh )Kyoto, 13 February 2023

Rigorous Function Calculi – 26 / 48

Algebraic equations
Solve a parametrised system of algebraic equations g(x, h(x)) = 0.
Requires an computational implementation of the implicit function theorem: asolution exists on a neighbourhood of x0 if D2 g(x0 , y0 ) is nonsingular.
A RIADNE has iterative solvers based on the interval Newton operatorĥ′ (x) = N [g, ĥ, h](x) = h(x) − D2 f (x, ĥ(x))−1 g(x, h(x))

and the Krawcykz operator isb
K[g, ĥ, h](x) = h(x) − M (x) g(x, h(x)) + (I − M (x)J(x))(ĥ(x) − h(x))
where Jb(x) = D2 f (x, ĥ(x)) and M (x) ≈ D2 f (x, h(x))−1 .

Kyoto, 13 February 2023

Rigorous Function Calculi – 27 / 48

Algebraic equations
Set g(x, y) = 4 + x − y 2 :x = RealVariable ( " x " ); y = RealVariable ( " y " )g = make_function ([ x , y ] ,4+ x - y * y )Look for solutions over [−1 : +1] with values also in [−1 : +1]:
xdom = BoxDomainType ([[ -1 ,+1]])
yrng = IntervalDomainType ([+1 ,+3])Compute the solution h to the scalar equation g(x, h(x)) = 0.
solver = IntervalNewtonSo lv er (1 e -8 ,6);h = solver . implicit (g , xdom , yrng );√Solution is y = 4 + x = 2 + x/4 − x2 /64 + · · · .

Kyoto, 13 February 2023

Rigorous Function Calculi – 28 / 48

Differential equations
The flow φ of ẋ = f (x, a) can be computed for t ∈ [0, h] using the Picard operator:

Rtφ(x, t, a) = x + 0 f (φ(x, τ, a), a)dτ.

Use an initial bound for the flow over a domain D , typically a box B such thatD + hf (B) ⊂ B.
Alternatively, compute Taylor expansion of φ(x, t), which satisfies the recurrence:∂φ|α|+k+1∂(f ◦ φ)|α|+k
=
α
k+1
∂x ∂t
∂xα ∂tk
Then taking cα,k =∂ |α|+k φ∂xα ∂tk x ,tc

φ̂(x, t) =

X

0

α

cα,k (x−xc ) (t−t0 )|α|6m∧k6n

Kyoto, 13 February 2023

and Cα,k =k±∂ |α|+k φ, we have∂xα tk B×[0,h]X|Cα,k −cα,k | |x−xc |α |t−t0 |k .

|α|6m∧k6n|α|=m∨k=n

Rigorous Function Calculi – 29 / 48

Differential equations
Consider the differential equation ẋ = f (x) = x for initial conditions in [−1 : +1]:
f = ValidatedV e ct o r Mu l ti v a ri a te F un c t io n . identity (1);dom = BoxDomainType ([[ -1 ,+1]]);For t ∈ [0, 1/2], the flow remains in the set [−3 : +3].
h =1/ two ;bbx = BoxDomainType ([[ -3 ,+3]]);
Compute the flow by evaluating Taylor series coefficients:
tolerance =1 e -8; order =12
integrator = TaylorSeriesInt eg ra to r ( tolerance , order );flow = integrator . flow_step (f , dom , h );Kyoto, 13 February 2023

Rigorous Function Calculi – 30 / 48

Differential equations
Define a more complicated problem:
x = RealVariable ( " x " )f = make_function ([ x ] ,[ - x + sin ( x )])dom = BoxDomainType ([ [1 ,2] ])h = Dyadic ( " 0.125 " )Solve using Picard iteration:
tolerance =1 e -9
integrator = TaylorPicardIn te gr at or ( tolerance )flow = integrator . flow_step (f , dom , h )print ( flow )Computed solution:
VectorFunctionPatch(result_size=1,dom=[{1.0:2.0},{-0.03125:0.03125}], rng=[{0.96542248:2.0350336}])[ {0.3272*x0^6*x1^5 -2.593*x0^5*x1^5 +8.408*x0^4*x1^5 -14.52*x0^3*x1^5 +14.17*x0^2*x1^5 -7.419*x0*x1^5+1.629*x1^5 +0.04681*x0^8*x1^4 -0.5617*x0^7*x1^4 +2.778*x0^6*x1^4 -7.574*x0^5*x1^4 +12.80*x0^4*x1^4-13.90*x0^3*x1^4 +9.505*x0^2*x1^4 -3.747*x0*x1^4 +0.6524*x1^4 -0.01165*x0^8*x1^3 +0.1205*x0^7*x1^3-0.4802*x0^6*x1^3 +0.9949*x0^5*x1^3 -1.317*x0^4*x1^3 +1.142*x0^3*x1^3 -0.6297*x0^2*x1^3 +0.2012*x0*x1^3-0.02840*x1^3 -0.0003479*x0^9*x1^2 +0.004808*x0^8*x1^2 -0.02333*x0^7*x1^2 +0.04241*x0^6*x1^2-0.02796*x0^5*x1^2 +0.07998*x0^4*x1^2 -0.06351*x0^3*x1^2 +0.03335*x0^2*x1^2 -0.01044*x0*x1^2+0.001481*x1^2 +0.00002472*x0^8*x1 -0.0003107*x0^7*x1 +0.0003195*x0^6*x1 +0.007722*x0^5*x1+0.0008026*x0^4*x1 -0.1674*x0^3*x1 +0.0004201*x0^2*x1 -0.0001451*x0*x1 +0.00002249*x1 +1.000*x0+/-0.000000000470} ]Kyoto, 13 February 2023

Rigorous Function Calculi – 31 / 48

Function sets
A powerful way of defining sets is using functions to express constraints, and tomap points from a parameter domain into the space of interest:

R = {x | f (x) ∈ C};

S = {h(p) | p ∈ D | g(p) ∈ E}.
Intersection and image for constrained image sets S are expressible using
function composition:

S ∩ R = {h(p) | p ∈ D | g(p) ∈ E ∧ [f ◦ h](p) ∈ C};
f (S) = {[f ◦ h](p) | p ∈ D | g(p) ∈ C}.

Constrained image sets are overt and compact, meaning that they have verifiableintersection and subset relations with open sets.
?

Reduces to nonlinear feasibility problems of the form {p ∈ D ∧ g(p) ∈ C} = ∅.

Kyoto, 13 February 2023

Rigorous Function Calculi – 32 / 48

Hybrid system — water tank
The water height h in a tank with continuous outflow and a valve-restricted inflow
needs to be controlled between hmin and hmax .

√By Torricelli’s law, ḣ = −a h + b α where α ∈ [0, 1] is the aperture of the inlet
valve, a and b are physical constants.
The valve can be opened or closed at a speed of 1/T and it is controlled so that it
starts to open as soon as h 6 hopen and starts to close as soon as h > hclose .
The automaton starts in location opening, the value of α increases with speed
1/T and the water level follows Torricelli’s law. As soon as α = 1 (the valve isfully open) the urgent transition to open is taken. The valve is kept open until
h = hclose , when the automaton switches to location closing and the valve starts
whenclosing with speed −1/T . The urgent transition to closed is activated √α = 0, and the water level decreases following the dynamics ḣ = −a h until
h = hopen and the transition to opening is taken.

Kyoto, 13 February 2023

Rigorous Function Calculi – 33 / 48

Hybrid system — water tank

Kyoto, 13 February 2023

Rigorous Function Calculi – 34 / 48

Hybrid system — water tank
A computation of one evolution loop starting from {(opening, 0, 0)} using Ariadne
is shown.
The result is a rigorous and accurate over-approximation of the exact reachable
set.

(Computation time ≈ 6s on a 2.4 GHz Intel Core 2 Duo with 4 Gb of memory.)

Kyoto, 13 February 2023

Rigorous Function Calculi – 35 / 48

Introduction
Motivation
TheoryPractice
Problems• Naming• Semantics• Parameters• Dispatching• Improvements• Extensions• Verification• Documentation

Open Problems

Conclusions

Kyoto, 13 February 2023

Rigorous Function Calculi – 36 / 48

Naming conventions
Good naming schemes make code easier to write and understand!

• Should Function<Real(Real)> be Effective by default?

• Should Function be required to be total, or allowed to be partial?

• Should Function be parametrised by the type of the domain set?

• Should we have separate ContinuousFunction, SmoothFunction,

ElementaryFunction classes?

Kyoto, 13 February 2023

Rigorous Function Calculi – 37 / 48

Choices for semantics
What does Validated information mean? The former is stronger than the latter:• Information giving an (open or compact) set x̂ of possible values of x : X.

• A finite prefix w ∈ Σ∗ of a name p ∈ Σω of x : X under a representationδ : Σω → X.

How can we mix Effective and Validated information in computations?

• Usually, one can choose a sensible Validated answer.

• There are corner cases (like evaluating a constant function) where this fails!

Should we consider a separate Symbolic information for problem specification
without an Effective implementation?

Kyoto, 13 February 2023

Rigorous Function Calculi – 38 / 48

Parametrising computations
Builtin functions like mul and sin, and functionals like flow of a differential
equation, have many possible implementations.
How can we select which implementation to use within a particular computation?
The more complicated solvers have individually-tunable accuracy parameters.
How can we tune computational parameters in solvers hidden within a computation?

• Since the final answer may be a very simple object e.g. a Kleenean, we do
not have access to these parameters at the top layer.

• Currently, computing both Kleenean and Real is controlled by an Effort
parameter. This could correspond to the bit-length of the computed name.

• Currently, data types specify their own “core” operations (e.g.

Bounds<FloatMP> uses rounded MPFR operations for elementary functions.
• In metric spaces only, controlling the Accuracy makes sense.
• Maybe solver parameters should be mappings of Effort to the actual values.
• Ultimately, designing efficient algorithms relies on estimating the trade-off
between accuracy and computational complexity.

Kyoto, 13 February 2023

Rigorous Function Calculi – 39 / 48

Double-dispatching
Most types have multiple equivalent representations. How do we combine these?
Object-oriented languages provide polymorphic dispatching to implement genericclasses.

• Generic classes, such as ValidatedReal can dynamically store objects ofdifferent types, such as Bounds<Dyadic> and Ball<FloatMP,FloatDP>.

• Dispatching (a fixed set of) unary operators in C++ is easily accomplishedusing virtual functions.

• Dynamically dispatching binary functions requires double dispatching• Unfortunately, there’s no good solution to double-dispatching in C++.
• In A RIADNE, we provide wrappers which can semi-automaticallydoubly-dispatch certain pairs of arguments.

Currently, the Number classes support polymorphic dispatching, including
double-dispatched binary arithmetic, while the Real class only used symbolic
operations. Is this distinction worth keeping?
Kyoto, 13 February 2023

Rigorous Function Calculi – 40 / 48

Clarity and efficiency improvements
There are many different idioms used in the code. These should be cleaned-up
and made consistent.
Modern C++ has support for concepts. These are being integrated into thecode-base.
Some of the core algorithms need efficiency improvements, notably the differential
equation solvers.

Kyoto, 13 February 2023

Rigorous Function Calculi – 41 / 48

Extensions to functionality
Implement types for piecewise-continuous and measurable functions.
Implement concrete function classes based on Chebyshev and Fourier
expansions.
Implement set-valued functions.

Kyoto, 13 February 2023

Rigorous Function Calculi – 42 / 48

Verification of implementation
We should also implement core functionality of an A RIADNE-like tool in a language
which allows verification of the algorithms and their implementation.

• Together with the groups of N. Müller, M. Konečný, M. Ziegler, H. Thies,

S. Park, and A. Simpson/A. Bauer we have been looking at designing a verified
language for effective (exact) real computation.

• N. Müller has been attempting to verify iRRAM code.

• M. Konečný has been experimenting with an Agda implementation of AERN
functionality.

• H. Thies, S. Park, H. Tsuiki and M. Konečný have been implementing TTE
functions and sets in Coq.

• P.C. M. Niqui and N. Revol have verified the operations of Taylor model
calculus in Coq (2011).

Kyoto, 13 February 2023

Rigorous Function Calculi – 43 / 48

Documentation
The main A RIADNE documentation is made with Doxygen.
It’s easy to make appallingly bad documentation with Doxygen!
Nevertheless, we’ve tried to+ Separate pages decribing the the theoretical basis of the tool and thealgorithms.

+ High-level overview of modular structure and functionality of each module.
+ Links to annotated tutorials.
− Unfortunately the documentation is also out-of-date in places.
We really need specific information from users to help improve the documentation!

+ Consider using Jupyter notebooks for Python tutorials. [Holger’s suggection!]

Kyoto, 13 February 2023

Rigorous Function Calculi – 44 / 48

Introduction
Motivation
TheoryPractice
Problems
Conclusions• Summary• Future Work

Kyoto, 13 February 2023

Concluding Remarks

Rigorous Function Calculi – 45 / 48

Summary
Function calculi can be developed based on concepts from computable analysis
and techniques from rigorous numerics.
A RIADNE is a software tool for reachability analysis of nonlinear hybrid systems
based on such a rigorous function calculus.

• The functional calculus includes support for interval arithmetic, linear algebra,
automatic differentiation, function models with evaluation and composition,
solution of algebraic and differential equations, constraint propagation andnonlinear programming.

•

Kyoto, 13 February 2023

Rigorous Function Calculi – 46 / 48

Future Work
Improvements to the design, implementation and efficiency of A RIADNE.
Implementation of function calculi for discontinuous functions in A RIADNE.
Formal verification of a computable, rigorous function calculus in Coq.

Kyoto, 13 February 2023

Rigorous Function Calculi – 47 / 48

Future Work
Improvements to the design, implementation and efficiency of A RIADNE.
Implementation of function calculi for discontinuous functions in A RIADNE.
Formal verification of a computable, rigorous function calculus in Coq.

Your ideas for collaboratative research at RIMS!

Kyoto, 13 February 2023

Rigorous Function Calculi – 47 / 48

Introduction
Motivation
TheoryPractice
Problems
Conclusions

That’s all, folks!

Kyoto, 13 February 2023

Rigorous Function Calculi – 48 / 48