Gabriel Sarazin

Amandine Marrel

Sébastien Da Veiga

Vincent Chabridon

Vincent Chabridon What

What is hidden behind the Sobolev kernels involved in the HSIC-ANOVA decomposition?

Speaker: Gabriel Sarazin

What is hidden behind the Sobolev kernels involved in the HSIC-ANOVA decomposition?

Gabriel Sarazin CEA, DES, IRESNE, DER, Cadarache F-13108 Saint-Paul-Lez-Durance, gabriel.sarazin@cea.fr Amandine Marrel CEA, DES, IRESNE, DER, Cadarache F-13108 Saint-Paul-Lez-Durance, Institut de Mathématiques de Toulouse, F-31062, Toulouse, amandine.marrel@cea.fr Sébastien Da Veiga Safran Tech, Modeling & Simulation Rue des Jeunes Bois, Châteaufort, F-78114 Magny-Les-Hameaux, sebastien.daveiga@gmail.com Vincent Chabridon EDF R&D 6 Quai Watier, F-78401, Chatou, vincent.chabridon@edf.fr In many industrial fields, physical phenomena are commonly modeled by numerical simulators. These simulation codes can take a high number (dozens, if not hundreds) of uncertain parameters as input variables. Since computer-based simulation experiments are time-consuming, the permissible number of runs is most often very limited. In this context, the benefits brought by a preliminary global sensitivity analysis are twofold. On the one hand, the screening step allows to select a smaller amount of influential inputs in order to reduce dimensionality or to support the construction of a surrogate model. On the other hand, the influence of each input variable is quantified and the associated rankings provide a relevant piece of information that can be used in the sequential building process of a metamodel as in [4]. In view of these two purposes, total-order Sobol indices are often regarded as the most adapted indicators. Unfortunately, unlike their first-order counterparts, their estimation requires a too large number of simulations.

To remedy this difficulty, a new class of sensitivity measures built upon the theory of reproducing kernel Hilbert spaces (RKHS) has emerged over the last decade [1]. Each input variable is assigned a continuous positive definite kernel and is equipped with the related RKHS. The same is done for the output variable. The distance between the joint input-output bivariate distribution and the bivariate distribution obtained under independence is measured through the distance of their respective embeddings in the tensorized RKHS, which leads to the Hilbert-Schmidt Independence Criterion (HSIC) [3]. Since HSIC indices can be expressed by means of kernel-based statistical moments, their estimation can be achieved at minimal computational cost and without estimating the joint input-output distribution. However, HSIC indices are often blamed for their lack of interpretability because their sum is not equal to 1. As a consequence, they do not fit into the advantageous framework that goes with analysis of variance (ANOVA). To bridge this gap, it has been recently proved that the use of orthogonal kernels enables an ANOVA-like decomposition for HSIC indices [2]. When the input variables are uniformly distributed, Sobolev kernels (parametrized by a smoothness parameter r) match all the required conditions. They allow to define higher order HSIC indices (especially total-order HSIC indices) and therefore to rigorously separate main effects and interactions.

The main objective of this work is to provide new insights into Sobolev kernels which have indeed been very few studied so far. First, we demonstrate that Sobolev kernels are characteristic. This ensures that the nullity of Sobolev-based HSIC indices is equivalent to independence. Then, another notable contribution of this work is the identification of one explicit feature map for Sobolev kernels, which may help the user understand the spectrum of dependency patterns that can be captured by this novel HSIC index. As with most kernels, the canonical feature map (resulting from Aronszajn's theorem) is uninformative. An additional pitfall comes from the fact that Sobolev kernels are not shift-invariant, which prevents the use of Bochner's theorem and the following characterization in terms of Fourier transform. Instead, another kind of feature map stems from Mercer's theorem which asserts that any continuous kernel on a compact set may be rewritten as an infinite sum that only implies the eigenvalues and eigenfunctions of the corresponding kernel integral operator. This key result discloses a feature map that sends the unit interval into the Hilbert space of real-valued square-summable sequences. The eigenvalue problem arising from this new feature map is tackled in two different ways. Firstly, we investigate the benefits of kernel feature analysis where eigenvalues and eigenfunctions are estimated from Gram matrix simulation. It allows to visualize eigenfunctions and sometimes to come up with closed-form candidate functions. Secondly, we demonstrate that the eigenvalue problem is equivalent to a Cauchy problem consisting of a linear homogeneous ordinary differential equation with constant coefficients and a sufficient number of boundary conditions. Then, two situations deserve a specific study. When r = 1, an exact analytical solution is available. On the contrary, when r = 2, the Cauchy problem cannot be completely solved. A solution in the form of a linear combination of analytical functions is obtained but the coefficients cannot be retrieved. An asymptotic approximation is then derived to accurately estimate small eigenvalues while robust numerical resolution is achieved by use of a semi-analytical adhoc calibration algorithm.

[1] Sébastien Da Veiga. Global sensitivity analysis with dependence measures. Journal of Statistical Computation and Simulation, 85(7):1283-1305, 2015.

[2] Sébastien Da Veiga. Kernel-based ANOVA decomposition and Shapley effects: application to global sensitivity analysis. arXiv preprint arXiv:2101.05487, 2021.

[3] Arthur Gretton, Olivier Bousquet, Alex Smola, and Bernhard Schölkopf. Measuring statistical dependence with Hilbert-Schmidt norms. In International conference on algorithmic learning theory, pages 63-77. Springer, 2005.

[4] Amandine Marrel, Bertrand Iooss, and Vincent Chabridon. Statistical identification of penalizing configurations in high-dimensional thermalhydraulic numerical experiments: the ICSCREAM methodology. arXiv preprint arXiv:2004.04663, 2020.


Djamel Merad

Tracking 3D free form object in video sequence

Keywords: Tracking, 3D free form object, Skeletonization, Graph matching

In this paper we describe an original method for the 3D free form object tracking in monocular vision. The main contribution of this article is the use of the skeleton of an object in order to recognize, locate and track this object in real time. Indeed, the use of this kind of representation made it possible to avoid difficulties related to the absence of prominent elements in free form objects (which makes the matching process easier). The skeleton is a lower dimension representation of the object, it is homotopic and it has a graph structure. This allowed us to use powerful tools of the graph theory in order to perform matching between scene objects and models (recognition step). Thereafter, we used skeleton extremities as interest points for the tracking.

Introduction

Object tracking is an important field of artificial vision. It is a difficult problem involved in tasks as motion evaluation or automatic visual control. The possible applications are from medical imagery to security, including road traffic analysis, etc.

Most of the algorithms allowing 3D pose estimation from a 2D images sequence are based on the same principle. The first step consists in recognizing the type of target in order to choose the best model from a database (initialization). The second step, performed during the video stream processing, consists in target tracking.

An initial pose of the object is known either by an initialization process before tracking, or by a prediction process using the results obtained from the previous image. This initial pose allows to project the target's model onto the image's plane. Starting from this projection, it is possible to know if the pose is accurate or not. A comparison measure between the 2D image and the projected model allows a pose adjustment (fitting).

The general principle of tracking is common to all the methods. Nevertheless they significantly vary according to their matching method or prediction techniques. Among the different existing tracking methods, we will cite the one of Lowe [START_REF] Lowe | Three-dimensional object recognition from single two-dimensional images[END_REF] [12], of Gennery [START_REF] Gennery | Tracking known threedimensionnal objects[END_REF] [6] and of Harris [START_REF] Harris | Rapid, a video rate object tracker[END_REF] which we can qualify as generic methods. These three tracking methods give very good performances and use algorithms of different complexities.

Some other methods were developed since. Some model based approaches use the automatic visual control to track the salient edges of a model [START_REF] Drummond | Real-time visual tracking of complex structures[END_REF] while some other methods use primitives as points, lines, or ellipses [START_REF] Comport | A real-time tracker for markerless augmented reality[END_REF] [START_REF] Marchand | Robust real-time visual tracking using a 2d-3d model-based approach[END_REF]. However, the presented methods were been developed for polyhedral models having straight edges, which facilitate their detection. That is the reason why these methods can not be applied to free form object tracking.

To track free form objects, some methods use a priori and a posteriori information. These methods either will index several images (shots) of the target object depending on some points of interest and store them in a database, or will build a model starting from this camera shot. The objects are tracked using "patches" together with the indexed images for a homographical matching [START_REF] Vacchetti | Stable realtime 3d tracking using online and offline information[END_REF]. Another way to manage free form object tracking is to use a poor sparse metric environment model in order to perform real time model tracking using the camera [START_REF] Skrypnyk | Scene modelling, recognition and tracking with invariant image features[END_REF]. Unfortunately, these systems need a learning step on several shots of the real object, they can not use the 3D model only.

Most of the effective tracking methods were developed for polyhedral objects. It would be enough to extract a polyhedral primitive of our free form object and to adapt the existing tracking algorithms. We succeeded in it by using skeletons for the free form object representation. Indeed, a skeleton is an object representation in an inferior dimension. It has the benefit to preserve most of the topological and geometrical information of the silhouette. The skeletonization or the object representation by a set of segments allows us to bypass the problems related to free form objects.

In image processing and computer vision in general, there are many applications for 2D and 3D ob-jects skeletons (encoding, compression etc.). Further to our interesting results obtained with this formalism in free form object recognition and localization [START_REF] Merad | A skeleton based method for efficient 3d object localization: Application to teleoperation[END_REF] [18], we explored its application to track this kind of objects.

More precisely, we adapted the Lowe method [START_REF] Lowe | Robust model-based motion tracking through the integration of search and estimation[END_REF] to track a 2D skeleton using a 3D skeleton, therefore to track a 3D object in a 2D shot. The points of interest we used are the skeleton extremities. The initialization step corresponds to the target pose computing, without any a priori knowledge of the scene. The initialization system is described in section [START_REF] Comport | A real-time tracker for markerless augmented reality[END_REF] and detailed in [START_REF] Merad | A skeleton based method for efficient 3d object localization: Application to teleoperation[END_REF] [START_REF] Merad | A 3d free form object localization using skeletons: Application to teleoperation[END_REF]. Once the initial position of our object known, it is possible to start tracking.

The 3D skeleton is projected on the image plane, starting from its initial position, by using previously calculated rotation and translation matrices. The algorithm uses the fact that there is a small displacement between two following frames in a video. Therefore the positions of 2D skeleton points in the two images will be very close to each other. A comparison measure between the 2D skeleton of the image and the projected 3D skeleton of the model is used for a pose adjustment. This step is presented in section (3) and the obtained results are presented in section (4).

Initialization

In this section we present the identification and the localization of 3D free form objects, starting from a 2D image of scene, where partial occlusions may appear.

We analyzed the architecture of 3D free form more popular recognition systems and we noticed the problems to which they are confronted with. Some of these difficulties are: to build the model for complex objects, primitives extraction allowing a correct matching and the lack of salient elements allowing a robust pose estimation. We also can cite complexity and occlusion problems which have not a satisfactory solution.

The approach we developed (figure 1), is based on using a 3D skeleton of an a priori known model and the on-line extraction of the 2D skeleton of an image of the object in front of the camera. The use of this formalism allowed us to solve the problems presented before: pose estimation and robustness to occlusions and auto-occlusions. Our method is essentially composed of three steps.

• Primitives modelization and extraction: search for an appropriate representation and primitive detection in the image/model. G1 ←-SkeletonToGraphTransform(S); 4.

hypothesis ←-Indexation(G1, {G2}); 5.

for each h from {hypothesis} do 6.

h ←-HypothesisChecking(G1, h); 7.

if h = null then 8.

OverlayImageWithModel; 9. exit; 10. else 11.

h ←-NextHypothesis; 12.

display("No corresponding model for this image");

Modelization and primitives extraction

The skeletonization represent an important step of our method. Since this concept appeared as form descriptor, many skeletonization algorithms have been proposed in the literature. The different skeletonization techniques can be classified in two categories. The discrete methods, as thinning, "grassfire" potential fields, the distances maps [START_REF] Kegl | Piecewise linear skeletonization using principal curves[END_REF] [15] [START_REF] Nilsson | Finding the minimal set of maximum disks for binary objects[END_REF] [25] and the continuous methods, based on Voronoi diagrams [1] [4].

The discrete domain algorithms are very popular, mainly because of their simplicity. Unfortunately, Figure 3: Hull object and 2D skeleton the obtained skeletons are often sensitive to image rotation and the connectivity is not always preserved and often needs a post-processing. In the continuous domain, the skeleton is a sub-graph of a Voronoi diagram. The points of the object shapes represent a sampling of the continuous surface. The Voronoi diagram is built starting from this sampling in order to extract the skeleton. This method has several advantages. Only the shape points are needed, which allows to considerably reduce the number of points to be processed. The obtained skeleton is connected and topologically equivalent, because the shape representation of the objects implicitly preserve these aspects. Not using a discrete coordinates grid allows us to use the Euclidean distance, therefore the skeleton is not sensitive to image rotation. Another benefit is the skeleton has a graph structure (without post-processing), easy to build.

Unfortunately there are also some drawbacks for this kind of methods. The main problem is the choice of a sub-graph that represents the nearest approximation of the skeleton. Another problem of the continuous domain comes from the construction of the Voronoi diagram. The diagrams presented in this article are always limited to non-aligned and no co-circular seeds. For real images we have a great number of such cases. This forces us to move the seeds, leading to a bad localization of the skeleton. Furthermore, the skeletons are very sensitive to noise induced by the sampling. Therefore an extra-processing has to be performed to prune the unwanted edges.

Our proposed hybrid method compensate for these shortcomings [START_REF] Merad | Skeletonization of two-dimensional regions using hybrid method[END_REF]. In this method we combined the two techniques: the one based on distance map and the one based on Voronoi diagram, in order to benefit from their strong points. This skeletonization algorithm has been tested on different images. The experimental results showed that it guarantees the complete connectivity of the skeleton, the homotopy, the geometric invariance and also a good localization. The figure 3 shows the results of this technique on real images.

In order to obtain an efficient matching, it is es-sential to have a close representation of the model and of the image. In our approach, the skeletonization results in two graphs named 3D and 2D. The 3D graph is obtained from the 3D homotopic skeleton which is itself computed from the 3D model. This represents an off-line operation. An on-line skeletonization of 2D scene image (second line of figure 2) allows to build the 2D graph (line 3 of figure 2). Our method encodes these graphs in such a way that each skeleton graph store topological and geometrical information of the initial form. Indeed, each node of the graph represents a skeleton element that is a 2D segment for the 2D graph and a 3D segment for the 3D graph. Each edge in the graph represents a link between the different elements of the skeleton (figure 4).

Matching

The recognition consists in establishing an isomorphism between the 3D graph from the database and the 2D graph obtained from the image. The problem is to establish a good quality measure for graph matching. In our application, this task is difficult because the measure has to evaluate the similarity degree between similar structures of two sub-graphs, in order to manage the occlusion and shadow problems. Therefore, it has to allow, in the recognition step, to quickly choose the best model from the database ( l i n e4i nfi g u r e2 )a n da l s ot op e r f o r mao n et oo n e node matching, in order to locate and validate this h y p o t h e s i s( l i n e s5t o1 2i nfi g u r e2 ) .

After studying several graph matching algorithms, we implemented an isomorphism-based method using the topological signatures of each node. This approach was used by Shokoufandeh [START_REF] Shokoufandeh | Applications of bipartite matching to problems in object recognition[END_REF], Siddiqi [START_REF] Siddiqi | Shock graphs and shape matching[END_REF] and Macrini [START_REF] Macrini | View-based 3-d object recognition using shock graphs[END_REF] to match two shock graphs. We introduced several modifications on this method [START_REF] Merad | A skeleton based method for 3d freeform object recognition[END_REF] in order to obtain a robust matching algorithm at the final nodes level and also to carry out the "polygamous" nodes i.e. a node can be matched at the same time with several nodes.

In order to topologically describe a tree we used graphs eigenspaces. We remember that each graph can be represented as a adjacency matrix of {0, 1}. The values 1 give adjacent nodes of the graph (and 0 on the diagonal). The eigenvalues (EV) of the adjacency matrix of the graph store some significant structural properties of the graph (tree).

Explicitly, let T be a tree of maximum degree Δ(T )andT 1 ,T 2 ...,T s the sub-trees of its root. For each sub-tree T i , the degree of its root is δ (T i ). To obtain this signature, we compute the eigenvalues of each adjacency matrix of each sub-tree T i .L e tS i be the sum of δ (T i ) eigen values of T i . Then the ordered elements S i , become the components of the vector χ of dimension Δ (T ) called topological signature and assigned to the tree's root. If the elements number S i is lower than Δ (T ), then the vector is filled with 0. This method is recursively repeated in order to assign a vector to each node of the tree.

As we already mentioned, a trees isomorphism can not exist between the image skeleton and the model skeleton because of the occlusions or/and the noise. The solution consists in finding a maximal cardinality and minimal weight matching in a bipartite graph covering the nodes of the two skeletons. The bipartite graph is a graph where each edge is pondered by the topological distance. For more details see [START_REF] Merad | A skeleton based method for efficient 3d object localization: Application to teleoperation[END_REF] and [START_REF] Merad | A 3d free form object localization using skeletons: Application to teleoperation[END_REF].

We remember that our algorithm allows, in a first time, to index the object's skeleton in a database -which corresponds to the recognition step. In a second time, it allows a one-by-one matching of the nodes needed in the localization and verification step.

Projection and verification

In the precedent step, due to the noise, there are several matching hypothesis. To verify a hypothesis validity we will estimate the rigid transformation between the recognized (presumed) model and the camera. The model projection on the image allows to accept or reject this hypothesis. Each matched couple corresponds to a segment couple, one of the 2D skeleton and the other of the 3D curvilinear skeleton. We project the primitives of the 3D object model on the image plane. If the projected primitives do not coincide with the image primitives, according to a pre-established threshold, then the hypothesis is rejected and we treat the next hypothesis and so forth. If the points of the 3D skeleton projected on the im-age do not coincide with the points of the 2D skeleton, then the object is not recognized.

This verification module is based on the measurement of this error (1) which is the distance in the image (in pixels) between, each time, a 2D skeleton point and a 3D skeleton point projected on the image.

error = 1 N N i=1 (û i -u i ) 2 +(v i -v i ) 2 (1)
Where : N is the number of matching points, (u i ,v i ) the coordinates of a point image, (û i , vi ) estimation of the object's point after projection on the image plane.

The 3D pose estimation consists in finding the rigid transformation (R, T ) minimizing some calculated error (as the sum of error squares) of the one of two collinearity equations (in the image space or in the object space). The two methods generally used to solve this problem are the Gauss-Newton and the Levenberg-Marquardt methods [START_REF] Lowe | Fitting parametrized threedimensional models to images[END_REF].

We used the method of Lu and al [START_REF] Lu | Fast and globally convergent pose estimation from video[END_REF] named Orthogonal Iteration (OI) algorithm. Contrary to classical methods, used to solve the optimization problems on the whole, the OI algorithm cleverly exploits the specific structure of the 3D pose estimation problem.

To estimate the object's pose, this algorithm uses an appropriated error function defined in the object's space. The error function is rewritten in order to accept an iteration based on the classical solution of the 3D pose estimation problem, called absolute orientation problem.

This algorithm gives exact results and converges quickly enough, therefore it is very interesting for real-time applications.

The figure 5 presents some initialization results. From left to right, there are the real object and the 3D model recognized by our algorithm and overl a y e do nt h er e a li m a g e . T h e r ea r e1 0m o d e l si nt h e database and the tests showed that our methods can correctly recognize and localize several 3D objects even in presence of occlusions or auto-occlusions and shadows (due to the strength of our subgraph matching algorithm ). For more results see [START_REF] Merad | A skeleton based method for efficient 3d object localization: Application to teleoperation[END_REF] and [START_REF] Merad | A skeleton based method for 3d freeform object recognition[END_REF]

3T r a c k i n g

The initialization step corresponds to target pose computing, without a priori knowledge of the scene. Once the initial position of our object known, the tracking process can start. {R, T}← -InitialValues ; 3.

while tracking not finished do 4.

i ←i +1; 5.

{E1}←-3DSkeletonEndpointsProjection; 6.

{E2}←-2DSkeletonComputation( Ii); 7.

for each Endpoint e1 of E1 and e2 of E2 do 8.

d ←-DistanceComputation(e1,e2); 9.

{(2D, 3D)}← -CoupleChoice(d min); 10.

{R,T}← -NewPoseComput ({(2D, 3D)}); The original position we used is the one computed on the first image during the initialization step (line 2 of figure 6). The 3D skeleton is then projected on the image plane, starting from this position, by using the rotation and translation matrices calculated before (line 5 of figure 6). The algorithm is based on the fact that there are small movements between two video frames. The positions of 2D skeleton points will be very close from one frame to another.

We call p ij the positions of the extremities j of the 2D skeleton in the frame i, numbered from 1 to n i (their number can vary with the frame). We call q k the positions of the projections on the image plane of the 3D skeleton extremities, built using rotation and translation matrices R and T respectively.

We suppose the initialization was correct, therefore at the time of the projection of the 3D skeleton, the points q k which were used to compute the position, were near their matching points p ij .I nt h e next frame (number (i + 1)), the points p (i+1)j are different from points p ij . Anyway, the deviation of the object being small, there is also a small variation for the 2D skeleton position. This way the points q k which, by construction, were very close to some points p ij , will be very close to some points p (i+1)j . The movement of our object being very small, the old and the new 2D skeleton are very close one of the other, therefore theirs extremities are also very close.

By measuring the pixel differences separating the points p (i+1)j and q k , and sorting in ascending order, we can form a number of pairs (p (i+1)j ,q k ). Whereas during initialization matching between the extremities of the 2D and 3D skeleton required heavy computing, due to the fact that the position was unknown, it is almost immediate here by considering t h ep o i n t s3 da tt h eo r i g i no fq k (lines 6 to 9 in figure 6). The orthogonal iteration algorithm we used for initialization is applied on the pairs (2DPoint , 3DPoint). The 3D points being obtained by replacing q k by the original points of the 3D skeleton. This provides us a new value for R and T , which could be used in new computations on the next frame (line 10 of figure 6).

We improved the Lowe method by integrating the step of the graph matching in the tracking loop. For example, in the case of an abrupt and irregular movement or of an important occlusion, the model loses the object. While applying, our algorithm of graph matching we succeed to quickly re-localize this object. We remind that in this step we already know which 3D graph corresponds to the 2D graph, in that case, it is sufficient to make a one-to-one nodes matching of two graphs. This step is very fast because the skeletons graphs are small. This procedure is not called during all the tracking process but only when the projection error is greater then threshold. The algorithm of the figure 6 becomes as the algorithm of the figure 7. Tracking was implemented in C++ in order to realize real-time computing. There is a file in specific format allowing the passage from the initialization step, realized in Matlab, to tracking. This file contains the rotation and translation matrices corresponding to the calculated position, as well as the model of the 3D skeleton. The tests were realized on synthetic images in order to confirm our hypothesis. Using synthetic images allows to compare and estimate projection errors more precisely, because we know the exact pose used to generate our video. Moreover, that let us to generate various types of movement in various configurations. Let us specify that although the used object is a virtual one, its tracking in a sequence of frames is made realistic (difficult) by simulating shades, reflections, noises, self-occlusions. The videos were built with POV-Ray software to generate the images and Virtual Dub to realize the video editing. To read the video stream we used Intel's OpenCV library. Its main advantage is that it integrates some image processing functions. The skeletonization is the heart of our application, since it is during this operation that the points used to tracking will be extracted. Contrary to the initialization part, where we emphasized the extraction of a robust skeleton by a hybrid method, for the tracking part, we privileged the speed (computing time). For this reason the skeleton is calculated starting from the Delaunay diagram described in [START_REF] Georges | Triangulation de delaunay et maillage : Applications aux éléments finis[END_REF]. This method is particularly fast, for a video with 19 images per seconds, the extraction of the skeleton is carried out in real time on a Pentium M 1600 MHz. We performed several tests and the following paragraphs will only present a part. The figure 8 shows the tracking result on the cactus object, in a 68 frames sequence. This figure represents, from left to right and from up to down, the frames 1, 8, 16, 24, 32, 40, 48, 54, 67 respectively. We present here the most complicated case. We simulated a complex movement with a strong rotation around axis Z (optical axis of the camera). The task becomes difficult because of the disappearance of a part of the object. The tracking is carried out at the rate of 19 Hz by a wireframe model (in magenta). The points of interests i.e. the extremities of the skeletons are in red for the 3D points, and in green for the 2D points. The graph in the figure 9(a) represents the reconstruction errors (in pixels) for this example. On the whole, the tracking is carried out correctly, it happens sometimes that the model loose the object and recapture it in the next sequence -image 2 of the figure 8 (frame 8) (due to the integration of the graph matching procedure). We can also observe this on the graph of the figure 9(a).

Another example is presented in the figure 9(b) which represents the reconstruction errors of the first However, in some videos, the algorithm is not able to efficiently track the object (figure 13). This is essentially due to the lack of a sufficient number of matchings occurring in the case of occlusions or selfocclusions. The 3D pose parameters are not correctly computed. We remember that the used orthogonal iterations algorithm needs at least four good noncoplanar matchings. However, after few frames our algorithm recaptures the object and tracks it efficiently due to the integration of the graph matching procedure.

Conclusions

In this article we presented a model based method for 3D free form object tracking. The results of adapting our skeleton recognition method to object tracking were encouraging and justify our approach. Nevertheless we are aware that some improvements are needed for efficient tracking. In our future works we will focus on a 2D robust and fast skeletonization method. We also project to compensate the weakness of the Lowe method by adding a prediction step in our algorithm.

Figure 1 :

 1 Figure 1: 3D free form object recognition by using skeletons

Figure 2 :

 2 Figure 2: 2D/3D identification using skeletons

Figure 4 :

 4 Figure 4: Graph construction from the skeleton

Figure 5 :

 5 Figure 5: Cactus object recognition among several objects

Figure 6 :

 6 Figure 6: Tracking algorithm using skeletons

Figure 7 :

 7 Figure 7: New tracking algorithm

Figure 8 :

 8 Figure 8: Object tracking sequence

  (a) During tracking in -R-(b) During tracking in -RT-

Figure 9 :

 9 Figure 9: Reconstruction error in pixels

Figure 10 :

 10 Figure 10: Frame 27

Figure 12 :Figure 13 :

 1213 Figure 12: Frame 78


Sofiane Bessaï
email: sofiane.bessai@etudiant.univ-perp.fr

Ben Dorra

Khalifa

Hanane Benmaghnia
email: hanane.benmaghnia@univ-perp.fr

Matthieu Martel
email: matthieu.martel@univ-perp.fr

Fixed-Point Code Synthesis Based on Constraint Generation ‹

Keywords: Fixed-point arithmetic, code synthesis, precision tuning, linear programming, static analysis

Fixed-point arithmetic is a well-known alternative to floatingpoint arithmetic on embedded systems. It is used to reduce some computation costs in terms of speed and power consumption on certain platforms, e.g. medical devices, cars, and robots. In this article, we present POPiX, a novel fixed-point program synthesis tool based on static analysis. The originality of our method is to solve a system of constraints generated from the program source code. Thus, the solution of our constraints gives the new fixed-point formats while accomplishing the accuracy required by the user. Basically, POPiX takes as input an imperative program running in floating-point arithmetic and synthesizes a new program coupled to a fixed-point library relying on integers only. We evaluate POPiX on a collection of floating-point benchmarks coming from FPBench. Results demonstrate the efficiency of our analysis by achieving memory savings up to 75% with energy savings up to 3.5ˆ.

Introduction

Floating-point arithmetic is the dominant approximation to represent a large spectrum of real numbers. Although it offers better precision, programmers do not always need the high level of accuracy offered by the largest floating-point formats. In addition, owing to its complex internal circuitry and the increased memory requirements, floating-point arithmetic can be exorbitant in terms of speed and power consumption on certain platforms such as mobile phones, video game consoles, and digital controllers. To bridge this gap and since many embedded architectures can be implemented using very low bit-width numbers, the solution is to deploy the fixed-point arithmetic as an alternative to the floatingpoint one as it can be efficiently realized using integer arithmetic. Fixed-point numbers in a certain format maintain a fixed divisor (so the name fixed-point).

With this kind of arithmetic, the number of bits splits into two parts, respectively named integer part and fractional part with a radix point that lies between them. Besides, many fields have revived their interest in fixed-point arithmetic when searching for cost effective hardware processors with less design effort. For instance, machine learning algorithms and models have been recently implemented using fixed-points with little accuracy loss [START_REF] Gupta | Deep learning with limited numerical precision[END_REF][START_REF] Lin | Fixed point quantization of deep convolutional networks[END_REF]. The conversion of code designed for embedded systems from floating to a fixed-point equivalent version is a long-established problem addressed in the litterature [START_REF] Cattaneo | Embedded operating system optimization through floating to fixed point compiler transformation[END_REF][START_REF] Cherubin | Dynamic precision autotuning with TAFFO[END_REF][START_REF] Darulova | Synthesis of fixed-point programs[END_REF][START_REF] Deest | Toward scalable source level accuracy analysis for floating-point to fixed-point conversion[END_REF][START_REF] Martel | Trade-offs of certified fixedpoint code synthesis for linear algebra basic blocks[END_REF]. Nevertheless, rarely we found tools (except for [START_REF] Cherubin | Dynamic precision autotuning with TAFFO[END_REF]) that deal with adjusting the floating-point formats in the original program before the conversion pass. Practically, programmers do not always need the high level of accuracy in the floating-point formats. However, manually adapting the precision of the variables may require considerable programming skill and application domain expertise. Given this consideration, automating the task of adjusting the program variables precision to improve its performance characteristics, before conversion, can help programmers to achieve performance benefits. Generally, this process of adjustment, also called precision tuning, involves a user requirement of accuracy and a semantic analysis of the program. The benefit of this tuning phase is to provide an optimized mixed-precision programs and pieces of information indicating the most suitable fixed-point data formats in the converted program.

The goal of this article is to propose POPiX: a new tool to transform a given numerical floating-point program into semantically equivalent one that exploits fixed-point computations with integers only. The key idea of this work relies on a semantic modelling of the numerical errors propagation throughout the floating-point program. In order to achieve the conversion, POPiX combines two fundamental steps. The first step consists in generating an Integer Linear Problem (ILP) from the original program in order to obtain the minimal formats (number of bits before and after the radix point) which fulfill the accuracy requirements. Basically, this is done by reasoning on the most significant bit and the number of significant bits of the values. The ILP problem can be optimally solved in one shot by a classical linear programming solver (LP) with no iteration. To the best of our knowledge, this is the first work that interests in statically synthetizing fixed-point code using an ILP formulation of the program. At the end of the tuning phase, the second step collects each information provided by the former step. The tool internally calls a fixed-point library to convert the associated indications into ones that exploit fixed-point computation with the number of bits required for each of the integer and the fractional parts. We evaluate POPiX on a set of benchmarks coming from the FPBench community. The POPiX source code, and all the data and results presented in this article are publicly available at: https://github.com/sbessai/popix. The remainder of this article is organized as follows. In Section 2 we present a motivating example describing our approach. Section 3 provides all the technical details about our new developed fixed-point code synthesis framework. Experimental results are presented in Section 4. We discuss related work in Section 5, and conclude in Section 6.

Overview

Before we dive into the technical details of our tool, we present in this section an overview of our method using the example of a FIR low-pass filter code given in Figure 1. The starting point of our analysis is to assume that all the variables are in a given IEEE754 precision (here we use single precision) and that a range is given for the inputs of the program. In addition, the statement require nsb(res,8) is a postcondition added by the user to specify that res must have 8 significant bits at the end of the execution. POPiX first performs a range determination by dynamic analysis for all the program variables at each control point. Based on semantic equations, POPiX generates an ILP problem from the program source code annotated with the results of the range analysis and the accuracy requirement. This yields a system of constraints:

C " $ ' ' ' & ' ' ' %
nsbp+q ℓ 548 ą" nsbpresq ℓ 549 , nsbpresq ℓ 549 ą" nsbpifq ℓ 551 , nsbpresq ℓ 511 ą" nsbpresq ℓ 545 , nsbpresq ℓ 545 ą" nsbp+q ℓ 548 `8 `carrypq ´8 nsbpresq ℓ 545 ą" nsbp+q ℓ 548 `8 `carrypq ´8, nsbpresq ℓ 548 ą" nsbpresq ℓ 549 , nsbpaux0q ℓ 547 ą" nsbp+q ℓ 548 `6 `carrypq ´8, nsbpaux0q ℓ 543 ą" nsbpaux0q ℓ 547 For instance, the system C of Equation (1) below describes the constraints generated for the addition and assignement statements for Line 23 of Figure 1 (left hand side). Some notations can be highlighted for the system of Equation [START_REF] Adjé | Fast and efficient bit-level precision tuning[END_REF]. First, each variable of our program is assigned to a unique control point ℓ P Lab in order to determine easily their number of significant bits. Second, the function carrypq is used to compute if a carry bit can occur through the operation (returns 0 or 1). Concerning scalability, we generate a linear number of constraints and variables in the size of the analyzed program (« 500 for the FIR low-pass filter code). The solution to our system of constraints gives the minimal number of bits needed with an accuracy guarantee on the results (highlighted in blue in the left hand side of Figure 1). If we take back Line 23 under discussion, the pair |8, 10| denotes that the unit in the first place of variable res is 8 whereas it has 10 significant bits. More details about the nature of constraints that we generate for the language of our input programs was detailed in [START_REF] Adjé | Fast and efficient bit-level precision tuning[END_REF].

Based on the tuning results, POPiX synthesizes the C code given in the right hand side of Figure 1. First, it selects the best format (int16 t, int32 t, etc.) for each variable (this is called mixed-precision). For example, at lines 6 and 7, vectors x, y and z are defined as int16 t variables while the vector derivCoeffs is defined as int8 t. The data type selected by POPiX for each variable is the minimal one enabling us to encode the fixed-point value following the formats coming from the ILP solution. For example, the variable res has 10 significant bits (lines 29 and 31 of Figure 1) and can consequently be encoded into int16 t. POPiX determines the initial formats xM, Ly of the variables occurring in the code and synthesizes the alignments needed to change the formats, before performing some operation. For example, the shifts performed at lines 29 and 30 are done in order to align the operand of the addition of Line 31. Similarly, the shift of Line 28 is done to obtain the right format for the result of the multiplication of Line 27. Currently, the fixed-point operation are generated sequentially and some additional optimizations could be done, for example by using only one shift for lines 27 and 29.

Floating to Fixed-Point Programs Synthesis

POPiX workflow is based on two frameworks as depicted in Figure 3: a developed fixed-point library (Section 3.1) and a precision tuning framework (Section 3.2). In the rest of this section, we explain how the combination of these features are achieved and which benefits they provide.

Fixed-Point Arithmetic

Since fixed-point operations rely on integer operations, computing with fixedpoint numbers is highly efficient for embedded systems with small memories and simpler CPUs. However, this arithmetic is more difficult to handle for the developer. There exists some fixed-point libraries such as Libfixmath 3 , Fixmath 4 Fig. 2: Fixed-point representation of a in a format xM, Ly . and FPM 5 , but we have developed our own library in order to dispose of all the features we need such as mixed-precision, elementary functions, etc.

A fixed-point number is represented by a sign s P t0, 1u, an integer value V P N in base 2 and a format xM, Ly as shown in Figure 2. The number of bits before, respectively after, the radix point is M P Z (respectively L P N). The value of a fixed-point number is obtained by multiplying the integer value V by the sign s and a scaling factor 2 ´L as follows:

a " p´1q s ˆV ˆ2´L .

(

) 2
Example 1 The fixed-point number a " 3 ă2,1ą corresponds to the value 1.5. Using Equation (2), we obtain a " p´1q 0 ˆ3 ˆ2´1 " 1.5.

Let us note that the number of bits M of the integer part already presented in Figure 2 is computed through the unit in the first place (ufp) defined by @x P F, ufppxq "

# minti P Z : 2 i`1 ą |x|u " tlog 2 p|x|qu if x ‰ 0, 0 if x " 0 . (3)
Hence, the number of bits M before the radix point is given by

M " ufpp|a|q `1 . ( 4
)
Let W be the number of bits used to encode a. The number of bits L of the fractional part of a is

L " W ´M ´1 . ( 5
)
The difficulty of the fixed-point representation is to manage the format xM, Ly manually against the floating-point representation which manages it automatically, thanks to the exponent. Let ˚be a fixed-point elementary operation with P {', a, b, c} and let us consider the fixed-point numbers a, b and c such that c " a ˚b. For the addition and subtraction, the resulting format of c is given by

xM c , L c y " xmax pM a , M b q `1, W c ´max pM a , M b q ´1y . ( 6
)
For the multiplication and division, the resulting formats of c are

xM c , L c y " xM a `M b , W c ´M a ´M b ą . ( 7
)
and The algorithms of these elementary operations are detailed in [11,[START_REF] Martel | Trade-offs of certified fixedpoint code synthesis for linear algebra basic blocks[END_REF]. Let us consider a fixed-point elementary function f P {abs, sqrt, sin, cos, arctan} and the fixed-point number c " f paq. For the square root function, the value of the result can be approximated by the digit recurrence iteration algorithm defined in [START_REF] Martel | Trade-offs of certified fixedpoint code synthesis for linear algebra basic blocks[END_REF]. For f P {sin, cos, arctan}, Taylor's formula is used to approximate the result. For example, the corresponding formula of sine is sinpaq « a ' pppa b aq b aq c 6q. The order of the Taylor series development depends on the number of significant bits needed for the result. In the following, we describe how we compute the fixed-point numbers occuring in our programs.

xM c , L c y " xM a `Lb , W c ´M a ´Lb y . (8)

Constraint Generation by Static Analysis

POPiX presents a novel static technique based on a semantic modelling of the propagation of the numerical errors throughout the code. In practical terms, our approach depends on two integer quantities: iq The ufp of the values (see Equation ( 3)); iiq A user requirement denoting the final accuracy wanted for the outputs. Hereby, the term accuracy refers to the number of significant bits required by the user on a variable of the program, denoted by nsb. Formally, let x be the approximation of x in finite precision and let εpxq " |x ´x| be the absolute error. So, if nsbpxq " k, for x " 0, then we have

εpxq ď 2 ufppxq´k`1 . ( 10
)
An ILP problem can be generated from the program source code which can be optimally solved by a LP solver (we use GLPK 6 , in practice). Concerning our resulting data types, the key feature of our method consists in finding directly the minimal number of bits needed at each control point of the original program. Next, these precisions can be approximated to the upper number of bits corresponding to an existing format int16 t, int32 t, etc. By way of illustration, if a variable x has nsbpxq " 18 bits, then x is tuned to the int32 t format. After solving the ILP problem, POPiX collects a new key information concerning the optimized precisions (along with the ufp integer quantity already computed by range analysis) in order to fully specify the fixed-point formats xM, Ly introduced in Section 3.1. Finally, we call our fixed-point library to synthetize a fixed-point version of the program with only integer numbers. Through this technique, it is possible to achieve memory savings up to more than double with a precision cost that depends on the original program being optimized and energy savings up to 3.5ˆ(see Section 4).

Implementation Details POPiX has been developed in JAVA and C++ and uses the ANTLR tool v4.7.17 to parse the input programs. To be able to guarantee that no overflows will occur in computations, we perform a range analysis by launching the execution of the program a certain number of times in order to determine dynamically an under-approximation of the range of variables by using the ufp of the values. In the future, we plan to use a static analyzer. Nevertheless, POPiX uses the simple imperative language below.

x P Id ℓ P Lab d P {+, -, ˆ, ˜} math P {sin, cos, tan, arcsin, log, . . .} Expr Q e : e :

:= c#p | x | e ℓ 1 1 d e ℓ 2 2 | mathpe ℓ 1 q | sqrtpe ℓ 1 q Cmd Q c : c ::" c ℓ 1 1 ; c ℓ 2 2 | x " e ℓ 1 | while b ℓ 0 do c ℓ 1 1 | if b ℓ 0 then c ℓ 1 1 else c ℓ 2 | create vectorpv, sq | create matrixpm, r, cq | require nsbpx, nq
We denote by Id the set of identifiers and by Lab the set of control points of the program used to assign to each element e P Expr and c P Cmd a unique control point ℓ P Lab. Fortunately, POPiX is able to handle loops, conditionals and arrays. The declaration of vectors is expressed by the statement create vector(v,s), while s denotes the size of the vector v. The declaration of a matrix m is expressed by the statement create matrix pm, r, cq , while r and c denote respectively the number of rows and columns of the matrix. The statement require nsb(x,n) indicates the minimal nsb n that a variable x must have at a control point. The rest of the grammar is standard. Note that the usual mathematical elementary functions are supported.

Cost functions Cost functions are given as optimization objective to the linear solver. Depending on which cost function is used, different criteria may be considered for the tuning phase. POPiX currently handles the following cost functions: CF1 Optimizes the sum of the number of significant bits of all the variables at each control point of the program. CF2 Optimizes the sum of accuracies of only the variables assigned in the program. Compared to CF1, this cost function minimizes the size of the variables and the number of bits needed for the operators to store the intermediary results. Let us note that CF2 is used in the experiments of Section 4. CF3 Minimizes the maximal accuracy needed in the program, i.e. the worst accuracy at some control point of the tuned program. This function is usefull to make a program fit in a certain format (for example all variables in 16 or 32 bits.) CF4 Optimizes only the sum of the accuracies of the arithmetic operators and elementary functions. This function is relevant from an hardware point of view, for example to limit the size of the operators in FPGAs [START_REF] Gao | Numerical program optimization for high-level synthesis[END_REF]. CF5 Minimizes the number of type conversions. Indeed, type conversions introduced by the mixed-precision tuning may slow down the execution of the programs and one may prefer a compromise between memory savings and execution time. This function addresses this problem. Note that these cost functions are modified when dealing with arrays: the tool multiplies the precision by the number of elements and this process is done only once for each array instead of several times for each use of arrays.

Experimental Evaluation

In this section, we conduct some experiments to show the effectiveness of our code synthesis method presented in Section 3.

Experimental Metrics

In our experiments, our goal is to evaluate the benefits of POPiX in terms of mixed-precision, memory savings and energy consumption which are important metrics to validate our synthesis method. We also measure the time of analysis spent by POPiX and the execution time of the fixed-point progam with respect to the floating-point program in which we assume that all variables are in single precision (32 bits) before the analysis. For our benchmarks, we use applications from FPBench 8 , a synthetic benchmark for floating-point performance. We run each program with three accuracy requirements arbitrarily chosen by the user: 4, 8 and 16 bits which bound the relative error of the result. All the results we report in this section where gathered on two machines: Ubuntu 20.04 LTS, with an 2.7GHz i7 core and 16 GB of RAM and Ubuntu 20.04 LTS, with a CPU AMD Ryzen 5 3500u and 5.7 GB of RAM. Let us state that the reason we use the Intel machine is to exploit the Jouleit 9 tool in order to estimate the power consumption of the CPU, RAM and integrated GPU.

Results Analysis Table 1 shows the mixed-precision configurations obtained after analysis in terms of number of variables or operations that we may tune into int8 t, int16 t and int32 t and consequently the memory savings in terms of number of bits. The second left-most column of Table 1 headed "call" refers to the number of elementary functions in the code and the next column headed "op" denotes the number of elementary operations. In this experiment, we assume that 100% is the percentage of all variables initially in single precision. Clearly, the memory savings compared to the initial number of bits for the majority of the original programs is considerable reaching 75% for "CRadius" program for a requirement of 4 bits. For instance, the "carbonGas" program has a total of 13 variables all in single precision before analysis. In the synthesized code we obtain 7 variables tuned into int8 t and 6 variables in int16 t achieving a gain in number of bits of 63,5%. Concerning "azimuth"program, our analysis failed to infer mixed precision for a user accuracy requirement of 16 bits whereas it 4 depicts the energy consumed by the execution of the benchmarks for a requirement of 4 bits. We observe that the fixed-point version codes require significantly less energy than the floating-point codes. For instance, the energy saved on CPU and DRAM reaches « 43% for "carbonGas" program and more than 75% for "jetEngine" program. Let us note that this observation is also valid for the remaining user requirements with slight variations of savings. Finally, we present the results in terms of speed for each of our benchmarks in Table 2. We denote respectively by "t f loat ","t synthesis " and "t f ix " the time of execution of floating-point programs, the total synthesis time of POPiX and the execution time of fixed-point programs all given in milliseconds. We visualize that the time spent by POPiX for the majority of benchmarks is negligible not exceeding 342 ms for the "carbonGas" program (« 30 LOCs). Although our synthesis method is fast (few seconds), we observe that the execution time remains the same for the floating and fixed-point codes with negilible slow-down for some benchmarks.

Related Work

In recent years, many authors have investigated the possibility of automating fixed-point code synthesis. A similar approach to our work is the TAFFO tool proposed by Cherubin et al. [START_REF] Cherubin | Dynamic precision autotuning with TAFFO[END_REF]. TAFFO is a static precision tuning tool that Another solution for code conversion was introduced by Cattaneo et al. [START_REF] Cattaneo | Embedded operating system optimization through floating to fixed point compiler transformation[END_REF]. Their method relies on a self-contained compiler transformation pass implemented within LLVM to perform the conversion. Their tool was especially dedicated to MIOSIX, a real time operating system targeting embedded system.

The goal of the dissertation of Jha [START_REF] Jha | Towards Automated System Synthesis Using SCIDUCTION[END_REF] is to give an algorithm for optimal fixed-point expressions synthesis based on inductive synthesis. Two years later, Darulova et al. [START_REF] Darulova | Synthesis of fixed-point programs[END_REF] proposed a fixed-point program synthesis methodology based on expression rewriting and genetic programming. Their algorithm uses abstract interpretation to estimate the error bound of a fixed-point implementation. However, the latter two techniques provide pessimistic bounds for non-linear expressions and are limited to straight-line programs. Coversely, Aslan et al. [START_REF] Semih Aslan | A high-level synthesis and verification tool for fixed to floating point conversion[END_REF] developed a tool that takes an n-bit fixed-point input and creates an m-bit floating-point output with IEEE754 and custom formats.

In the context of polynomials, linear filters and signal processing algorithms, the members of the DEFIS project [13] presented many approaches for fixedpoint code synthesis. To mention a few, the idea described in [START_REF] Deest | Toward scalable source level accuracy analysis for floating-point to fixed-point conversion[END_REF] works by infer- ring high-level convolution operations from the original source code, and modeling them as part of the program representation. In addition, Najahi et al. [START_REF] Martel | Trade-offs of certified fixedpoint code synthesis for linear algebra basic blocks[END_REF] presented an automated approach to synthesize codes in fixed-point arithmetic for some linear algebra basic blocks. They take a mathematical description of the problem as well as the range of the input variables and generate fixed-point code. Lopez [11] addresses the transformation of linear filters and controllers into hardware operators using fixed-point arithmetic. His main contribution is a complete error analysis, with respect to the internal word-lengths and the formulation of the word-length optimization as a convex non-linear integer optimization problem solved using appropriate heuristics. An extension of this work to the full class of linear time invariant algorithms has been proposed in [START_REF] Volkova | Towards reliable implementation of digital filters[END_REF].

Conclusion

In this article, we have presented a new method for fixed-point code synthesis respecting an accuracy requirement imposed by the user. Our method is based on a static analysis of the code implemented by means of system of constraints which gives the minimal format needed to encode each value. Experimental results show the performance of the codes synthesized in terms of execution time, memory and energy savings on a set of benchmark related to embedded systems.

In the future we would like to validate our method by considering architectures more commonly used in embedded systems. Also, we aim at generating hardware instead of software fixed-point implementations, using FPGAs. Targeting FPGAs has two justifications: this type of hardware is becoming more and more popular today and it presents the advantage of allowing fully custom designs. Finally, for adoption reasons in real-world applications, we aim at extending POPiX in order to handle full C programs via an integration to LLVM.

Fig. 1 :

 1 Fig. 1: Left: tuned program generated with a pair |ufp, nsb| for each variable as highlighted in blue. Right: Program C generated with fixed point formats.

Fig. 3 :

 3 Fig. 3: POPiX workflow.

Fig. 4 :

 4 Fig. 4: Measurement of the energy consumption (CPU and DRAM) of the floating-point and fixed-point version of our benchmarks.

Table 1 :

 1 Mixed fixed-point formats in 4, 8 and 16 bits for the synthesized program with the precentages of the number of bits saved. reaches 4.2% and 2.8% respectively for requirements of 4 bits and 8 bits. Another observation is that for a requirement of only 4 bits, 17 variables are tuned into int32 t. A possible explanation of this result is the call to the elementary functions in the code (call = 7) which can use intermediate variables that request greater precision than the user accuracy requirement. Figure

	8 https://fpbench.org/
	9 https://github.com/powerapi-ng/jouleit

Table 2 :

 2 Execution time measurements obtained during the experiments.

			4 bits		8 bits		16 bits
	Program	t f loat	t synthesis	t f ix	t synthesis	t f ix	t synthesis	t f ix
	azimuth	1.91	340	1.8	301	2.6	233	3.2
	carbonGas	0.17	342	0.29	233	0.31	201	0.46
	CRadius	0.10	240	1.31	192	1.33	186	1.36
	CTheta	0.38	203	0.56	223	0.57	273	0.57
	doppler1	0.16	205	0.53	249	0.54	225	0.57
	doppler2	0.24	188	0.28	225	0.29	207	0.30
	doppler3	0.17	243	0.30	195	0.32	207	0.34
	instantCurrent	1.01	264	1.99	265	2.36	273	2.73
	jetEngine	0.34	294	1.18	264	1.24	276	1.82
	LeadLagSystem	0.53	352	0.71	394	0.86	346	1.09
	CX	0.43	178	0.29	208	0.37	175	0.52
	CY	0.39	203	0.41	197	0.55	197	0.68
	triangle12	0.17	199	1.57	208	1.63	203	2.53
	turbine1	0.24	220	0.31	228	0.35	269	0.49
	turbine2	0.28	246	0.49	222	0.61	212	1.01
	turbine3	0.48	247	0.47	253	0.49	217	0.83

https://code.google.com/archive/p/libfixmath/

http://savannah.nongnu.org/projects/fixmath/

https://github.com/MikeLankamp/fpm

https://www.gnu.org/software/glpk/

https://www.antlr.org/

‹ This work is supported by La Région Occitanie under Grant GRAINE -SYFI: https://www.laregion.fr.