Frédéric Loulergue
email: frederic.loulergue@nau.edu

Jolan Philippe
email: jolan.philippe@imt-atlantique.fr

Automatic Optimization of Python Skeletal Parallel Programs

Keywords: Algorithmic skeletons, Program transformation, Python

Skeletal parallelism is a model of parallelism where parallel constructs are provided to the programmer as usual patterns of parallel algorithms. Highlevel skeleton libraries often offer a global view of programs instead of the common Single Program Multiple Data view in parallel programming. A program is written as a sequential program but operates on parallel data structures. Most of the time, skeletons on a parallel data structure have counterparts on a sequential data structure. For example, the map function that applies a given function to all the elements of a sequential collection (e.g., a list) has a map skeleton counterpart that applies a sequential function to all the elements of a distributed collection. Two of the challenges a programmer faces when using a skeleton library that provides a wide variety of skeletons are: which are the skeletons to use, and how to compose them? These design decisions may have a large impact on the performance of the parallel programs. However, skeletons, especially when they do not mutate the data structure they operate on, but are rather implemented as pure functions, possess algebraic properties that allow to transform compositions of skeletons into more efficient compositions of skeletons. In this paper, we present such an automatic transformation framework for the Python skeleton library PySke and evaluate it on several example applications.

Introduction

Context and Motivation Most computing devices are now parallel architectures, and more and more data is produced and analyzed. However, writing programs for parallel architectures remains difficult compared to writing sequential programs. Parallel programming may increase programs performances but also hinders programming productivity.

The skeletal approach proposed by Murray Cole [START_REF] Cole | Algorithmic Skeletons: Structured Management of Parallel Computation[END_REF] consists in defining a program using computational patterns that have parallel implementations. In other words, programmers do not have to think about parallel aspects of their program, but only about how to write their program using already implemented patterns. Skeletons are often implemented as higher-order functions, providing both flexibility and expressivity for the developers. Generally, a skeleton on a distributed datastructure is associated with a sequential function for a corresponding sequential datastructure. Such a high-level approach makes parallelism more affordable, even for nonexpert parallel developers. However, most of the libraries handling parallel data-structures have been developed in order to be used by advanced users: the use of low-level languages or low-level primitives for parallelism, etc. For example, SkeTo [START_REF] Emoto | An Automatic Fusion Mechanism for Variable-Length List Skeletons in SkeTo[END_REF], SkePu [START_REF] Enmyren | SkePU: A Multi-Backend Skeleton Programming Library for Multi-GPU Systems[END_REF], OSL [18], or Muesli [START_REF] Ciechanowicz | The Münster Skeleton Library Muesli -A Comprenhensive Overview[END_REF] are skeleton libraries written in C++. Their main advantage is performance. However, C++ itself is a complex language compared to Python and its dynamic typing and high-level features.

PySke [START_REF] Philippe | PySke: Algorithmic skeletons for Python[END_REF], whose name is composed by Python and Skeleton, aims at providing computational patterns for lists and trees in both sequential and parallel, making the parallel aspects of programs as abstract as possible. PySke targets users without a deep background in computer science. There are two main challenges for such developers when using a skeleton library that provides a wide variety of skeletons: the choice of skeletons to compose the application, and the way these skeletons are composed. These application design choices may have a big impact on the performance of the application. Nonetheless, compositions of skeletons, especially on immutable data-structures, obey some algebraic laws that can be used for optimizing skeleton compositions.

Contributions In [START_REF] Philippe | Towards automatically optimizing PySke programs (poster)[END_REF] we evaluated such optimizations. However, the transformations were not automated but hand-written. In this paper, we propose a fully automatic optimizer of PySke programs based on program transformation, and we evaluate it on some example applications.

Outline The paper is organized as follows. In Section 2, we discuss work related to our contribution. We give an overview of PySke in Section 3. Section 4 is devoted to the design of the new program transformation feature of PySke. We evaluate the performance of optimized example programs in Section 5. We conclude and discuss future work in Section 6.

Related Work

The transformation rules we use to optimize programs come from the Bird-Merteens Formalism (BMF) tradition [START_REF] Bird | Algebra of Programming[END_REF][START_REF] Bird | An introduction to the theory of lists[END_REF]. For example, the map function/skeleton -that applies the same function to all the elements of a list -satisfies the following property: The composition of two map can be transformed into a single use of map. Considering two functions f and g, (map f ) • (map g) can be transformed into map (f • g), where • is function composition. This transformation removes the allocation and transversal of an intermediate list, thus optimizes the initial program.

To our knowledge, PySke is the only algorithmic skeleton library for Python. But there are many skeleton libraries or eDSL for other programming languages as well as programming languages specifically designed around the concept of algorithmic skeletons. In this section, we discuss only libraries and languages where program transformation is used as a means to optimize parallel programs.

For skeletons libraries, there are two main categories:

-Libraries where the optimization is done at compile time using meta-programming techniques either provided directly by the host language, or extensions, or abusing some other language features;

-Libraries where the optimization is done at execution time, either by encoding the skeleton expressions in a data structure similar to an abstract syntax tree and perform the optimization on this structure, or by using the introspection/reflection features of the host language.

In the first category, the Delite [START_REF] Sujeeth | Delite: A compiler architecture for performance-oriented embedded domainspecific languages[END_REF] framework for Scala can be considered as a skeletal parallelism approach. It features a set of data structures and mostly classical skeletons on them: map and variants, reduce and variants, filter, sort; and one less usual skeleton: group-by, as Delite has dictionaries as one of its supported data structures. Delite provides compile-time optimization through staged programming. Delite targets heterogeneous architectures CPU/GPU but only shared memory architectures. Several C++ libraries [START_REF] Emoto | An Automatic Fusion Mechanism for Variable-Length List Skeletons in SkeTo[END_REF][START_REF] Falcou | Quaff: Efficient C++ Design for Parallel Skeletons[END_REF]18] are also in this category: they abuse the template feature of C++ to provide compile-time optimization. Note that in these libraries, the transformations are quite ad-hoc, and the main goal is to remove intermediate data-structures and loop traversals. The map transformation rule given above is actually such a transformation. But more complex or high-level transformations are usually not handled by these libraries that often only optimize the sequential parts of the programs but stop their optimizations as soon as there are communications. To our knowledge, there is no compile-time staging framework for Python thus such an approach cannot be used for PySke without writing a specific compiler (possibly a source-to-source compiler).

In the second category, Accelerate [START_REF] Chakravarty | Accelerating Haskell array codes with multicore GPUs[END_REF][START_REF] Clifton-Everest | Streaming irregular arrays[END_REF] -a Haskell framework for programming GPUs -features classical data-parallel skeletons (map and variants, reduce, scan and permutation skeletons) on multi-dimensional arrays and streams. Optimization of the GPUs programs are done at runtime. PySke also belongs to this category: we optimize compositions of skeletons at runtime, and we rely on a data-structure representing these compositions to perform the transformations. In this category, another approach is to rely also on the reflection features of the host language. Lithium [START_REF] Aldinucci | An Advanced Environment Supporting Structured Parallel Programming in Java[END_REF] is a structured parallelism framework that leverages Java reflection to optimize programs. Introspection and reflection also exist in Python, and we do rely on such features in the extension of PySke considered in this paper.

Finally, there are skeleton languages. The main advantage of such an approach is that the way optimizations are performed is not limited to what the host language allows. The main drawbacks are that specific compilers should be designed and implemented, and new languages do not have as large libraries as mainstream languages do, which hinders their adoption. An intermediate approach is to have a specific language for writing the composition of skeletons and use a usual sequential language for writing the arguments to these skeletons. This is less flexible than having a completely new language but mitigates the drawbacks mentioned. P3L [START_REF] Pelagatti | Structured Development of Parallel Programs[END_REF] associated with the FAN transformation framework [START_REF] Aldinucci | Towards parallel programming by transformation: the FAN skeleton framework[END_REF][START_REF] Bacci | Skeletons and transformations in an integrated parallel programming environment[END_REF] is a representative of this approach. While being very efficient, this kind of approach is much more complex to deploy than a Python library and requires a significant effort from the programmer.

An Overview of PySke

PySke is a library for Python currently implemented on top of MPI and mpi4py [START_REF] Dalcin | Parallel distributed computing using Python[END_REF]. Note that the programming model of PySke is independent of the underlying communication library.

PySke offers a global view of programs. A PySke program is written (and read) as a sequential program but it operates on parallel data structures. This aspect is very different from the SPMD paradigm of MPI where most of the time a program is actually parametrized by the process identifier (returned by the method Get_rank in mpi4py), and the global parallel program should be understood as the parallel composition of instantiations -for all possible process identifiers -of this parametrized sequential program. This "par of seq" structure is more complicated to deal with than the "seq of par" structure that global view offers [8].

For readers familiar with MPI, PySke can be thought as a library of collectives. There is a major difference: arguments to MPI collectives have regular C types but the collection on all processors of the values of these sequential types may be thought as a parallel data structure whereas in PySke there are classes dedicated to parallel data structures. In PySke, the type of a value indicated whereas this value is sequential or parallel. In MPI there is no parallel type, and it may be difficult (and it is an undecidable problem in general) to know if a value is sequential or should be thought as being part of a distributed data structure. When possible, PySke offers the same methods for both a sequential data structure and the corresponding parallel data structure. The code in Figure 1 computes the variance of a discrete random variable data. This code is valid when data is either an instance of SList, i.e. a sequential list, or an instance of PList, i.e. a parallel list distributed on all the processors of the parallel machine running the program.

This example shows two classical skeletons: map that applies the same function (here f) to all the elements of a data-structure, and reduce that uses a binary associative operation (here add as defined in the Python module operator) to "sum" all the elements of a list using the binary operation. For their implementations in PList, map does not require any communication to be executed and reduce does require some communication: first, the partial sums are computed on each processor, then these partial sums are sent to all processors (total exchanged), and finally, the final sum is computed.

The user of PySke can see a PList as a usual list (SList extends the Python lists with additional methods). The implementation of PList is however an MPI one and therefore follows the SPMD approach. On each processor, a PList is actually composed of several fields: the content (a sequential list specific to each processor), the global size (same value on all processors), the distribution (the same list on each processor, this list contains the lengths of all local contents), and the index of the first element of the local list in the global list (for example at processor 2, this value would be 9 if processor 0 has a local list of 5 elements, and processor 1 a local list of 4 elements). All the methods provided by PList ensure that the content of these fields stay what we have just described: This may require communications for some skeletons which intuitively do not need them to perform the intended computations on the content of the parallel list.

For example, the filter skeleton -that takes as argument a predicate p and returns a parallel list where only the values satisfying the predicate p are kept -does not require any communication to compute the content of the returned parallel list. It is enough to perform a sequential filter on each of the local contents. However, if the skeleton only performed this local filtering, the global size, the local index, and the distribution would no longer represent the actual global size, local index, and distribution.

The current set of list skeletons contains [START_REF] Philippe | PySke: Algorithmic skeletons for Python[END_REF]: reduce, map and variants (mapi, map2, zip), scan and variants, get_partition, flatten, and balance.

For example, SList([1,2,3]).scan(add) is [0, [START_REF] Aldinucci | An Advanced Environment Supporting Structured Parallel Programming in Java[END_REF][START_REF] Baader | Term Rewriting and All That[END_REF][START_REF] Bird | Algebra of Programming[END_REF]. get_partition makes the way the data-structure is distributed (or partitioned) visible in the value itself. For example if the global view of a parallel list pl is [1, 2, 3, 4, 5, 6] and the distribution on 4 processors is [2, 2, 1, 1], then the global view of pl.get_partition() is [ [START_REF] Aldinucci | An Advanced Environment Supporting Structured Parallel Programming in Java[END_REF][START_REF] Aldinucci | Towards parallel programming by transformation: the FAN skeleton framework[END_REF], [START_REF] Baader | Term Rewriting and All That[END_REF][START_REF] Bacci | Skeletons and transformations in an integrated parallel programming environment[END_REF], [START_REF] Balland | Effective strategic programming for Java developers[END_REF], [START_REF] Bird | Algebra of Programming[END_REF]]. flatten is the inverse operation (and requires communications).

These two skeletons can be used to implement a filter skeleton:

def filter(self, p): return self.get_partition(). map(lambda l: l.filter(p)). flatten()

After a call to filter the resulting parallel list may be unbalanced, i.e., some processors may contain much more elements than some others. The balance skeleton redistributes values such that the parallel list is evenly distributed, i.e., any processor has at most one more element than each other processor.

PySke also features a set of skeletons on trees. There are three tree data-structures: binary trees, linearized trees, parallel trees. The two first structures are sequential structures, while the last one is a parallel data structure. All of them represent binary trees.

Figure 2 is the code necessary to count the number of elements of a parallel tree that satisfy a predicate p. The code differs for a value of class BTree of sequential binary trees (variable bt) and a value of class PTree of parallel trees (variable pt). This is due to the fact that in order to be executed in parallel the reduction of a tree using an operator ⊕, this operator should satisfy a property called the closure property that allows expressing the operator using 4 auxiliary functions. This property is in a way what corresponds to associativity for lists. map on trees needs two arguments: a function to apply to the leaf values and a function to apply to the values at the nodes. There def id(x): return x def sum(x, y, z): return x + y + y def f(x): return 1 if p(x) else 0 count_bt = bt.map(f,f).reduce(sum) count_pt = pt.map(f,f).reduce(sum, id, sum, sum, sum)

PySke with Automatic Optimization

In order to support automatic optimization of PySke applications, the user programming interface needs to slightly change with respect to what we presented in the previous section. We present these changes in Section 4.1. The remaining sub-sections are devoted to the design of the automatic optimization mechanism of PySke which is rule-based. It therefore relies on concepts from term rewriting systems [START_REF] Baader | Term Rewriting and All That[END_REF]: we give a short overview in Section 4.2, and explain how we implemented support for such systems in a generic way in Python in Section 4.3. This generic framework is then instantiated to deal with PySke skeletons (Section 4.4).

User Programming Interface

We illustrate the differences between the previous API and the new one through the example of Figure 3: the computation of the dot product of two vectors represented as two parallel lists.

dot_product1 is the version written using only the features presented in Section 3. It uses the data structure PList (imported as PL in this example) of the module pyske.core.list.plist and its skeletons.

dot_product2 is the version written using the new API: the PList class of the module pyske.core.opt.list. First note that dot_product2 can still take as input usual parallel lists: we just need to use a wrapper that transforms a PL into a PList (Lines 13-14). If dot_product were to take as argument parallel lists of type PList then no wrapping would be necessary. The only other difference with the previous API is that it is necessary to call the method run to launch the optimization of the skeleton composition and the execution of the optimized version. Note that to implement dot_product2 we can reuse dot_product1: this makes the transition to the new version easier.

uncurry is a higher-order function that transforms a function taking as input two formal parameters into a function taking only one formal parameter, this parameter being a pair. curry does the inverse transformation. Compositions of such functions are optimized too. The performances of this example are discussed in Section 5.

Term Rewriting

Basically, PySke programs are skeleton expressions. Therefore, they can be represented as terms, and we can use rules to transform these terms. Our optimization mechanism is thus a term rewriting system [START_REF] Baader | Term Rewriting and All That[END_REF].

Terms A signature Σ is a set F of function symbols with a function ar : F → N that for each function symbol f gives its arity. We call symbol functions of arity 0 constants. We assume a countable set X of variables such that X ∩ F = ∅. The set T of terms over X and Σ is the smallest set such that:

for all variable x ∈ X , x ∈ T , for all function symbol f ∈ F, number n such that n = ar (f ), and for all terms t 1 , . . . , t n then f (t 1 , . . . , t n ) ∈ T .

The set of variables of a term t can be computed by the following function V:

V(x) = {x} V(f (t 1 , . . . , t n )) = ∪ n k=1 V(t k ) A closed
term is a term t that does not contain any variable, i.e. V(t) = ∅. A term is said linear if each x ∈ V(t) appears at most once in t.

A substitution σ is a partial function from X to T . It can be extended to T as follows: for any t ∈ T , either t ∈ X and we can directly use σ, or there exist a function symbol f of arity n, and n terms t 1 , . . . , t n such that t = f (t 1 , . . . , t n ). We can then define: σ(t) = f (σ(t 1 ), . . . , σ(t n )).

σ ⊥ denotes the substitution undefined everywhere, and for a substitution σ, a variable x and a term t, the update σ[x → t] is the substitution defined by:

σ[x → t](y) = σ(y) if y = x t if y = x
Rules A rule is a pair of terms, denoted by l → r, such that V(r) ⊆ V(l). A rule is linear if both l and r are linear. Assuming t is a closed term and p is a linear term, we define the function match(t, p) as follows:

       match(t, x) = σ ⊥ [x → t] match(f (t 1 , . . . , t n ), g(t 1 , . . . , t m )) = merge(match(t 1 , t 1 ), . . . , match(t n , t n )) if f = g and n = m match(f (t 1 , . . . , t n ), g(t 1 , . . . , t m )) = ⊥ otherwise
The function merge returns ⊥ if one of it arguments is ⊥, otherwise it merges the substitutions in arguments. As p is supposed to be linear, each substitution maps variables that are different from the variables mapped by each other substitution.

The application of a rule l → r at the root of a closed term t proceeds as follows:

1. we check if the left-hand side of the rule matches t, i.e. we compute σ = match(t, l), 2. if σ = ⊥ then it means l does not match t, so the rule cannot be applied, 3. if σ is a valid substitution, the rule can be applied and then we obtain the term σ(r).

A rule can be applied to sub-terms of a term t. In this case, if the matching succeeds, it is the matched sub-term that is replaced by σ(r).

Example Let us consider the following signature: F = {map, reduce, pl, incr, +, •} with ar (map) = 2, ar (reduce) = 2, ar (pl) = 0, ar (incr) = 1, ar (+) = 2, ar (•) = 2, and we use + and • as infix operations. They represent respectively the addition and function composition.

Let us consider the following term: t = reduce(map(map(pl, incr), incr), +)

that can be thought as a representation of the following PySke program:

pl.map(incr).map(incr).reduce(add)

and the following rule:

map(map(l, f ), g) → map(l, g • f )
where l, f , and g are variables.

The rule cannot be applied to the root of t, but the left-hand side of this rule matches the sub-term t = map(map(pl, incr), incr) of t:

match(map(map(pl, incr), incr), map(map(l, f ), g)) = [l → pl, f → incr, g → incr]
Applying the obtained substitution to the right-hand side of the rule gives:

map(pl, incr • incr)
and replacing t in t by this new term, we finally obtain:

reduce(map(pl, incr • incr), +)

This skeleton expression is likely to be more efficient because it avoids to create an intermediate list in-between the two calls to map. We use such rules in PySke to optimize skeleton compositions.

Term Rewriting in Python

Our framework first provides features that are close to the theoretical term rewriting system framework presented in the previous sub-section. These are provided in one Python module: terms.

The terms module offers two main classes: Var that just extends str and that corresponds to variables in the theoretical framework, and Term than implements nonvariable terms. This latter class has two main attributes: function and arguments that basically correspond to the function symbol and arguments in a theoretical term.

We however take advantage of the flexibility of Python type system. The function attribute can be either a string, or a Python function. The arguments are represented as a list. Its elements can be either values of type Term, strings (in this case representing a class), or Python values of arbitrary types. The two latter cases can be thought as constants in the theoretical framework.

Substitutions in the theoretical framework are implemented as Python dictionaries. Applying a substitution to a term t (in Python a value either of type Var or of type Term) is implemented as a Python function subst(t, s) where s is a dictionary.

One important feature of the theoretical framework is the match function. It is implemented in Python as a method of the class Term that takes as argument a value either of type Var or of type Term. This function either returns a dictionary in case the matching succeeds, or None if the matching fails (None thus corresponds to ⊥ of the theoretical framework).

Finally rules are implemented as a Python namedtuple, so essentially records. These records contain of course two terms, left and right, but also a name for the rule, and a type. This type is actually a sub-class of Term as we will explain in Section 4.4.

The module terms contains also a function to apply a rule at the root of a term, and a function inner_most_strategy(t) that repeatedly tries to apply all the available rules to every sub-terms of t until it is no longer possible to apply any rule. The strategy to apply the rules is discussed in the next section.

As an example, Figure 4 presents the encoding of the rule presented in the previous section. This example illustrates the fact that the function attribute of a term can be either a string (most of the cases) or a Python function. This is the case for compose (Line 4) which is a function defined as:

def compose(f, g): return lambda x: f(g(x))
Similarly the elements of arguments, the second parameter to the constructor of Term can be either variables, terms, or arbitrary Python values (this latter case is not present in the example).

Rule-Based Skeleton Compositions Optimization

We have a framework to represent generic terms, and rules to transform such terms.

There remain two main questions to use such a framework to automatically optimize skeleton compositions:

Rule(left=Term('map', [Term('map', [Var('PL'), Var('f')]), Var('g')]), right=Term('map', [Var('PL'), Term(compose, [Var('f'), Var('g')])]), name="map map", type=_List) Fig. 4: A Rewriting Rule in Python 1. How to represent PySke skeleton compositions as instances of Term in a way that is mostly transparent to the user with respect to the former PySke API? 2. How to evaluate/execute such terms?

The design principles for a solution to (1) are:

inherit from Term, and use Python introspection capabilities to obtain a concise implementation.

Thus essentially, for each class that provided a data structure and associated skeletons in PySke we define a "wrapper" class that features the same methods than the initial class, but these methods build terms instead of executing parallel code.

A class such as Plist has static methods (and we identify the class constructor to a static method in our framework) and methods. For both kind of methods, the representation is such that:

the method name is represented as the function symbol of the term (in class Term, it is attribute function), the class or object that is the target of the method call is represented as the first argument of the term (first element of the list arguments), the remaining of the method arguments are represented as additional elements in the list arguments.

Python offers many introspection features. In particular it is possible to capture any call to non-static methods that are not defined in a class. Therefore to define a PList wrapper class that inherits from Term we have basically to define static methods that correspond to the initial class static methods, and a generic way to catch calls to nonstatic methods to build the corresponding terms.

The initial PList class features the init(f, size) static method. The wrapper PList(Term) class thus features the following static method: @staticmethod def init(f, size):

return PList('init', ['PList', f, size])

As indicated before, the value for the attribute function is a string that is the name of the corresponding static method in the initial PySke PList class, and the arguments list, in addition to the arguments to init that are f and size, contains as first element the name of the class (here PList) that contains the init static method. As explained before, for non-static methods, we do not re-implement each method but rather rely on a feature of Python that allows to capture any call to non-defined methods in a class (technically we redefine the __getattr__ attribute). We refer to the source code 3 for more technical details. An instance of this generic code for method map follows:

def map(self, f): return PList('map', [self, f])
Using this wrapper class PList is therefore exactly the same than using the initial PList but nothing is executed: instead a term representing the skeleton expression is built. The only difference are: the raw method that allows to embedded a value of the initial class into the wrapper class as illustrated in Section 4.1, and launching the optimization and execution of such a term.

Indeed, in fine it is necessary to execute the skeleton composition. This is done using a method called run(). It proceeds in two steps:

1. first the term is transformed using rules such as the rule presented in Figure 4, 2. then the optimized term is executed.

The optimization and execution features are not specifically implemented for each wrapper class. They are actually offered by two methods of the class Term:

opt() transforms a term using all the available rules by calling the transformation function inner_most_strategy that applies all the available rules as many times as possible on all the sub-terms of the current term; eval() executes a term. We refer to the code for details, but basically Python allows to retrieve attributes of a class (including methods) by their names represented by a string, and to apply any function or method to a list of values in such a way that this list is considered as the effective parameters of the function or method. We use both this features to retrieve the target class and methods in the initial PySke classes from objects of the wrapper classes.

Strategy The transformation rules are applied using an innermost strategy: the rules are first applied (when possible) to the leaves of the tree representing the composition of skeletons, and then in an upwards fashion. When two rules may be applied to the same sub-expression of the skeleton expression, priorities associated to rules are used.

As future work we plan to extend this default behavior as discussed in Section 6.

Applications and Experiments

Currently, our framework contains a dozen of optimization rules:

optimization of a composition of maps, optimization of compositions of map and reduce (using and internal map_reduce skeleton), optimization of compositions of map and zip to map2, optimization of compositions of map and reduce using boolean operators as their arguments (a generalized De Morgan rule), optimization of compositions of higher-order functions such as uncurry and curry (for example, the composition of these two functions is the identity function, and the identity function is a neutral element for function composition).

We evaluated our framework on three examples:

1. the dot product presented in Figure 3, 2. the conjunction of the negation of a list l of Boolean values:

l.map(operator.not_).reduce(operator.and_, True)

3. the computation of the normalized average of vectors: a parallel list of vectors (with 10 components) is first normalized (the norm of each vector is computed and the vector is multiplied by a scalar that is the inverse of its norm), then the average of these lists is computed:

smul(1 / l.length(), l.map(normalize).reduce(vadd, vzero))

The experiments were conducted on a shared memory machine (256 Gb), with two Intel Xeon E5-2683 v4 processors each having 16 cores at 2.10 GHz. Each example has been run 30 times. We used the following software: CentOS 7, Python 3.6.3, mpi4py 3.0.2, OpenMPI 2.6.4. For the dot product example, we run four versions, on a parallel list of 5 × 10 7 elements:

direct corresponds to dot_product1, wrapper builds the term representation of the skeleton composition, then executes it without transformation. This version allows to check the overhead introduced by the building of the intermediate representation and the overhead introduced by using such a representation for execution, optimized is the same than wrapper but the skeleton composition is optimized using all available rules, finally hand_written is the version using directly the PList but the composition is optimized manually.

The results are presented Figure 5. First of all, there is no noticeable overhead due to the wrapper classes with respect to the use of the direct implementation. That is true both for the non-optimized code (direct and wrapper) and the optimized code (optimized and hand-optimized). Secondly, the automatic optimization of PySke provides, in this case, a significant speedup of about 2.5.

For the list of Boolean values example, the average speed-up of the automatically optimized version with respect to the direct implementation is 85%. For the vector average after normalization, the average speed-up is more modest at 16%. In all cases, the speed-up is there for only a slight change in the program with respect to the former API. q q q q q q 0 2 4 6 0 10 20 30

Number of processors

Average time (s) type q direct hand-written optimized wrapper q q q q q q 1.0 Skeletal parallel programming is a productive way to write parallel programs. With PySke we propose a library simple to use for non-expert programmers. However even with algorithmic skeletons, two difficulties remain: the choice of skeletons when the library is large, and the way these skeletons are composed to form an application. These choices have a large impact on the performance of the parallel program. To further eases the development of efficient parallel programs with PySke, we propose automatic program transformation based on high-level transformation rules. Evaluation of this new mechanism in PySke shows it can provide significant improvements to skeletal programs.

For the moment the transformation rules are hand-written by the developers of PySke. We plan to have a more flexible approach by extending our framework so it can read transformation rules from a file. This would allow expert users of PySke to add more domain specific transformation rules. The transformation rules will be written in a format similar to what can be found in term rewriting tools such as ELAN [START_REF] Van Den Brand | The ELAN environment: a rewriting logic environment based on ASF+SDF technology -system demonstration[END_REF]. Another extension is to provide ways to change the strategy of application of these rules. A very flexible solution would be to allow expert users to write their own strategies as it is possible in embedded term rewriting tools [START_REF] Balland | Effective strategic programming for Java developers[END_REF].

Fig. 1 :

 1 Fig. 1: Variance in PySke

Fig. 2 :Fig. 3 :

 23 Fig. 2: Count on Trees in PySke

Fig. 5 :

 5 Fig. 5: Dot Product Example

available at https://github.com/pyske/PySke


Frédéric Loulergue
email: ‡frederic.loulergue@nau.edu

Jolan Philippe
email: jolan.philippe@imt-atlantique.fr

New List Skeletons for the Python Skeleton Library

Keywords: high-level parallel programming, structured parallelism, algorithmic skeleton, Python, lists

Algorithmic skeletons are patterns of parallel computations. Skeletal parallel programming eases parallel programming: a program is merely a composition of such patterns. Dataparallel skeletons operate on parallel data-structures that have often sequential counterparts. In algorithmic skeleton approaches that offer a global view of programs, a parallel program has therefore a structure similar to a sequential program but operates on parallel data-structures. PySke is such an algorithmic skeleton library for Python to program shared or distributed memory parallel architectures in a simple way. This paper presents an extension to PySke: new algorithmic skeletons on parallel lists. This extension is evaluated on an application.

I. INTRODUCTION

Parallel architectures are now everywhere, from smartphones to super-computers. But while parallel architectures are mainstream, it is not the case for parallel programming and skilled parallel programmers. Python is a language more and more popular: in popularity rankings (based on internet search queries) and in lines of code in public repositories it is ranked as one of the first three programming languages in 2019. It is also more and more popular as an introductory teaching language, both in computer science curricula [START_REF] Guo | Python is now the most popular introductory teaching language at top u.s. universities[END_REF] but also in informatics curricula. A recent survey of the Python Software Foundation and JetBrains 1 also shows that in 2018 58% of developers used Python for data analysis (from 50% in 2017). This ratio is even higher for the 84% of developers who use Python as their main programming language.

As the size of data to analyze grows, it is necessary to provide the programmers who are not expert in parallel programming, and often not computer scientists, ways to leverage the parallel architectures they have for data analysis in Python.

In the various approaches for high-level parallel programming, algorithmic skeletons [START_REF] Cole | Algorithmic Skeletons: Structured Management of Parallel Computation[END_REF], [START_REF] Pelagatti | Structured Development of Parallel Programs[END_REF], [START_REF] González-Vélez | A survey of algorithmic skeleton frameworks: high-level structured parallel programming enablers[END_REF] are quite popular even though the most popular frameworks belonging to this family of approaches are not presented as skeletal parallelism frameworks. Basically, an algorithmic skeleton is a pattern of parallel computation. There are task-and data-parallel skeletons. Often data-parallel skeletons operate on parallel data-structures that have sequential counterparts, and therefore the skeletons themselves have sequential counterparts. For example the map skeleton is often a skeleton that applies a sequential function to all the elements of a parallel collection 1 https://www.jetbrains.com/research/python-developers-survey-2018 and it is also a function on sequential lists in functional programming languages. The most popular skeleton framework is Google's MapReduce [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF] although the skeletons it provides differ from more traditional map and reduce skeletons.

PySke [START_REF] Philippe | PySke: Algorithmic skeletons for Python[END_REF] is a library of algorithmic skeletons for Python, on parallel lists and parallel trees. PySke provides a global view of programs, i.e. a PySke program is written as a sequential program but it operates on parallel data-structures. It is simpler than the Single Program Multiple Data (SPMD) style of most MPI programs. When possible (it is not possible for all skeletons on trees), PySke provides the same methods for both a sequential and a parallel data-structure.

The contributions of the present paper are:

• new expressive algorithmic skeletons on lists,

• evaluation of the new skeletons on applications. The remaining of the paper is organized as follows. We discuss related work in Section V. Section II gives an overview of PySke both from the user and implementer perspectives. In Section III, we present new skeletons on the parallel list data-structure. We experiment with these new skeletons in Section IV and conclude in Section VI.

II. AN OVERVIEW OF PYSKE

A. A User's Perspective

PySke is a library for Python currently implemented on top of MPI and mpi4py [START_REF] Dalcin | Parallel distributed computing using Python[END_REF]. However, PySke's programming model is independent of the communication library.

PySke offers a global view of programs. A PySke program is written and read as a sequential program but operates on parallel data structures. This aspect is very different from the SPMD paradigm of MPI where most of the time a program is actually parametrized by the process identifier (returned by the method Get_rank in mpi4py). The global parallel MPI program should be understood as the parallel composition of instantiations -for all possible process identifiers -of this parametrized sequential program. This "par of seq" structure is more complicated to deal with than the "seq of par" structure that global view offers [START_REF] Bougé | The data parallel programming model: A semantic perspective[END_REF].

For readers familiar with MPI, PySke can be thought as a library of collectives. Nevertheless there is a major difference between MPI and PySke. Arguments to MPI collectives have regular C types but the collection all the values of these sequential types may be thought as a parallel data structure whereas in PySke there are classes dedicated to parallel data structures. In PySke, the sequential/parallel nature of values is def _max0_copy(num): return max(0, num), num def _max_sum(pair_a, pair_b): a_m, a_s = pair_a b_m, b_s = pair_b max_ = max(a_m, a_s + b_m) sum_ = a_s + b_s return max_, sum_ def mps(input_list: List): max_, _ = input_list .map(_max0_copy) .reduce(_max_sum, (0, 0)) return max_ Fig. 1. A First PySke Example therefore visible in their types. No parallel type exists in MPI, and it may be difficult (and it is an undecidable problem in general) to know whether a value is sequential or should be thought as being part of a distributed data structure.

PySke contains two main kinds of data structures: lists and trees. We refer to [START_REF] Philippe | PySke: Algorithmic skeletons for Python[END_REF] for a presentation of the tree data structures. In the present paper, we focus on list data structures. With respect to the first version of PySke, we know have an interface List that contains the documentation for all the methods on lists. It also features detailed type hints that we consider to be very useful documentation items. Both SList, the class of sequential lists, and PList, the class of parallel lists, implement this interface. In the remaining, a method on a sequential data structure is called a primitive, while it is called an algorithmic skeleton, or skeleton, on a parallel data structure. It is important to notice that PList are immutable data structures. All skeletons and primitives return new lists or scalar values. For each primitive we provide a skeleton, and vice versa. This allows for the rapid prototype of PySke applications, even without mpi4py installed. While SList is strictly speaking not an immutable data structure, because it inherits from list for convenience, the recommended style is to use it as an immutable data structure.

The program example of Figure 1 illustrates the polymorphic nature of PySke. The variable input_list can be either SList or PList. In the former case, the execution can be sequential only. In the later case, it uses MPI for parallel processing.

mps returns the largest sum of all the sums of the prefixes of the input list. For example, the prefix list returning the largest sum, is [6, -3, 0, 4] in the following execution2 : >>> mps(PList.from_seq([6, -3, 0, 4, -4])) 7 from_seq builds a parallel list from a sequential list. [3,2,2,2] [3,2,2,2] [3,2,2,2] [3,2,2,2]

>>> SList([1, 2, 3]).map(lambda x: x + 1) [2, 3, 4] >>> SList([1, 2, 3]).reduce(add, 0) 6 >>> SList([]).reduce(add, 0) 0
Basically, to compute the maximum prefix sum, we compute a pair containing the maximum prefix sum and the sum. The intuition behind the definition of _max0_copy is that if we consider an element as a singleton list. The maximum prefix sum for this list is either 0 if the element is negative (the prefix is then the empty list), or the element itself. _max_sum is defined as follows: pair_a represents the maximum prefix sum (a_m) and the sum (a_s) of a list a, and pair_b the same for a list b. To compute the maximum prefix sum of the list obtained by concatenating a and b, we need to take the maximum of the maximum prefix sum for a and the maximum prefix for a followed by b. The maximum prefix sum value for this second list is a_s + b_m.

There are several variants of the map skeleton:

1) mapi where the function argument to mapi is applied to each index and associated value in the list, 2) map2 that takes two arguments in addition to the current object: a function and a list, and applies the function to elements of both lists; these two lists should have the same distribution (as explained in the next section); map2i is a variant of map2, 3) zip creates a list of pairs from two lists having the same distribution.

There are several variants of the prefix sum skeleton scanl. An example of scanl follows: >>> par_list = PList.init(str, 4) >>> par_list.scanl(concat, '').to_seq() ['', '0', '01', '012'] init is the main way to build parallel list in PySke. It takes as argument a function used to initialize the values of the list, and the size of the list. to_seq() transforms a parallel list into a sequential list. par_list.to_seq() is ['0', '1', '2', '3'].

B. An Implementer's Perspective

Users can think about their PList in a global way, i.e. as a sequential list of type SList. However these lists can be indeed distributed on different machines, in a transparent machine. The underlying implementation of PySke follows the SPMD paradigm and is currently implemented using MPI.

For an example list, Figure 2 shows both the global view presented to the user of the library, and the way this view is realized following the SPMD paradigm (for 4 processors). Each row in the SPMD view corresponds to a private attribute in the Python implementation. Names in italic represent attributes that contain the same values on all processors. "content" is specific to each processor and contains the specific part of the global list that this given processor has. Note that the type of this attribute in SList. "start index" is the index in the global size of the first element of the local list. In the global list, the value 6 has index 3. But this value is the first element of the local list at processor 1. Therefore the value of "start index" at processor 1 is 3. The "distribution" attribute contains the list of all local sizes.

init builds an evenly distributed list, i.e. a list where the local sizes differ at most by 1, and the possibly larger local lists are on processors of low identifier. The list presented in Figure 2 is such an evenly distributed list. from_seq does not build an evenly distributed list: all the elements at are processor 0.

map, scanl and their variants return lists that have the same distribution as their input. When a PySke skeleton manipulates more that one list, all these lists should have the same distribution.

C. Distribution Changing Skeletons

Some skeletons exist to change the distribution of a list. get_partition makes the distribution of the lists visible in the structure itself. For example, get_partition on the list of type PList in Figure 2 yields the PList [[0, 2, 4], [START_REF] Philippe | PySke: Algorithmic skeletons for Python[END_REF][START_REF] Bougé | The data parallel programming model: A semantic perspective[END_REF], [START_REF] Cosmo | A "minimal disruption" skeleton experiment: seamless map & reduce embedding in OCaml[END_REF][START_REF] Chakravarty | Accelerating Haskell array codes with multicore GPUs[END_REF], [START_REF] Emoto | An Automatic Fusion Mechanism for Variable-Length List Skeletons in SkeTo[END_REF][START_REF] Légaux | Development Effort and Performance Trade-off in High-Level Parallel Programming[END_REF]] (global view). Now each processor contains only one element (local size 1), but this element is a list. The distribution of this new list is

[1, 1, 1, 1].
flatten is the inverse operation of get_partition. Communications are needed to compute the new values of all the attributes but "content".

The filter skeleton that filters out elements in a list that does not satisfy the given predicate can easily be implemented using these skeletons: def filter(self, predicate): return self.get_partition() .map(lambda l: l.filter(predicate)) .flatten() Note that filter is not recursive: this code is the method of PList, while in the application of map, filter is the sequential filter of class SList. Also the distribution of the list returned by filter may be unbalanced: all the elements of one processor may have been filtered out while the other processors keep their contents. One very important feature of filter is that it is implemented only using other PySke skeletons. We call such skeletons, derived skeletons. They are skeletons that PySke users could define. Most other skeleton libraries require skeletons such as filter to be defined at the implementer's level. Having a set of skeletons expressive enough to be able to define derived skeletons makes PySke more maintainable than other libraries. It also allows to have intermediate level users who do not need to know the implementation details of PySke but can still provide derived skeletons to less skilled users.

An example of use follows:

>>> par_list = PList.init(lambda x: x, 10) >>> par_list.filter(is_even).to_seq() [0, 2, 4, 6, 8]

The balance skeleton simply makes a list evenly distributed. This skeleton perform communications to exchange content between processors.

III. NEW LIST SKELETONS

This section present new list skeletons for PySke.

A. The distribute Skeleton

The distribute skeleton is a generalization of the balance skeleton. The balance skeleton assumes that the target distribution is the default distribution when parallel lists are created: each processor contains at most one more value than each other processor. However, such a distribution may not be appropriate.

First, the size of the elements in the parallel list may be different (for example if we have a list of matrices of different sizes) and/or the computational resources needed to deal with the elements in the list depends on the values or on the size of each element. In this case a non evenly distribution of data may lead to more balanced computations.

Second, some skeletons that are applied to both the current parallel list object and an additional parallel list object require that the distributions of the parallel lists are identical. Rather than balancing both parallel list, it may be more efficient to redistribute one of the lists so its distribution matches the distribution of the other list.

The skeleton distribute takes an argument: a sequential list that represents the target distribution. This list should have a length equal to the number of processors, the values contained in this list should be positive (possibly zero), and the sum of these values should be equal to the global size of the parallel list the distribute method is applied to.

This skeleton proceeds as follows: first it computes the global indices of the first and last element for each processor, both in the source distribution and in the target distribution. This gives either an interval of indices, or the empty interval ∅. Then each processor intersects its source global interval with the target global intervals of all processors (including itself). These global intervals are offset to local intervals using the start index of the source parallel list. Finally these intervals are used to get the slices of the local contents of the parallel list to send to other processors. After a call to the MPI

B. New Derived Skeletons

The distribute skeleton allows to implement communication oriented derived skeletons as shown in Figure III-B. Such skeletons are often present in other libraries, but there are implemented at a low level: in SPMD using directly MPI calls. In our case, the implementation of these derived skeletons could have been written by a PySke user without knowing the implementation details of PySke. This has two main advantages:

• the implementation is more maintainable,

• the PySke implementation is also more portable. As we mentioned earlier, the programming model of PySke is not tied to MPI. For example, a shared memory implementation of PySke is possible. With a high-level implementation of scatter and gather, it is not necessary to re-implement them for the shared memory version of PySke. The implementations of Figure III-B would work with any new implementation of the nonderived algorithmic skeletons. Informally, gather gathers all the elements of the current list to the given processor. scatter scatters all the elements hold by the specified processor to all processors. The returned list is evenly distributed. scatter_range is a variant of scatter. Instead of specify a processor identifier as the source of the elements to scatter, these elements are specified by their indices in the list, as a range.

The code of Figure III-B only uses already described skeletons. The distribution information is actually defined as an object of Distribution. This class inherits from 'list', and it provides the static method balanced that creates an evenly distributed distribution for a list of the given size. return value is not None selected = self.mapi(select)

.filter(not_none) size = selected.length() distr = Distribution.balanced(size) return selected.distribute(distr) Fig. 4. Derived Algorithmic Skeletons par.procs() is not a skeleton but is part of a small set (in module par) of utility functions. procs() simply returns the list of available processors.

C. The permute Skeleton

The permute skeleton aims at reordering the elements of a list without changing its distribution. It takes as argument a bijective function on the set of the list indices.

For example:

>>> SList([1, 2, 3]).permute(lambda x: 2-x) [3, 2, 1]
The implementation of permute on PList proceeds as follows: First, each element in the local lists are tupled with their new index in the global list, and the processor that holds this index, Then, all these tuples are then grouped by processor of destination. The obtained lists of pairs (new index, value) are exchanged using MPI alltoall. Finally, the local contents are updated using these lists of pairs.

The steps manipulating lists are implemented using SList primitives.

1 def fft(input_list: PList[float]) -> PList[complex]: 2
size, nprocs = len(input_list), len(par.procs()) 3 log2_s, log2_p = int(math.log2(size)), int(math.log2(nprocs)) 4 result = input_list.map(complex) 5

for j in range(0, log2_p): 6 permutation = result.get_partition()\ 7

.permute(partial(_bit_complement, log2_p -j -1))\ 8

.flatten() 9 result = permutation.map2i(partial(_combine, size, log2_s, j), result) 10 for j in range(log2_p, log2_s): 11 permutation = result.get_partition()\ 12

.map(lambda l: l.permute(partial(_bit_complement, log2_s -j -1)))\ 13.flatten() 14 result = permutation.map2i(partial(_combine, size, log2_s j), result) 15 return result To compute the processor that holds an index, it is first necessary to compute the scanr variant of the prefix sum of the distribution of the parallel list. In this variant, the result list has the same length as the input list, but its first element is the first element of the initial list. scanr takes only one argument in addition to the input list: a binary operation.

For example, if we consider the list of Figure 2, the prefix sum scanr returns the list [START_REF] Pelagatti | Structured Development of Parallel Programs[END_REF][START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF][START_REF] Dalcin | Parallel distributed computing using Python[END_REF][START_REF] Cosmo | Skeletal Parallel Programming with OcamlP3l 2.0[END_REF]. Then for a given index i, we search the value v in this list such that i is smaller than v and i is greater or equal to the value preceding v in the list. For example, for index 6, this value is 7. The index of 7 in the list being 2, the processor holding the index 6 is processor 2.

IV. EXPERIMENTS

We experiment on a PySke parallel implementation of the Fast Fourier Transform (Figure 5). The algorithm needs log n steps to complete, where n is the size of the input list. Both the number of processors and n should be a power of two. Basically, the computation needs to combine successively elements at indexes that correspond to a butterfly network. As the list is distributed on processors, the algorithm has therefore two main phases. In the first phase (lines 5-9), it is necessary to permute whole partitions of the current parallel list. In the second phase (lines [START_REF] Cosmo | A "minimal disruption" skeleton experiment: seamless map & reduce embedding in OCaml[END_REF][START_REF] Loulergue | Implementing Algorithmic Skeletons with Bulk Synchronous Parallel ML[END_REF][START_REF] Chakravarty | Accelerating Haskell array codes with multicore GPUs[END_REF][START_REF] Matsuzaki | A Library of Constructive Skeletons for Sequential Style of Parallel Programming[END_REF][START_REF] Emoto | An Automatic Fusion Mechanism for Variable-Length List Skeletons in SkeTo[END_REF], it is necessary to permute elements only in the same local partition. As PySke provides the same skeletons for both parallel lists and sequential lists, and that in a parallel list the local partitions are sequential lists, both phases are actually very similar. _bit_complement computes the correct indices that correspond to the butterfly network. _combine actually computes the new complex number based on two complex numbers that are linked in the butterfly network, and a complex number on the unit circle that only depends on the size of the input list, and the indexes of both complex numbers.

The experiments have be conducted on a shared memory machine with two Intel Xeon E5-2683 v4 processors each q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 6 is the mean of a series of 30 measures. The size of the input list was 2 18 floating point numbers. The speed-up is a relative speed-up with respect to the timing for 1 processor. PySke achieves a good relative speed-up in this case.

V. RELATED WORK Algorithmic skeletons are, of course, naturally combined with functional languages. In this context, an algorithmic skeleton is simply a higher-order function implemented in parallel. There are several skeleton approaches available for functional programming languages. For OCaml, OCamlP3L [START_REF] Cosmo | Skeletal Parallel Programming with OcamlP3l 2.0[END_REF] and its successor Sklml offer a set of a small set data and task parallel skeletons. Both rely on imperative features of OCaml. parmap [START_REF] Cosmo | A "minimal disruption" skeleton experiment: seamless map & reduce embedding in OCaml[END_REF] is a lightweight skeleton library that provides only parallel map and reduce on shared memory machines. Algorithmic skeletons for OCaml have also been implemented using the Bulk Synchronous Parallel ML library [START_REF] Loulergue | Implementing Algorithmic Skeletons with Bulk Synchronous Parallel ML[END_REF]. The proposed library is similar to PySke, but PySke has a much richer set of skeletons. For Haskell, Accelerate is a skeleton library that targets GPUs only. The initial proposal [START_REF] Chakravarty | Accelerating Haskell array codes with multicore GPUs[END_REF] featured classical data parallel skeletons (map and variants, reduce, scan and permutation skeletons) on multi-dimensional arrays. Both these languages are statically typed languages, more difficult to learn, and less used than Python.

PySke is related to SkeTo [START_REF] Matsuzaki | A Library of Constructive Skeletons for Sequential Style of Parallel Programming[END_REF], [START_REF] Emoto | An Automatic Fusion Mechanism for Variable-Length List Skeletons in SkeTo[END_REF], OSL [START_REF] Légaux | Managing Arbitrary Distributions of Arrays in Orléans Skeleton Library[END_REF], [START_REF] Légaux | Development Effort and Performance Trade-off in High-Level Parallel Programming[END_REF], and Muesli [START_REF] Ciechanowicz | The Münster Skeleton Library Muesli -A Comprenhensive Overview[END_REF], [START_REF] Ciechanowicz | Enhancing Muesli's Data Parallel Skeletons for Multi-core Computer Architectures[END_REF], [START_REF] Wrede | Fish school search with algorithmic skeletons[END_REF]. All these libraries are C++ libraries. They are, of course, more efficient than PySke, but using them is less flexible, and requires a higher level of expertise. PySke contains all the skeletons of OSL but the BSP homomorphism skeleton [START_REF] Légaux | Managing Arbitrary Distributions of Arrays in Orléans Skeleton Library[END_REF] and the exception forwarding skeleton [START_REF] Légaux | OSL: an algorithmic skeleton library with exceptions[END_REF]. However, now that base infrastructure of PySke is solid, adding these two skeletons will not be a complex task. PySke contains additional skeletons on lists with respect to OSL.

Muesli provides a permutePartition skeleton on arrays that is limited to permuting whole partitions. The permute skeleton of PySke is more fine grained but used in combination with get_partition, it can implement the same behavior than Muesli's permutePartition as shown by the FFT example. Muesli also offers gathering and scattering skeletons: they are not derived skeletons as our scatter and gather skeletons are.

Both compared to OSL and Muesli, PySke has the advantage of providing skeletons on trees. This is a feature shared with earlier versions of SkeTo, but not with recent ones. PySke is therefore unique in this respect. SkeTo and Muesli provide skeletons on 2D matrices in addition to 1D arrays.

VI. CONCLUSION AND FUTURE WORK

Programming with algorithmic skeletons is a productive way of writing parallel programs. PySke is a simple library for non-expert programmers, but thanks to its rich set of skeletons, it is very expressive. Experiments shows performances are scalable.

We plan to continue adding new skeletons for the existing PySke data-structures, and new data structures such as multidimensional arrays.

As the set of skeletons offered by PySke grows, PySke users may find it difficult to choose the right skeleton, and to rightly compose them. In order to ease the work of PySke users, we plan to design and implement a mechanism for skeleton composition optimization by program transformation. This mechanism will be based on concepts of term rewriting systems [START_REF]Term Rewriting Systems[END_REF], and will be performed at runtime. Preliminary results on non-automatic transformations showed the potential of such an approach [START_REF] Philippe | PySke: Algorithmic skeletons for Python[END_REF].

Fig. 2 .

 2 Fig. 2. Parallel Lists: API View vs. Implementation View

Fig. 3 .

 3 Fig. 3. Example of Execution of the distribute Skeleton

Fig. 5 .

 5 Fig. 5. Fast Fourier Transform

We present the code and results as given in the Python interpreter, omitting imports. >>> represents the prompt.


J L Ramírez

A Rubiano

N Jouandeau

L Gallimard

O Polit

Requirements for artificial muscles to design robotic fingers

This work is part of the ProMain project that concerns the modeling and the design of a soft robotic hand prosthesis, actuated by artificial muscles and controlled with surface Electromyography (EMG) signals. In a first stage, we designed a robotic finger based on the equivalent mechanical model of the human finger. The model takes into account three phalangeal joints, flexion and extension movements are studied. The robotic finger has three Degrees of Freedom (DoF). The finger is designed to be under-actuated and driven by tendons, i.e. only one servo motor actuates the whole finger, and the motor is coupled to the finger mechanism through two flexible wires. As the aim is to design a robotic hand prosthesis that uses artificial muscles, we propose and carry out two experiments to characterize the specifications of the actuator. The first experiment measures the pinch force of the human finger, and the second measures the achieved force using our robotic finger and five different servo motors. It allows us to enhance experimental results with the mathematical model of the finger, to identify the requirements of the artificial muscle.

Introduction

In the domain of robotic prostheses, the main challenge is the design of well sized mechatronic limbs and smart controllers that should help people to achieve desired movements. Up to now, several robotic hands have been designed to accomplish dexterous manipulation tasks, imitating the behavior of the human hand. One example is the KH Hand [START_REF] Mouri | Developments of new anthropomorphic robot hand and its master slave system[END_REF], which is anthropomorphic, uses DC motors and has five fingers and 15 Degrees of Freedom (DoF). Likewise, the Shadow Hand [START_REF] Röthling | Platform portable anthropomorphic grasping with the bielefeld 20-dof shadow and 9-dof tum hand[END_REF] has five fingers and 20 DoF, and uses pneumatic actuators. Furthermore, the HIRO III hand [START_REF] Endo | Five-fingered haptic interface robot: Hiro iii[END_REF] has 15 DoF and is actuated by DC motors, and has an interesting application because it was designed to be a haptic device with force feedback in the fingertips.

However, despite their potential, the hands present in the state of the art have drawbacks that hinder them to be used as a hand prostheses. The main disadvantage is the rigid behavior. The human hand is versatile and due to its adaptability and flexibility it is able to achieve different tasks in an accurate way. Consequently, to design a robotic hand prosthesis, it is necessary to add adaptability and flexibility as a main design requirement.

The utilization of smart and soft materials has led to the development of new adaptive devices (known as soft robots [START_REF] Nurzaman | Soft robotics [tc spotlight[END_REF]) for physical rehabilitation and improvement of human skills. Particularly in the field of soft robotic hands, some works have been conducted: the UB Hand [START_REF] Ficuciello | Experimental evaluation of the UB Hand IV postural synergies[END_REF] and the Pisa/IIT hand [START_REF] Catalano | Adaptive synergies for the design and control of the Pisa/IIT SoftHand[END_REF], which are good examples of the development achieved. Even the DLR Hand II [START_REF] Reinecke | Experimental comparison of slip detection strategies by tactile sensing with the biotac R on the dlr hand arm system[END_REF] can be considered to has a soft behavior due to its driving mechanism that is based on tendons.

Concerning the smart materials for artificial muscles, significant studies are currently being pulled off to develop new actuation strategies, but these technologies are still far from implementation in anthropomorphic robotic hands [START_REF] Melchiorri | Development of the UB Hand IV: Overview of Design Solutions and Enabling Technologies[END_REF]. The first step in the development of artificial muscles for a robotic hand prosthesis is the determination of the actuator requirements. Consequently, ProMain project intends the development of a soft robotic hand prosthesis, with artificial muscles and a control system based on surface EMG signals.

The artificial muscle's requirements are closely related to the robot's features as the used mechanism, the desired level of softness, and the task to be performed. Therefore, it is necessary to develop a mathematical model to describe the robot and the task, allowing the analysis in a more general way. The mathematical models for robotic hands are based on the kinematic and dynamic analysis of its articulated chains. Usually, the Denavit-Hartemberg (DH) parameters [START_REF] Denavit | A kinematic notation for lower-pair mechanisms based on matrices[END_REF] are used to model any kind of articulated chains. The main drawback of the DH parameters utilization is the limitation of active rotations because they are only possible around zaxes of each reference, and a new frame must be added if an extra rotation is needed. Furthermore, the use of homogeneous matrices can cause singularities due to the presence of non-considered rotations introduced with soft materials.

Therefore, in this paper, we adopt the kinematic model DHKK-SRQ [START_REF] Ramirez | Hybrid kinematic model applied to the under-actuated robotic hand prosthesis promain-i and experimental evaluation[END_REF], that allows to know the complete kinematic of the finger, considering the soft behavior of the robot. The model uses the so-called DHKK convention [START_REF] Hugel | Automatic generation of humanoids geometric model parameters[END_REF] (which corresponds to the DH parameters, modified by Khalil and Kleinfinger) and express the soft rotations using quaternions. Concerning the task to be performed, the interaction with the object is very important, and it is usually modeled using the postural [START_REF] Ficuciello | Experimental evaluation of the UB Hand IV postural synergies[END_REF] and the adaptive synergies [START_REF] Catalano | Adaptive synergies for the design and control of the Pisa/IIT SoftHand[END_REF]. The DHKK-SRQ also takes into account the shape and inertial forces of manipulated objects. The proposed dynamic model uses the principle of virtual displacements and virtual works [START_REF] Ramirez | Morphological optimization of prosthesis' finger for precision grasping of little objects[END_REF]. Artificial muscle's requirements are linked with the robot mechanism, so we use the finger of the ProMain-I [START_REF] Ramirez | Hybrid kinematic model applied to the under-actuated robotic hand prosthesis promain-i and experimental evaluation[END_REF] hand as the prototype. This finger is the result of a morphological optimization process [START_REF] Ramirez | Morphological optimization of prosthesis' finger for precision grasping of little objects[END_REF] and consists of three phalangeal joints with flexion -extension, (i.e. 3 DoF). The finger is designed to be under-actuated and driven by tendons, so, only one servo motor actuates the whole finger, and the motor is coupled to the mechanism through flexible wires.

Two experiments are presented to identify the actuation characteristics, the first measures the human force and the second evaluates if classical servo motors fulfill the human hand requirements. The second experiment consists of two stages as follows:

1. Test the finger with five different servo motors, measuring the fingertip force and calculating the motor torque using the kinematic and dynamic model. 2. Measure the pinch force, using a test platform where two fingers are placed facing each other (in pinch position), to validate that the grasping force corresponds to the double of the fingertip force.

Experimental data are analyzed using the kinematic and dynamic models to define the characteristics of the artificial muscle. The obtained results indicate that new actuation technologies can fulfill the force, speed and displacement requirements of a robotic hand prosthesis. This progress is encouraging and permits to follow new directions in the research of smart materials for artificial muscles in robotic hands.

Furthermore, the layout of this paper is as follows. First in section 2, we introduce our finger prototype and its modeling. Then, section 3 shows the experimental set-up and the measures. Subsequently, section 4 exposes the data analysis and the definition of the artificial muscle requirements.

Under actuated Robotic Finger ProMain-I

Our prototype is a bio-inspired tendon-driven finger [START_REF] Ramirez | Morphological optimization of prosthesis' finger for precision grasping of little objects[END_REF] composed of three joints: the metacarpophalangeal (MP), the proximal interphalangeal (PIP) and the distal interphalangeal (DIP). All the joints have one DoF to perform flexion and extension. The finger is controlled by only one servo motor, and the drive mechanism uses two tendons for transmitting motion, one for the flexion and one for the extension, as shown in Figure 1. Considering that the tendons are fastened to the motor pulley and the fingertip, the clockwise rotation of the actuator produces the flexion, and the counterclockwise rotation produces the extension.

Due to the under actuation, the rotation angle of the PIP and DIP joints are linked with the rotation angle of the MP joint. The relation between the angles is approximated with θ 2 = 0.23θ 1 and θ 3 = 0.72θ 1 , where θ 1 is the MP joint angle, θ 2 is the PIP joint angle and θ 3 is the DIP joint angle. Furthermore, the parameters l 1 , l 2 and l 3 are the lengths of the proximal, medial and distal phalanges, as shown in Figure 1. For the kinematic analysis, the hybrid DHKK-SRQ model [START_REF] Ramirez | Hybrid kinematic model applied to the under-actuated robotic hand prosthesis promain-i and experimental evaluation[END_REF] is used, and the dy-namic model uses the virtual displacements and the virtual works approach [START_REF] Ramirez | Morphological optimization of prosthesis' finger for precision grasping of little objects[END_REF]. The kinematic and the dynamic models are presented in the following subsections.

kinematic model of the ProMain-I finger

The DHKK convention, allows the representation of open-loop and close-loop kinematic chains, and presents a convenient definition of the axes z i , which corresponds to the rotation axis of the i th joint. The angle of rotation around z i is denoted by θ i , and is applied using the transformation matrix that is described in Eq. ( 1). The matrix i-1 T i results of the application of: 1) a rotation α i around x i-1 , 2) a translation a i along of x i-1 , 3) a rotation θ i around z i and 4) a translation d i along of z i [START_REF] Hugel | Automatic generation of humanoids geometric model parameters[END_REF].

The parameters α i , a i , θ i and d i , are known as the DHKK parameters, a graphical representation of the parameters is shown in Figure 2.

i-1 T i =     cos θ i sin θ i 0 a i sin θ i cos α i cos θ i cos α i sin α i sin α i d i sin θ i sin α i cos θ i sin α i cos α i cos α i d i 0 0 0 1     (1)
Consequently, the kinematic of a robot composed of n joints is the matrix 0 T n , which is a composition of the orientation of the end effector 0 R n , and the position vector [ 0 P x n , 0 P y n , 0 P z n ], as shown in the following expression:

0 T n = n ∏ i=1 i-1 T i =     0 R n 0 P x n 0 P y n 0 P z n 0 0 0 1     ( 2
)
𝑎 𝑖 𝑥 𝑖 𝛼 𝑖 𝑑 𝑖 𝜃 𝑖 𝐽𝑜𝑖𝑛𝑡 𝑖-1 𝐽𝑜𝑖𝑛𝑡 𝑖 𝑃 𝑖-1 ′ 𝑃 𝑖 ′ 𝑃 𝑖-1
Figure 2: Graphical representation of DHKK parameters [START_REF] Ramirez | Hybrid kinematic model applied to the under-actuated robotic hand prosthesis promain-i and experimental evaluation[END_REF].

We consider a quaternion as a vector represented in the hypersphere S 3 [START_REF] Yan | Study of diagonalization on skew self-conjugate matrix in quaternion field[END_REF], which is centered at the origin of the frame. We represent a rotation of a vector r, with value γ, around the vector n as NRN, where N = [[ cos(γ/2), n sin(γ/2) ]] is the quaternion that corresponds to the vector n, R = [[ 0, r ]] is the quaternion that correspond to the vector r, and N is the conjugate of the quaternion N.

Likewise, in the 3D space, a group of rotations with angles α, β , θ applied to the vector r around the axes x, y, z, as shown in Figure 3, is expressed as: XYZRXYZ

where The hybrid model DHKK-SRQ unifies the DHKK parameters with multiple sets of quaternions, using an optimal analysis of the available sensor feedback. Thus, each i th joint (for i = 1, . . . , n where n is the number of joints) is considered as an element that has a hybrid (i.e. rigid and soft) behavior. Consequently, the i th joint is modeled, in a first step as a rigid element, with only the rotation θ i around z i . Subsequently, the rotations α i and β i around axis x i and y i are added using the hypersphere S 3 i centered in the joint frame. As a result, the model can apply rotations in all axes avoiding to add any extra reference frame.

X = [[ cos(α/2), x sin(α/2) ]] Y = [[ cos(β /2), y sin(β /2) ]] Z = [[ cos(θ /2), z sin(θ /2) ]]
Initially, the rotations are performed using the homogeneous matrices that are given by Eq. ( 1), and the kinematic of the rigid joints results from Eq. ( 2). Afterwards, we formulate the extra rotations (those that appear from the low stiffness joints) α i and β i using SRQ. Additionally, if SRQ is active, then the rotation θ i is formulated using a quaternion; so that, the set of quaternions for each i th joint is like the proposed in Eq. ( 3). The equivalent kinematic model of the finger is shown in Figure 4. The frame (x 1 , y 1 , z 1 ) and the hypersphere S 3 1 are associated to the joint 1 (MP). The frame (x 2 , y 2 , z 2 ) and the hypersphere S 3 2 are associated to the joint 2 (PIP). And the frame (x 3 , y 3 , z 3 ) and the hypersphere S 3 3 are associated to the joint 3 (DIP). The frame (x f , y f , z f ) corresponds to the fingertip position. Table 1 shows the DHKK parameters of the ProMain-I finger.

i i i α α α a a a d d d θ θ θ 1 0 0 0 θ 1 2 0 l 1 0 θ 2 3 0 l 2 0 θ 3 f 0 l 3 0 0
Table 1: DHKK parameters of the ProMain-I Finger.

Dynamic model of the ProMain-I finger

The equivalent dynamic model of the finger is shown in Figure 5. w 1 , w 2 and w 3 are respectively the weights of the proximal, medial and distal phalanges, and are placed at the coordinates (x 1 , y 1 ), (x 2 , y 2 ) and (x 3 , y 3 ). f R is the applied force that is equivalent to the reaction force. The proposed dynamic model uses the principle of the virtual displacements and virtual works [START_REF] Ramirez | Morphological optimization of prosthesis' finger for precision grasping of little objects[END_REF]. The virtual work δW is calculated for the external forces (e.g. weight, applied force and input torque) in Eq. ( 4) and the inertial forces (e.g. centrifugal forces) in Eq. [START_REF] Ficuciello | Experimental evaluation of the UB Hand IV postural synergies[END_REF].

δW e = Q T e δ r e (4)
where Q T e is the external forces vector and δ r e is the virtual displacement vector of the forces application points.

δW ı = M qT δ r ı ( 5
)
where M is the diagonal mass matrix composed with the masses m i and inertias J i .

The index i is used to denote the i th joint, and the index ı (i.e. without dot) is used to notice the independent coordinates. qT is the second derivate with respect to the time of the q vector shown in Eq. ( 6), and represents the acceleration vector. And δ r ı is the virtual displacement vector of the inertial frameworks.

q = [x 1 , y 1 , θ 1 , x 2 , y 2 , θ 2 , x 2 , y 2 , θ 1 ] (6)
The dynamic equilibrium is given by δ q T [M q-Q e ] = 0 when substituting Eq. ( 4) and Eq. ( 5) in δW e = δW ı . In order to solve the equilibrium equation, considering the movements restrictions, it is necessary to separate dependent and independent coordinates. The separation is performed using the transformation proposed in Eq. ( 7), as result we have the equilibrium defined by Eq. ( 8), which is separable.

δ q = Bδ q iı , B = -C -1 qd C qı I (7)
where C qd is the jacobian of dependent coordinates, C qı is the Jacobian of independent coordinates and I is the identity matrix.

δ q T iı B T [M q -Q e ] = 0 (8)
Solving Eq. ( 8), we obtain the dynamic function which gives the input torque τ 1 as function of the force f R and the kinematic q, q, q. The resulting expression is shown in Eq. ( 9), where we use the abbreviations C i := cos(θ i ) and S i := sin(θ i ).

τ 1 ( f R , q, q, q) = H 1 -4l 1 H 10 + H 11 (l 2 S 2 -4) 8l 2 S 2 -32 (9)
where

H 1 = 2l 2 1 θ1 (l 2 (m 2 + m 1 )S 2 -4m 1 -6m 2 )S 2 1 H 2 = -m 2 l 2 2 S 2 2 θ2 /2 H 3 = (( ẍ1 + 3 θ2 + 2g)m 2 + (2 θ2 + g)m 3 + m 1 g -2 f R )l 2 S 2 H 4 = (m 2 l 2 C 2 l 1 C 1 θ1 + l 2 2 θ2 (m 2 + 2m 3 )C 2 2 )/4 H 5 = l 2 (2m 3 ÿ2 + m 3 l 3 C 3 θ3 + m 2 ÿ1 )C 2 /2 + 2m 3 l 3 S 3 θ3 H 6 = (-6g -6 ẍ1 )m 2 + (-4 ẍ2 -4g)m 3 + 8 f R -4m 1 g + 2J 2 θ2 H 7 = (m 1 + 3m 2 /2)l 2 1 θ1 C 2 1 /4 H 8 = θ2 l 2 (m 2 + 2m 3 /3)C 2 + 2m 3 l 3 C 3 θ3 /3 H 9 = 4m 3 ÿ2 /3 + 2m 2 ÿ1 H 10 = (H 2 + H 3 + H 4 + H 5 + H 6 )S 1 H 11 = 8(H 7 + 3l 1 (H 8 + H 9 )C 1 /8 + J 1 θ1 )

Experimental set-up

With the aim of defining the actuator characteristics, we held two experiments. The first measures the human pinch force, that is compared with the force achieved by a servo motor through a second test. In the following subsections, both experiments are described.

First experiment: Measure of the human hand pinch force

The pinch force can be considered as the force applied by two fingers of the hand, usually the index and the thumb fingers. The applied force must be adapted to the object's weight, acceleration, surface texture, contour and structure [START_REF] Ramsay | A functional data analysis of the pinch force of human fingers[END_REF]. Consequently, the measure of the pinch force have to be customized to each problem [START_REF] Schreuders | Measurement error in grip and pinch force measurements in patients with hand injuries[END_REF], and that is why in this study we carried out an experiment suited to our requirements.

The force is measured in a group of five healthy subjects between 24 and 32 years old, all of them are males. The sensor is a hand dynamometer Vernier TM D-BTA, suitable to measure the pinch force, whose characteristics are: 1) accuracy of 0.6N, 2) resolution of 0.2141N and 3) operational range from 0 to 600N. The data is collected using a digital oscilloscope connected to a computer. Figure 6 shows the scheme of the experiment. Each subject is asked to apply the maximal pinch force during a period of 15 seconds. To avoid fatigue each subject perform only three trials. Results presented in Table 2 correspond to the mean value of the pinch force performed by each subject, and the corresponding standard deviation of each measure.

Single signal capture Dynamometric Measures

Mean signal analysis

Second experiment: Measure of the robotic finger pinch force

With the objective of establishing if a classic servo motor can fulfill the requirements of the human hand, we design two test platforms. The first one is a single-finger platform and allows to measure the kinematic and the fingertip force. The second one is a two-finger platform and consists of two fingers placed against them. It will be used to validate that the pinch force is equal to the double of the finger force. The CAD models of the test platforms are shown in Figure 7 and Figure 8.

The experiments are performed using five different servo motors, two standard and three serial, listed in Table 3. To measure the force, we used a resistive-based force sensor Flexiforce , that measures up to 5N. The sensor is calibrated in the range of 0.6N to 4.8N and is placed on a support (platform) that is located in the trajectory of the fingertip. Furthermore, to ensure the repeatability, each servo motor is tested five times performing the same movement. Figure 9 shows the fingertip force f R during three repetitions of flexion and extension using the HS-422 servo motor. The force applied by the MX-106R servo motor is higher than the sensor limit and also higher than mean human finger force, so its value is not taken into account. Table 4 shows the measured forces. Table 4: Mean value of the fingertip force (one-finger platform).

The servo motor that performs the force closest to the human finger force is chosen as actuator for the two-finger platform. Consequently, we took the HS-422 servo motor for the test with the two-finger platform. Consequently, the tests with the second platform are conducted for four different distances between fingers (50, 55, 60 and 65mm). For each distance the test is carried out five times. Table 5 shows the measured grip force using the two-finger platform. Considering that the finger performs flexion and extension in 2D, the kinematic is measured using a high-performance CCD camera Prosilica GE-2040, which tracks black markers placed on the finger joints and the fingertip. Figure 10 shows the scheme of the experiment. The kinematic data is correlated with the DHKK-SRQ and the dynamic model to identify the torque and displacement requirements of artificial muscles. Kinematic results are presented and analyzed in the next section.

Distance [mm] [mm] [mm] Mean pinch Force [N] [N] [N] Standard deviation [N] [N] [N] 50

Interface for controlling the platform

Requirements and characterization of the artificial muscle

The main features of an actuator are the force, strain and time response. Thus, we need to define these characteristics for the artificial muscle. Considering that our goal is to design a robotic hand that will be able to perform precision grasping movement, the actuator features can be established from measures of the human hand pinch. However, it is important to take into account that the robotic finger mechanism can modify the actuator requirements. Consequently, the proposed approach, to identify the artificial muscle requirements, is defined by three kind of measures: 1) the human pinch force, 2) the human response time and 3) the kinematic and dynamic behavior of the robotic finger ProMain-I.

Pinch force requirements: Considering that our goal is to define a reference value of the human pinch force, we have collected multiple samples from each subject and computed a mean force value for each one, as shown in Table 2, section 3. Thus, the pinch force is in the interval [4.80N, 6.74N], and the mean value is 5.94N. A graphical representation of the analysis is shown in Figure 11. Furthermore, we measured the time response of the subjects, and we obtained a mean value of 0.244s with a standard deviation of 0.06s. Kinematic tracking and resulting force of the ProMain-I finger: As result of the test, we get the position vectors of the joints, i.e. the vector [ 0 P x 1 , 0 P y 1 , 0] for the joint 1 , [ 0 P x 2 , 0 P y 2 , 0] for the joint 2 and [ 0 P x 3 , 0 P y 3 , 0] for the joint 3 . Likewise the vector [ 0 P x f , 0 P y f , 0] corresponds to the fingertip position. Considering that the movement is performed in the plan, 0 P z i is always zero. Figure 12 shows the results of the measured kinematic using the HS-422 servo motor. The standard deviation of the measures of the MP joint (whose position is always zero) is 0.4733mm. Likewise, the standard deviation of the measures of the joints PIP and DIP are 0.1848mm and 0.5598mm respectively. The standard deviation of the measures of the fingertip is 1.6069mm.

Furthermore, when the test is carried out with the two-fingers platform, one finger behaves as an obstacle to the other one. Thus, the experiment shows that the relations of the angles are different from the proposed ones given in section 2. This difference in the angles behavior is not an issue as we are tracking all the movements and updating the kinematic with the DHKK-SRQ model.

Concerning the fingertip force and the pinch force, the lower difference takes place when the distance between the fingers is 60mm. In that case, the force applied by each finger is the half of the total pinch force, that is 4.70N/2 = 2.35N. Thus, the difference with respect to the force measured using the same actuator (with a mean value of 2.69N) and the one-finger platform is 0.34N.

Analysis of results: Summarizing the obtained results, we can state that if we attempt to reproduce the pinch properties of the human hand, the actuator must fulfill the following features:

• The torque should be enough to produce a fingertip force of 2.96N with a time response of 0.244s. • The actuator should be suitable to produce a rotation of π/2 rad.

To calculate the required torque to apply a force of 2.96N using the proposed dynamic model, it is necessary to know the vector q, which corresponds to the dependent and independent coordinates of the robot and it is composed as shown in Eq. ( 6). Taking into account that the orientation of (x i , y i ) is the same of (x i , y i ), we use the relations x i = x i + l i cos θ i /2 and y i = y i + l i sin θ i /2 to calculate the vector q. The values x i , y i and θ i correspond to the measured kinematic information. As a result, applying the dynamic model, and using the kinematic DHKK-SRQ model, we found out that the torque τ 1 ( f R , q, q, q) = 74.5Nmm.

Conclusions

We have presented the requirements for artificial muscles to design robotic fingers. This finger is under-actuated and based on human morphology according to artificial muscles. This finger is a mechanical chain with 3 joints driven with flexible wires to produce more adaptive grasping movements. The proposed kinematic and dynamic models allow us to calculate the torque requirement to mimic human finger force. Experiments compare pinch and achieved forces for 5 different servo motors. According to the obtained results and the performed analysis, we identify that smart materials (e.g. ionic polymer metal composite IPMC or shape memory alloy SMA) can fulfill the requirements of artificial muscles to design robotic hand prostheses. Consequently, we envisage the implementation of artificial muscles, based on smart materials, in the next prototype of our robotic finger.

Figure 1 :

 1 Figure 1: Finger mechanism.

Figure 3 :

 3 Figure 3: Combination of rotations α, β and θ of π/2 rad around x, y and z. [10].

Figure 4 :

 4 Figure 4: kinematic model of the robotic finger.

Figure 5 :

 5 Figure 5: Dynamic model of the robotic finger.

Figure 6 :

 6 Figure 6: Experimental set-up to measure the human pinch force.

Figure 7 :Figure 8 :Table 3 :

 783 Figure 7: CAD Model of the single-finger platform.

Figure 9 :

 9 Figure 9: Results of the fingertip force using the HS-422 servo motor.

Figure 10 :

 10 Figure 10: Experimental set-up to track the kinematic and measure applied force of the robotic finger.

Figure 11 :

 11 Figure 11: Human Pinch Force.

Figure 12 :

 12 Figure 12: Results of the position tracking using the HS-422 servo motor.

Table 2 :

 2 The mean value of the human pinch force.

Table 5 :

 5 Mean pinch force (two-finger platform).

		4.02	0.02
	55	4.62	0.08
	60	4.70	0.05
	65	3.54	0.06

Acknowledgements Through this acknowledgment, we express our sincere gratitude to the Université Paris Lumières UPL for the financial support through the project PROMAIN. This work has been partly supported by Université Paris Lumières UPL and by a Short Term Scientific Mission funding from LEME-UPO-EA4416 / LIASD-UP8-EA4383. We also acknowledge Colciencias -Colombia and the Universidad Militar Nueva Granada for the financial support of the PhD students.