more recently discrete differential geometry has also been proposed as a tool for developing consistent exact numerical models of rods [START_REF] Bergou | Discrete elastic rod[END_REF].

Beam theories can be obtained employing a reduction process from the 3D continuum, or introducing specific assumptions for the kinematics, and for the strain energy, the latter models usually are referred to as structural models. In the paper, we will follow the latter strategy. The definition of the structural model implies therefore the specification of its configuration space, the set of the possible geometrical configurations V of the model endowed with a smooth chart (generally, an atlas). Smoothness indicates that the change of coordinates between any two charts defined on the set V is differentiable. The dimension of the configuration space characterizes the structural model at hand. The space of the admissible variations of the configuration, necessary for setting up the principle of virtual work, is given by the tangent bundle to V. This aspect is fundamental in the nonlinear analysis of space rods, since the tangent operator is needed in the implementation of the iterative method of solution of the nonlinear equilibrium equations.

In a geometrically exact beam model, the configuration of the space rod is mapped by the position of the centerline and by the finite rotations of three directors attached to the beam axis. In this way, spatial rod models represent a one-dimensional particularization of the general polar continuum, also known as Cosserat continuum. In the case when the rotations of the directors are independent, we get the Timoshenko beam or its generalization to the spatial case ( [START_REF] Simo | A finite beam formulation, The three dimensional dynamic I[END_REF][START_REF] Simo | A three-dimensional finite strain rod model, Part II: Computational aspects[END_REF][START_REF] Smolenski | Statically and kinematically exact nonlinear theory of rods and its numerical verification[END_REF][START_REF] Simo | On the dynamics in space of rods undergoing large motions. A geometrically exact approach[END_REF]). Analogously, a three-dimensional nonlinear polar model with drilling rotations was developed for shells in [START_REF] Eremeyev | Phase transitions in thermoelastic and thermoviscoelastic shells[END_REF][START_REF] Eremeyev | The nonlinear theory of elastic shells with phase transitions[END_REF][START_REF] Eremeyev | Thermomechanics of shells undergoing phase transition[END_REF][START_REF] Fox | A drill rotation formulation for geometrically exact shells[END_REF][START_REF] Ibrahimbegovic | Stress resultant geometrically nonlinear shell theory with drilling rotations-Part I. A consistent formulation[END_REF][START_REF] Ibrahimbegovic | Stress resultant geometrically nonlinear shell theory with drilling rotations -Part II. Computational aspects[END_REF][START_REF] Pietraszkiewicz | Extended non-linear relations of elastic shells undergoing phase transitions[END_REF], and without drilling rotation in [START_REF] Altenbach | On the effective stiffness of plates made of hyperelastic materials with initial stress[END_REF][START_REF] Bersani | Buckling of an elastic hemispherical shell with an obstacle[END_REF][START_REF] Simo | On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization[END_REF][START_REF] Simo | On a stress resultant geometrically exact shell model. Part II: the linear theory; Computational aspects[END_REF]. The space SO(3) of finite rotational degrees of freedom is, however, noncommutative, so that their increments cannot be simply added, for details see [START_REF] Argyris | An excursion into large rotations[END_REF][START_REF] Atluri | Rotations in computational solid mechanics[END_REF][START_REF] Rodrigues | Des lois géométriques qui régissent les déplacements d'un systéme solide dans l'espace et de la variation des coornnées provenant de ces déplacements considéres indépendment des causes qui peuvent les produire[END_REF][START_REF] Simo | Formulations of finite elasticity with independent rotations[END_REF]. This problem has been strongly addressed in computational mechanics by Crisfield and Jelenic [START_REF] Crisfield | Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation[END_REF] and by Jelenic and Crisfield [START_REF] Jelenic | Geometrically exact 3D beam theory: Implementation of a strain-invariant finite element for static and dynamics[END_REF] who, in order to overcome the problem, proposed to interpolate the local rotations with respect to a corotational frame. They also pointed out that the strain measures must remain unaffected by a superimposed rigid rotation of the configuration (objectivity) also in the interpolation scheme. Later alternative methods to the corotational scheme have been proposed (see [START_REF] Eugster | Director-based beam finite elements relying on the geometrically exact beam theory formulated in skew coordinates[END_REF][START_REF] Ibrahimbegovic | On the geometrically exact formulation of structural mechanics and its applications to dynamics, control and optimization[END_REF]). A representation based on the natural (Bishop) frame for shear undeformable space rods under Kirchhoff-Love hypotheses has been introduced by Langer and Singer [START_REF] Langer | Lagrangian aspect of the Kirchhoff elastic rod[END_REF] and subsequently used by Bergou et al. as starting point for an exact rod model built up on discrete differential geometry [START_REF] Bergou | Discrete elastic rod[END_REF].

As a consequence of the fact that the configuration space for the rod model is a non-Euclidean manifold, the structural coefficients of the space (i.e., the Christoffel symbols of the connection) do not vanish; therefore, the directional Gateaux derivative loses the covariant character. The Gateaux derivative is usually employed in continuum mechanics for obtaining the variation of the configuration, while the second Gateaux derivative of the total potential energy yields the tangent stiffness operator needed in incremental analyses. Since the latter does not represent the second covariant gradient, the tangent operator calculated in this way loses its symmetry properties.

The issue was raised by Simo [START_REF] Simo | The (symmetric) Hessian for geometrically nonlinear models in solids mechanics: Intrinsic definition and geometric interpretation[END_REF][START_REF] Simo | The Hamiltonian structure of nonlinear elasticity: The material and convective representations of solid, rods, and plates[END_REF], who first noted that the tangent operator evaluated using the Gateaux derivative is symmetric only at an equilibrium configuration, if the load is conservative. He also described a procedure for evaluating the Hessian of the total potential energy that does not require the explicit knowledge of the symbols of the connection that amounts to a straightforward symmetrization of the operator obtained performing the directional derivative. These facts are indicated in detail in [START_REF] Crisfield | Non-linear Finite Element Analysis of Solids and Structures[END_REF]. Later, Jelenic [START_REF] Jelenic | A kinematically exact space finite strain beam model-finite element formulation by generalized virtual work principle[END_REF] and others authors afterward used this conclusion. The evaluation of the symmetric tangent operator for an exact structural model is a necessary requirement in order to extend the formulation to beams composed by materials with local symmetries, as occurs in functionally graded beams (see for instance [START_REF] Eremeyev | Local symmetry group in the general theory of elastic shell[END_REF] for the use of material local symmetry groups in structural models, and [START_REF] Birsan | Deformation analysis of functionally graded beams by the direct approach[END_REF] for the theory of functionally graded beams).

In a general sense, in order to obtain the consistent tangent operator of the energy functional, it is needed to define a connection on the tangent bundle to the configuration manifold and then to calculate the second covariant gradient. It is symmetric only if the connection is torsion-free; therefore, the unique torsion-free connection consistent with the adopted metric has to be chosen (i.e., the Levi-Civita connection). This fact, consisting in the fundamental theorem of the Riemman's geometry, already present in [START_REF] Crisfield | Non-linear Finite Element Analysis of Solids and Structures[END_REF][START_REF] Simo | The (symmetric) Hessian for geometrically nonlinear models in solids mechanics: Intrinsic definition and geometric interpretation[END_REF][START_REF] Simo | The Hamiltonian structure of nonlinear elasticity: The material and convective representations of solid, rods, and plates[END_REF], allows to overcome the noncommutative nature of the rotation group.

The cited works refer to shear deformable rods, whose numerical implementation in a finite element (FE) scheme presents, as it is well-known, shear locking. Kirchhoff -Love rod models are able to avoid shear locking, so they can be effectively used in modeling bending and torsion-dominated nonlinear problem. However, membrane and torsional locking may be present, as it has been illustrated in [START_REF] Cuomo | Isogeometric analysis of space rods: Considerations on stress locking[END_REF] where a family of mixed finite elements for curved elements, based on isogeometric interpolation, has been introduced that avoids both membrane locking and the effect of flexure locking on the twisting energy terms. Recently, it has been shown in [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF] that using a hierarchic shell approach, it is possible to develop a numerical model free of shear locking, (obtaining a sort of switch element) affected only by membrane locking. On this line, the unshearable model used in this work can be considered an interesting base for developing efficient numerical models for rods. Kirchhoff -Love model offers several other advantages, since there is no need of assuming any constitutive behavior for the shear that would introduce significant errors in the theory. Among the many new interesting problems of mechanical engineering that can be effectively modeled in this way, we recall the simulation of living tissues, like cilia, of wind-induced loads on eolic turbines, the analysis of micro piezoelectric actuators, etc. [START_REF] Andreaus | Piezoelectric passive distributed controllers for beam flexural vibrations[END_REF][START_REF] Alessandroni | A revival of electric analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF][START_REF] Cazzani | On the dynamics of a beam partially supported by an elastic foundation: an exact solution-set[END_REF][START_REF] Dell'isola | Continuum modeling of piezoelectromechanical truss beam: An application to vibration damping[END_REF][START_REF] Dell'isola | Piezo-electromechanical (PEM) structures: Passive vibration control using distributed piezoelectric transducers[END_REF][START_REF] Del Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF][START_REF] Giorgio | Multimode vibration control using several piezoelectric transducers shunted with a multiterminal network[END_REF][START_REF] Maurini | On a model of layered piezoelectric beams including transverse stress effect[END_REF][START_REF] Porfiri | Modeling and design of passive electric networks interconnecting piezoelectric transducers for distributed vibration control[END_REF][START_REF] Porfiri | Circuit analog of beam and its application to multimodal vibration damping, using piezoelectric transducers[END_REF][START_REF] Vidoli | Modal coupling in one-dimensional electromechanical structured continua[END_REF].

The configuration manifold of a Kirchhoff -Love model is different from the one used in the Timoshenko model, since the bending rotations are not independent and are obtained from the position field. Therefore, the configuration space that in the former case is R 3 × SO(3) becomes R 3 × S 1 , so that the derivation of the Hessian outlined in Simo's paper cannot be directly applied and needs some adjustments.

In this work, we analyze in detail and obtain the tangent stiffness operator for the exact spatial Kirchhoff-Love rod model presented in [START_REF] Greco | B-spline interpolation for Kirchhoff-Love space rod[END_REF][START_REF] Greco | An implicit G 1 multi patch B-spline interpolation for Kirchhoff love space rod[END_REF] that is characterized by intrinsically enforcing the Kirchhoff-Love constraints through a split of the rotation operator that defines the unit triad attached to the centroid curve. Differently from standard approaches that adopts the Frenet triad, we consider a more general intrinsic triad that in general is not torsion-free. In particular, we introduce a first rotation operator, indicated with ∈ S 2 , defined as a rotation from the reference to the current unit tangents (i.e., this first unit operator is a finite rotation without drilling around the reference centroid axis). By means of the operator, it is defined a service unit triad that successively can be rotated, through the operator R ∈ S 1 , around the current unit tangent vector in order to define the current unit triad. The rotation operator is fully characterized by the displacement field. In this way, the strain energy of the rod, in particular the bending energy term, depends on the second covariant derivative of the centroid displacement, differently from what happens in shear deformable beams. The required generalization is obtained on the basis of arguments generalizing or paralleling those presented for instance in [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF][START_REF] Dell'isola | The relationship between edge contact forces, double forces and interstitial working allowed by principle of virtual power[END_REF][START_REF] Dell'isola | How conctact interactions may depend on the shape of Cauchy cuts in Nth gradient continua: Approach "à la D'Alambert[END_REF][START_REF] Ferretti | Modeling the onset of shear boundary layers in fibrous composite reinforcements by second gradient theory[END_REF]. The tangent operator obtained in this way is fundamental for the iterative incremental analysis; further, higher expansion of the potential is necessary for the study of the post-critical behavior, by means of perturbation methods, as shown in the many contributes given by Luongo and coworkers, [START_REF] Luongo | Multiple interaction and localization phenomenon in postbuckling of compressed thin-walled members[END_REF][START_REF] Luongo | On the perturbation analysis of interactive buckling in nearly symmetric structures[END_REF][START_REF] Luongo | Perturbation methods for nonlinear autonomous discrete-time dynamic systems[END_REF][START_REF] Luongo | Multiple scale bifurcation analysis for finite-dimensional autonomous systems[END_REF][START_REF] Luongo | Mode localization in dynamics and buckling of linear imperfect continuous structures[END_REF]. The method described in the paper can be usefully employed for these higher expansions.

Kinematics

In this section, a brief summary of the finite kinematics of a Kirchhoff-Love rod and of its tangent form is given, which is described in detail in [START_REF] Greco | B-spline interpolation for Kirchhoff-Love space rod[END_REF][START_REF] Greco | An implicit G 1 multi patch B-spline interpolation for Kirchhoff love space rod[END_REF]. Following this approach, the rod is defined by its centroid curve and by the directors lying in the plane orthogonal to the centroid line.

Reference configuration

The centroid line of the rod is mapped from an interval A by means of the position vector p 0 : A ⊂ R → R 3 . The tangent to the reference position is

t0 := d p 0 dS , t0 = 1, ( 1
)
where S is the arc length of the reference centroid curve. Let n0 : A → R 3 be a director orthogonal to t0 , and ν0 := t0 × n0 , so that the intrinsic unit triad of the rod in the reference configuration is given by { t0 , n0 , ν0 }.

The normal vector field n0 (S) is chosen arbitrarily. In particular, the natural frame, defined as a torsion-free frame, can be selected, or any other convenient choice is possible. In the paper, we refer to a generic reference frame, with an initial twist. The curvatures of the centerline are defined on the reference configuration as

χ n0 = - d t0 dS • ν0 , χ ν0 = d t0 dS • n0 , χ t0 = d n0 dS • ν0 (2)
The position of the generic point of the rod is given by *

p 0 = p 0 + ϑ n n0 + ϑ ν ν0 , ( 3
)
where ϑ n and ϑ ν are the quota functions, independent of the arc length.

Current configuration

The position vector for the points of the centerline is p : A → R 3 . Denoting by u : A → R 3 the displacement field of the centroid curve, the position vector is given by

p := p 0 + u (4)
and the tangent vector is

t := d p dS , t := t t , t := √ t • t. ( 5
)
The orientation of the current triad is obtained applying two subsequent rotations. The first rotation ( t0 , t)(S) brings the unit tangent vector { t0 (S) on t(S)}; that is, it is a rotation without drilling around the tangent to the reference centroid line. The axial vector of this rotation is t0 × t, and Rodrigues formula gives the representation:

( t0 , t) := t0 • t I + t0 × t × I + t0 × t ⊗ t0 × t 1 + t0 • t . ( 6
)
The rotation operator, , defines an intermediate unit triad on the orthogonal plane to the current centroid curve as

{ t, n , ν } := { t0 , n0 , ν0 }, (7)
An additional rotation operator, R( t, φ)(S), is then applied, given by a rotation φ(S) of the intermediate triad around the current tangent to the rod,

R( t, φ) := I + sin[φ] t × I + (1 -cos[φ]) t × t × I . ( 8
)
Notice that in this way, the Kirchhoff-Love constraints are automatically enforced. The spatial triad for the Kirchhoff rod is then

{ t, n, ν} := { t, R n , R ν } = { t0 , R n0 , R ν0 }. (9)
The curvatures of the spatial configuration are defined in a similar way to Eq. ( 2):

χ t = d n dS • ν, χ n = - dt dS • ν, χ ν = dt dS • n b = t × dt dS = χ n n + χ ν ν. ( 10
)
The spatial position of a generic point of the rod is obtained as

* p = p + ϑ n n + ϑ ν ν. (11)
The configuration space of the Kirchhoff rod is then given by

M := {ϕ = ( p, , R) : [0, L] → (R 3 × S 2 ) × S 1 }. ( 12
)
Notice that the dimension of the manifold M is 4, since the position p and the rotation can both be expressed as functions of the displacement field u.

Tangent kinematics

In this section, we derive the variation field for the configuration of the rod. If p is a 1-parameter admissible family of positions of the centerline, such that p =0 = p, the variation of the centerline position ũ : M → T p M is defined as

ũ := d p d =0 . ( 13
)
The variations of the rotations are defined in a similar way. Let Q = R be a curve of admissible rotations for the rod, with = ( t0 , t ) and R = R(φ , t ), then the variations are given by

Q := d Q d =0 = R + R ˜ , R := dR d =0 , ˜ := d d =0 . ( 14
)
whose spatial representation is given by (see [START_REF] Greco | An implicit G 1 multi patch B-spline interpolation for Kirchhoff love space rod[END_REF] for details)

Q Q T = R ˜ T R T + R R T , = t ⊗ t -t ⊗ t + φ t × I, ( 15
)
with

φ := dφ d =0 , t := d t d =0 . ( 16
)
The axial vector of the rotation ( 15) is the spin:

ω = φ t + ωb , ωb = ωn n + ων ν. (17)
The bending rotation ωb , accounting for the Kirchhoff constraints, is given by the following:

ωn = - 1 t d ũ dS • ν, ων = 1 t d ũ dS • n. (18)
The tangent space to a configuration ϕ is then

T ϕ M := { φ := ( ũ, ωb , φ) : [0, L] → (R 3 × SO(2)) × SO(1)}. ( 19
)

Equilibrium operator

Assuming a linear elastic constitutive behavior (as it is done in [START_REF] Greco | B-spline interpolation for Kirchhoff-Love space rod[END_REF][START_REF] Simo | A three-dimensional finite strain rod model, Part II: Computational aspects[END_REF]) and approximating the transport operator for the stresses along the normal directors with the identity, the elastic potential energy of the rod takes the form

ψ int (ϕ) = 1 2 L 0 E A ε 2 (ϕ) + E I n χ 2 n (ϕ) + E I ν χ 2 ν (ϕ) + G J t χ 2 t (ϕ) dS, ( 20
)
where the finite strain measures are respectively

ε(ϕ) := t • t -1 2 , χ t (ϕ) := d n dS • ν - d n0 dS • ν0 χ n (ϕ) := dν dS • t - dν 0 dS • t0 , χ ν (ϕ) := - d n dS • t - d n0 dS • t0 . ( 21
)
A conservative potential ψ ext is assumed for the external loads.

Weak formulation

The first Gateaux derivative of the total potential energy ψ int + ψ ext gives the weak formulation of the equilibrium equations. In the following, in order to simplify the notation, the symbol indicating the directional Gateaux derivative is substituted with

∂(•) = d(•) d =0 (22)
The internal virtual power is then

∂ ϕ ψ int (ϕ) = L 0 (N (ϕ)ε(ϕ) + M n (ϕ) χn (ϕ) + M ν (ϕ) χν (ϕ) + M t (ϕ) χt (ϕ)) dS ( 23
)
where the constitutive relations and the variations of the strain measures are

N (ϕ) = E Aε(ϕ), ε(ϕ) := dε(ϕ ) d =0 = d ũ dS • t M n (ϕ) = E I n χ n (ϕ), χn (ϕ) := dχ n (ϕ ) d =0 = - d 2 ũ dS 2 + φ b • ν M ν (ϕ) = E I ν χ ν (ϕ), χν (ϕ) := dχ ν (ϕ ) d =0 = d 2 ũ dS 2 -φ b • n M t (ϕ) = G J t χ t (ϕ), χt (ϕ) := dχ t (ϕ ) d =0 = d φ dS + 1 t 2 b • d ũ dS ( 24
)
in which the second covariant line gradient is defined as

d 2 ũ dS 2 = t d dS 1 t d ũ dS , = d dS d ũ dS - 1 t d t dS d ũ dS . ( 25
)
The variation of the strain measures can be cast in the equivalent form:

⎛ ⎝ ε t × χ χt ⎞ ⎠ = ⎛ ⎜ ⎝ t • d dS 0 -(I -t ⊗ t) d 2 dS 2 b b t 2 • d dS d dS ⎞ ⎟ ⎠ ũ φ (26)
Notice that the Gateaux derivative that we have performed coincides with the covariant derivative, since ψ int is a scalar.

Strong formulation

Applying the Green-Gauss-Ostrogradsky's formula to the internal power, we obtain the first directional derivative of the strain energy potential ψ along the generic tangent direction φ ∈ T ϕ M, whose stationarity conditions yield the equilibrium equations of the rod:

∂ ũψ, ∂ φ ψ ũ φ = {0, 0}, ∀ ũ, ∀ φ. ( 27
)
The strong form of the equations is as follows:

∂ ũ(ψ ) = L 0 - d dS N t + 1 t d dS t t × M + M t b t -f dS ∂ φ (ψ) = L 0 - d M t dS + t × M • b -m t dS. ( 28
)
with the boundary conditions

N eff (L) = f t (L), or ut (L) = 0, -N eff (0) = f t (0), or ut (0) = 0, T n (L) = f n (L), or un (L) = 0, ( 29
) -T n (0) = f n (0), or un (0) = 0, T ν (L) = f ν (L), or uν (L) = 0, -T ν (0) = f ν (0), or uν (0) = 0, where f t = t • f , f n = n • f , and f ν = ν • f ,
the effective axial stress resultant for the rod N eff and the two shear components are given by

N eff := t N + M n χ n t + M ν χ ν t , ( 30
)
-T n := 1 t d dS ( t M ν ) + M n χ t - M t t 2 χ n , T ν := 1 t d dS ( t M n ) -M ν χ t + M t t 2 χ ν , (31)
and

M t (L) = m t (L), or φ(L) = 0, -M t (0) = m t (0), or φ(0) = 0, t M n (L) = m n (L), or ωn (L) = 0, ( 32
) -t M n (0) = m n (0), or ωn (0) = 0, t M ν (L) = m ν (L), or ων (L) = 0, -t M ν (0) = m ν (0), or ων (0) = 0.

Tangent equilibrium operator

As illustrated in [START_REF] Simo | The (symmetric) Hessian for geometrically nonlinear models in solids mechanics: Intrinsic definition and geometric interpretation[END_REF], the tangent stiffness operator is defined as the second (covariant) gradient (the Hessian) of the internal energy on the tangent space to the generic configuration of the rod. In general, it is not easy to define a parallel transport operator on the state space, but if it is introduced a Riemmanian's metric by means of the Kozsul formula, it is possible to define the Levi-Civita connection for the configuration space induced by the assumed metric. We stress that the Gateaux derivatives does not represent the components of the (covariant) gradient in the tangent space of the structural model, i.e., the Gateaux derivative is not a tensor and plays the role of the partial derivative in the Riemann's geometry.

Second Gateaux derivative: loss of symmetry

In order to show the noncommutative nature of the second Gateaux derivative, in this section we derive the second Gateaux derivatives of the vectors of the intrinsic triad, which turn out to be nonsymmetric. As it will be evident from the results obtained, the cause of the nonsymmetry stands in the noncommutativity of the rotation operator. This observation will be of great importance in the subsequent developments.

Let φ and φ be any two elements of T ϕ M. Then, the second Gateaux derivatives of the cross section directors are given by

∂ φ ∂ φ t(ϕ) = ∂ φ ∂ φ t(ϕ), ∂ φ ∂ φ n(ϕ) = ∂ φ ∂ φ n(ϕ), ( 33
) ∂ φ ∂ φ ν(ϕ) = ∂ φ ∂ φ ν(ϕ),
where, clearly, symmetry holds only on for the second directional derivative of the unit tangent vector. In particular, the second Gateaux derivative of the unit vector n splits in

∂ φ ∂ φ n(ϕ) = ∂ ū∂ u + ∂ ū∂ φ + ∂ φ ∂ u + ∂ φ ∂ φ n(ϕ) (34)
where the first addend is given by

∂ ū∂ u( n(ϕ)) = d ū dS • t ⊗ n t 2 + n ⊗ t t 2 • d u dS t + - d ū dS • n ⊗ n t 2 • d u dS n + - d ū dS • ν ⊗ n t 2 • d u dS ν. (35)
From the Eq. ( 35), it appears that the components along the vector t and n are symmetric, while the component along ν is nonsymmetric. The remaining terms of the right-hand side of the Eq. ( 35) are given by

∂ ū∂ φ ( n(ϕ)) = - φ t d ū dS • ν t, ∂ φ ∂ u( n(ϕ)) = - φ t d u dS • ν t, (36)
∂ φ ∂ φ ( n(ϕ)) = -φ φ n,
and are symmetric.

The second Gateaux derivative of the unit vector ν splits in

∂ φ ∂ φ ν(ϕ) = ∂ ū∂ u + ∂ ū∂ φ + ∂ φ ∂ u + ∂ φ ∂ φ ν(ϕ) (37)
where the first term is given by

∂ ū∂ u( ν(ϕ)) = d ū dS • t ⊗ ν t 2 + ν ⊗ t t 2 • d u dS t + - d ū dS • n ⊗ ν t 2 • d u dS n + - d ū dS • ν ⊗ ν t 2 • d u dS ν. (38)
From the Eq. ( 38), it appears that the components along t and ν are symmetric, while the component along n is nonsymmetric. The remaining terms of the right-hand side of the Eq. ( 38) are given by

∂ ū∂ φ (ν(ϕ)) = φ t d ū dS • n t, ∂ φ ∂ u( ν(ϕ)) = φ t d u dS • n t, ∂ φ ∂ φ (ν(ϕ)) = -φ φ ν, (39)
and are symmetric.

In the following sections, we define the connection for the tangent space to the configuration space of the structural model considered, i.e., we define the connection for the Kirchhoff rod model; in particular, we adopt the unique torsion-free connection compatible with the assumed metric.

The Koszul formula: the Levi-Civita connection

Let {M, G ϕ } be a Riemannian manifold, equipped with the metric G ϕ : T ϕ M × T ϕ M → R. The fundamental theorem of Riemannian geometry states that there is a unique connection, ∇, (the Levi-Civita connection), which preserves the metric tensor, that is, such that for φ, φ, φ ∈ T ϕ M,

∇ φ G ϕ ( φ, φ) = G ϕ (∇ φ φ, φ) + G ϕ ( φ, ∇ φ φ), ∀ φ, φ, φ, ∈ T ϕ M, ( 40
)
and that is torsion-free, i.e., for all differentiable f :

M → R [ φ, φ] f = ∇ φ φ f -∇ φ φ f, ∀ φ, φ, ∈ T ϕ M. ( 41
)
The unique Levi-Civita connection can be obtained by the Koszul formula, i.e., ∀ φ, φ, φ

∈ T ϕ M 2 G ϕ (∇ φ φ, φ) = ∂ φ G ϕ ( φ, φ) + ∂ φ G ϕ ( φ, φ) -∂ φ G ϕ ( φ, φ) + +G ϕ ([ φ, φ], φ) -G ϕ ([ φ, φ], φ) -G ϕ ([ φ, φ], φ). ( 42
)
Now, we particularize the previous formulae to the Kirchhoff rod model presented. We adopt the following metric:

G ϕ ( φ, φ) := L 0 ū • u + ω( ū) • ω( u) + φ φ dS. ( 43
)
The choice is motivated by the fact that any variation of the position of the centroid line implies a rotation of the directors.

Observing that the metric defined by the definition [START_REF] Luongo | Bifurcation analysis of damped visco-elastic planar beams under simultaneous gravitational and follower forces[END_REF] does not depend explicitly on the configuration of the manifold ϕ, the first three terms on the right-hand side of the Koszul formula trivially vanish so that only the second line of Eq. ( 42) remains. It can be explicitly evaluated using the expression for the Lie parenthesis.

The Lie parenthesis is obtained evaluating the commutator for two independent elements of T ϕ M. Recall that to the variation ũ is associated a translation and a rotation of the directors represented by a skew-symmetric matrix (an infinitesimal rotation superposed to ). The translation commutes with any translation and rotation, so that, indicating by the skew-symmetric matrix associated with any rotation vector ω b and with the skew-symmetric matrix associated with any rotation φ t the commutator reduces to the following:

[ φ, φ] = axial ¯ ˙ + ¯ ˙ -˙ ¯ -¯ ˙ = ωb × ωb . ( 44
)
Using this result in Koszul's formula [START_REF] Luongo | Mode localization in dynamics and buckling of linear imperfect continuous structures[END_REF], and observing that the second and third terms of the second line of the expression vanish for the properties of the triple product, we have

2 G ϕ (∇ φ φ, φ) = [ ωb , ωb ] • ωb = ωb × ωb • ωb . ( 45
)
Therefore, the only non-null structure coefficient is

∇ ωb ωb := 1 2 [ ωb , ωb ] = 1 2 ωb × ωb = ωn ων -ων ωn 2 t. (46)
The first term in Eq. ( 50) is the second directional Gateaux derivative of the energy functional. It splits in the sum of two parts, the first known as the tangent material-stiffness, K M symmetric for an associated material model, and the second known as the tangent geometric-stiffness, G that turns out to be nonsymmetric, if the tangent map is noncommutative. In [START_REF] Simo | The (symmetric) Hessian for geometrically nonlinear models in solids mechanics: Intrinsic definition and geometric interpretation[END_REF], Simo demonstrated that splitting the geometric stiffness into its symmetric and skew-symmetric parts, the latter is equal to the opposite of the second term of the result [START_REF] Porfiri | Circuit analog of beam and its application to multimodal vibration damping, using piezoelectric transducers[END_REF], restoring the symmetry of the second covariant derivative. He proved this property for a Timoshenko rod, using the SO(3) group for the rotations. In this work, we prove that the same conclusion is also true for the configuration space introduced in (12), valid for a Kirchhoff -Love rod, with the metric [START_REF] Luongo | Bifurcation analysis of damped visco-elastic planar beams under simultaneous gravitational and follower forces[END_REF].

With an obvious meaning of the symbols, the tangent stiffness operator, from expression (50), is written as

K tan [ φ, φ](ψ(ϕ)) = K M [ φ, φ](ψ(ϕ)) + sym G[ φ, φ](ψ(ϕ)) + skew[G[ φ, φ](ψ(ϕ))] ∂ φ ∂ φ ψ(ϕ) -∂ ∇ φ φ (ψ(ϕ)). (51)
Having obtained the coefficients of structure, we are able to evaluate the last term in [START_REF] Reissner | On finite deformations of space-curved beams[END_REF]. From the property that the connection is torsion-free, it will be found that skew[G[ φ, φ](ψ(ϕ))] -∂ ∇ φ φ (ψ(ϕ)) = 0, so that the tangent stiffness operator reduces to

K tan [ φ, φ](ψ(ϕ)) = K M [ φ, φ](ψ(ϕ)) + sym G[ φ, φ](ψ(ϕ)) . ( 52
)

Material tangent stiffness: K M

The material tangent stiffness can be evaluated as

M[ φ, φ](ψ(ϕ)) = L 0 E A(∂ φ ε(ϕ))(∂ φ ε(ϕ)) dS + L 0 E I n (∂ φ χ n (ϕ))(∂ φ χ n (ϕ)) dS + + L 0 E I ν (∂ φ χ ν (ϕ))(∂ φ χ ν (ϕ)) dS + L 0 G J t (∂ φ χ t (ϕ))(∂ φ χ t (ϕ)) dS. (53)
Using the results [START_REF] Eremeyev | Thermomechanics of shells undergoing phase transition[END_REF] for the variations of the strain measures, we have the following: material axial stiffness

K ax M [ φ, φ](ψ(ϕ)) = L 0 E A d ū dS • t ⊗ t • d u dS dS (54)
material bending stiffnesses

K ν M [ φ, φ](ψ(ϕ)) = L 0 E I ν d 2 ū dS 2 -φ b • n ⊗ n • d 2 u dS 2 -φ b dS, (55)
and

K n M [ φ, φ](ψ(ϕ)) = L 0 E I n - d 2 ū dS 2 + φ b • ν ⊗ ν • - d 2 u dS 2 + φ b dS. ( 56
)
material torsion stiffness

K t M [ φ, φ](ψ(ϕ)) = L 0 G J t d φ dS + 1 t 2 d ū dS • b d φ dS + 1 t 2 d u dS • b dS. (57)
In the definitions ( 54), ( 55), [START_REF] Simo | On the dynamics in space of rods undergoing large motions. A geometrically exact approach[END_REF], and (57), the vectors t, n, ν, and b are referred to the current configuration ϕ.

Geometric tangent stiffness: G

The formal expression of the geometric tangent stiffness is

G[ φ, φ](ψ(ϕ)) = L 0 E Aε(ϕ)(∂ φ ∂ φ ε(ϕ)) dS + L 0 E I n χ n (ϕ)(∂ φ ∂ φ χ n (ϕ)) dS + + L 0 E I ν χ ν (ϕ)(∂ φ ∂ φ χ ν (ϕ)) dS + L 0 G J t χ t (ϕ)(∂ φ ∂ φ χ t (ϕ)) dS. (58)
Preliminarily, we explicitly calculate the second Gateaux directional derivative of the strain measures that appear in the formula [START_REF] Simo | The (symmetric) Hessian for geometrically nonlinear models in solids mechanics: Intrinsic definition and geometric interpretation[END_REF], starting from the first variations stated in Eq. ( 24) that are recalled here for convenience in a slightly more elaborate form, for an arbitrary element { u, φ} ∈ T ϕ M:

∂ φ ε(ϕ) = d u dS • t, ( 59
)
-∂ φ χ ν (ϕ) = - d 2 u dS 2 + φ t × dt dS • n = - d dS d u dS + 1 t d t dS d u dS + φ t × dt dS • n, (60)
∂ φ χ n (ϕ) = - d 2 u dS 2 + φ t × dt dS • ν = - d dS d u dS + 1 t d t dS d u dS + φ t × dt dS • ν, (61)
∂ φ χ t (ϕ) = d φ dS + 1 t d u dS • t × d t dS . ( 62
)

Second Gateaux derivative of ε

The second derivative of the axial strain-rate measures is given by

∂ φ ∂ φ ε(ϕ) = d u dS • d ū dS , ( 63
)
and this is trivially symmetric.

Second Gateaux derivative of χ ν

The calculation of the second Gateaux derivative of χ ν ((60)) splits in two parts,

∂ φ ∂ φ χ ν (ϕ) = ∂ ū∂ φ χ ν (ϕ) + ∂ φ ∂ φ χ ν (ϕ) ( 64
)
where the first defines the action of the variation ū only, while the second defines the action of the variation of the rotation φ on the strain-rate measure (60) along the generic tangent direction φ ∈ T ϕ M. In particular, one has the following:

∂ ū∂ φ χ ν (ϕ) = - d 2 u dS 2 • t ⊗ n t • d ū dS - d u dS • n ⊗ t t • d 2 ū dS 2 + + d n dS • t d ū dS • n ⊗ n t • d u dS + dν dS • t d ū dS • ν ⊗ n t • d u dS , ( 65
)
and

∂ φ ∂ φ χ ν (ϕ) = φ ν • d 2 u dS 2 + φ ν • d 2 ū dS 2 + φ φ d n dS • t . ( 66
)
Expression (65) is not symmetric, while the second term, given by (66), is symmetric in its arguments. The result (65) is split in its symmetric and in its skew-symmetric parts:

∂ ū∂ φ χ ν (ϕ) = sym ∂ ū∂ φ χ ν (ϕ) + skew ∂ ū∂ φ χ ν (ϕ) (67)
where the symmetric part is given by the sum of the first three terms in the right-hand side of Eq. ( 65), while the skew-symmetric part is defined by

skew ∂ ū∂ φ χ ν (ϕ) = dν dS • t 1 2 d ū dS • ν ⊗ n t 2 • d u dS - d u dS • ν ⊗ n t 2 • d ū dS + = - dν dS • t t • ∇ ω ω.
(68)

Second Gateaux derivative of χ n

The second derivative of χ n is evaluated in a manner analogous to the previous case:

∂ φ ∂ φ χ n (ϕ) = ∂ ū∂ φ χ n (ϕ) + ∂ φ ∂ φ χ n (ϕ). ( 69
)
The two terms are respectively given by

∂ ū∂ φ χ n (ϕ) = d 2 u dS 2 • t ⊗ ν t • d ū dS + d u dS • ν ⊗ t t • d 2 ū dS 2 + - dν dS • t d ū dS • ν ⊗ ν t • d u dS - d n dS • t d ū dS • n ⊗ ν t d u dS , ( 70
)
and

∂ φ ∂ φ χ n (ϕ) = φ d 2 u dS 2 • n + φ d 2 ū dS 2 • n -φ φ dν dS • t . ( 71
)
As in the previous case, the second term (71) is symmetric in its arguments, while the first, given by (70), is not symmetric and is split in its symmetric and in its skew-symmetric parts as follows:

∂ ū∂ φ χ n (ϕ) = sym ∂ ū∂ φ χ n (ϕ) + skew ∂ ū∂ φ χ n (ϕ) , ( 72
)
where the symmetric part is given by the sum of the first three terms of the result (70), while the skew-symmetric part is

skew ∂ ū∂ φ χ n (ϕ) = d n dS • t 1 2 - d ū dS • n ⊗ ν t 2 • d u dS + d u dS • n ⊗ ν t 2 • d ū dS = - d n dS • t t • ∇ ω ω. (73)

Second Gateaux derivative of χ t

Starting from [START_REF] Vidoli | Modal coupling in one-dimensional electromechanical structured continua[END_REF], the second directional derivative of the torsion curvature is

∂ φ ∂ φ χ t (ϕ) = ∂ ū∂ φ χ t (ϕ), ∂ φ ∂ φ χ t (ϕ) = 0. (74)
where

∂ ū∂ φ χ t (ϕ) = - 1 t 2 d t dS • n d ū dS • t d u dS • ν + d u dS • t d ū dS • ν + + 1 t 2 d t dS • ν d ū dS • t d u dS • n + d u dS • t d ū dS • n + + 1 t 2 d u dS • ν d 2 ū dS 2 • n + d u dS • ν d ū dS • t d n dS • t + - 1 t 2 d 2 ū dS • ν d u dS • n + d u dS • n d ū dS • t dν dS • t . ( 75
)
Also, in this case, the first two terms of the Eq. ( 75) are symmetric in the their arguments, while the last two terms are not symmetric and are split in the symmetric and in the skew-symmetric parts, given by

∂ ū∂ φ χ t (ϕ) = sym ∂ ū∂ φ χ t (ϕ) + skew ∂ ū∂ φ χ t (ϕ) (76)
where the first is the sum of the first two terms in the right-hand side of the Eq. ( 75), while the skew-symmetric part, after some calculation, reduces to

skew ∂ ū∂ φ χ t (ϕ) = t • d dS ∇ ω ω . ( 77
)

The tangent stiffness operator for the Kirchhoff rod model

Collecting the results obtained, the geometric tangent operator defined by Eq. ( 58) can be expressed as

G[ φ, φ](ψ(ϕ)) = L 0 N (ϕ)∂ φ ∂ φ ε(ϕ) dS + + L 0 M n (ϕ) sym{∂ φ ∂ φ χ n (ϕ)} + skew{∂ φ ∂ φ χ n (ϕ)} dS + + L 0 M ν (ϕ) sym{∂ φ ∂ φ χ ν (ϕ)} + skew{∂ φ ∂ φ χ ν (ϕ)} dS + + L 0 M t (ϕ) sym{∂ φ ∂ φ χ t (ϕ)} + skew{∂ φ ∂ φ χ t (ϕ)} dS, (78)
The explicit form of the skew-symmetric parts is

L 0 M t t • d dS ∇ ω ω + t × M • b ( t • ∇ ω ω) dS. ( 79
)
Integrating by part, one has

skew G[ φ, φ](ϕ) = L 0 - d M t dS + t × M • b ( t • ∇ ω ω) dS. (80)
The equation in the bracket is the equilibrium equation dual to the torsional velocity of rotation given in Eq. [START_REF] Fox | A drill rotation formulation for geometrically exact shells[END_REF].

In order to obtain the tangent stiffness, we still have to add the last term in Eq. ( 51) that represents the variation of the energy functional in the direction of the coefficients of structure. Since the only non-null connection symbol is given by the expression (47) that has only component along the t direction, then the term -∂ ∂ φ φ (ψ) reduces to

-∂ ∇ φ φ (ψ) = - L 0 - dM t dS + t × M • b ( t • ∇ ω ω) dS. ( 81
)
This term cancels out with (80), and then, the symmetry of the covariant Hessian is recovered.

Conclusions and developments

The consistent tangent stiffness operator for a Kirchhoff -Love rod model has been obtained performing the second covariant gradient of the elastic potential energy, assuming elastic behavior for the cross section. This has been done explicitly evaluating the Levi-Civita connection for the tangent space of the configurations, introducing a suitable metric that accounts for the rotations of the directors. The rotation operator for the Kirchhoff -Love rod has been split in two parts that induce a corresponding split in the rotation on the tangent space of the configurations. The only one coefficient of the defined connection different from zero in the case of Kirchhoff -Love rods is a torsional rotation, defined as the derivative of the bending rotation in the direction of the bending rotation.

The tangent operator evaluated in this way is symmetric for any pairs of admissible variations of the configuration of the rod. The result is essential for developing consistent numerical schemes for the nonlinear analysis of rods, where the consistent symmetric tangent operator is necessary in order to guarantee quadratic rate of convergence, and for buckling analysis.

The strong as well as the weak forms of the equilibrium equations have also been obtained, together with the boundary conditions that define the admissible boundary loads for the nonlinear Kirchhoff -Love beam model. The tangent kinematics has been stated, from which it is possible to obtain the discrete deformation operator in a numerical finite element model.

The results obtained in this paper will be the basis for developing an efficient numerical model for the nonlinear analysis of curved space rods, using the isogeometric interpolation proposed in [START_REF] Greco | B-spline interpolation for Kirchhoff-Love space rod[END_REF][START_REF] Greco | An implicit G 1 multi patch B-spline interpolation for Kirchhoff love space rod[END_REF]. The efficiency of the method, already examined for the linear case, will be tested when large deformations are attained, causing degradation of the properties of the tangent stiffness matrix. An interesting problem that will be investigated is the applicability of the method presented to the case of nonconservative forces and to the related buckling analysis of systems of rods with perturbation methods, as done in [START_REF] Luongo | Bifurcation analysis of damped visco-elastic planar beams under simultaneous gravitational and follower forces[END_REF][START_REF] Luongo | Double zero bifurcation of non-linear viscoelastic beams under conservative and nonconservative loads[END_REF][START_REF] Paolone | Flexural-torsional bifurcations of a cantilever beam under potential and circulatory forces: Part I Nonlinear model and stability analysis[END_REF][START_REF] Paolone | Flexural-torsional bifurcations of a cantilever beam under potential and circulatory forces: Part II Post-critical analysis[END_REF], where the stability of a cantilever elastic beam under the action of a follower tangential force and a bending conservative couple at the free end was analyzed. Similar issues arise in the analysis of deformable curved beams under the action of wind loads.

Applying the vector [START_REF] Paolone | Flexural-torsional bifurcations of a cantilever beam under potential and circulatory forces: Part I Nonlinear model and stability analysis[END_REF] to the unit triad { t, n, ν}, we get the skew-symmetric part of the second variation of the unit triad that, as evidenced by Eqs. [START_REF] Jelenic | A kinematically exact space finite strain beam model-finite element formulation by generalized virtual work principle[END_REF], [START_REF] Luongo | Multiple interaction and localization phenomenon in postbuckling of compressed thin-walled members[END_REF], is associated with the bending rotations only. This fact is graphically represented in Fig. 1 that shows two subsequent bending rotations around two orthogonal axes occurring in reversed order a → b → c and a → d → c. The red arrow denotes the tangent to the rod axis, and the blue or green arrow denotes the normal vector. The final configurations differ only by a torsional rotation.

Second covariant gradient: the Hessian

We are now in position for evaluating the second covariant gradient of the internal energy potential ψ int (ϕ). Differentiating the first gradient of ψ int , given by

we have the following:

where the apex indicates the contravariant basis. With simple algebra, we get: )


Joseph E Salah

Marie-Louise Duboin

Christine Mazaudier

Ionospheric electrodynamics over Saint

published or not. The documents may come   L'archive ouverte pluridisciplinaire


Creation of a BI environment from scratch with open source software.

A practical case -Thierry Dautcourt 1 1 Inria, 54600 Villers-lès-Nancy, Thierry.Dautcourt@inria.fr This presentation examines a practical case and the methods used to build a BI (Business Intelligence) environment, in a Research's Institute lacking initial BI maturity. We describe here some conditions for success and the choice of open source software to build the technical BI environment.

Context

The French Institute for Research in Computer Science and Automation (INRIA) is a public research body fully dedicated to computational sciences. It is organized into 210 research teams split across eight Inria centres. In 2010, our management requested our own BI environment. The main purpose of the BI environment was to produce « automatically » the major performance indicators related to our quadrennial contract signed with the Government (Ministry). At the time, these major indicators were produced by aggregating heterogeneous spreadsheet files and with a lot of manual work for collecting, aggregating and reconciling the data.

How to proceed ? Where to begin ? The field of contract indicators covers a very large domain from human resources, and financial to scientific production. In addition, the quality of data was uneven and the culture of data management rather low. The approach was to work on three axes: identifying one strategic indicator, finding people to work with and building a technical environment for people to access analytic data.

Focus on one strategic indicator

The first axis was to focus on one major indicator among many. We choose one whose values were not computed. This indicator is the rate of deposit 1 in the open access system of scientific publications produced by researchers. As Inria is strongly involved in the «Berlin Declaration 2 » the goal was to provide full open access to knowledge in the Sciences and Humanities. For this measure, we relied on two sources of data: the publications present in the open archive system (HAL 3 ), and the existing publications in the team's annual report activity 4 .

Finding the people to work with

The second axis was to find « open minded people » to work with. The subject of publications was an opportunity, as people working in scientific and technical Information are usually working in a network and are very involved in data quality and data cleaning processes.

Create a BI environment with open sources

The third axis is the technological environment. We didn't have a budget for software licences. Instead we used three open source software packages, Talend Data Integration for the ETL (Extract, Transform and Load) and SpagoBI for the BI web server. Data are stored in MySQL databases.

Talend Data Integration 5 (community version) was used for building the data warehouse. This tool is very accurate for incremental conception. The different phases of transformation are mostly done without any lines of programming but with connecting high-level components in sequence. However, this does not eliminate method. Evolution over time is easy to manage and costs of modification are small (1 or 2 days), depending on the size of modification.

The ETL is used for extracting, homogenising, denormalizing, and matching data in a data warehouse. We must also manage errors issued by harmonizing data of two sources, an open archive and an annual report. The goal was to build a multidimensional dataset (cube) on subject "publications". SpagoBI 6 is the second tool for providing a web access to data (and errors for correction), for exploring and consulting analytical views. This is a full open source Business Intelligence platform that uses many tools as analytical engines. Our first major interest was about OLAP (On-Line Analytical Processing) cube.

The dimensions used for analysis are the organizational dimension (team project and research centre), time dimension (year) and other dimensions specific to the subject area, such as the typology of publications dimension. The measure, rate of deposit, is computed accounting for the overlap between publications in the team's annual activity report and publications in open access archive HAL.

The others conditions to succeed (and not least!)

The platform allows for monitoring the deposit of publications over time, in open archive. But this is not enough. Additional effort focused on two directions: 1) Services to researchers for providing a dynamic html frame on publications ready to integrate in team or personnel web pages 7 ; 2) An accompanying program from scientific and technical information staff, helping researchers to deposit publication in the HAL archive. The indicator helped highlight "low deposit teams", and to uncover the reasons for their poor commitment.

For many years, to support research teams in the production of their annual activity report, a document skeleton is produced using data from our IS (information system). For the 2013 annual activity report, for the first time, all publications were automatically fed in the team annual report skeleton. This is also the case for other data. Each year, more data are added automatically from the data warehouse in the "activity report skeleton" such as team participants, funding's, International, European and National relationships with Industry and University or Research Centre.

In 2014, our BI environment integrates about 12 analytics cubes and 40 views. Increasingly, the data warehouse became a key part of our IS. The business approach and agile construction of data warehouse help the dialog between business requestors and IT people for building IS. The exposition and reuse of data by a larger population increases the robustness of the Institute's data and along the way, quality indicators.

  Figure 1 presents the measure by two dimensions: years and research organisation, here INRIAs 8 research centers.


G Marceau-Caron
email: gaetan.marceau-caron@inria.fr

M Schoenauer
email: marc.schoenauer@inria.fr

Racing Multi-Objective Selection Probabilities

Keywords: Multi-Objective Evolutionary Optimization, Hoeffding Races, Uncertaintly Handling, Noisy Multiobjective Optimization

In the context of Noisy Multi-Objective Optimization, dealing with uncertainties requires the decision maker to define some preferences about how to handle them, through some statistics (e.g., mean, median) to be used to evaluate the qualities of the solutions, and define the corresponding Pareto set. Approximating these statistics requires repeated samplings of the population, drastically increasing the overall computational cost. To tackle this issue, this paper proposes to directly estimate the probability of each individual to be selected, using some Hoeffding races to dynamically assign the estimation budget during the selection step. The proposed racing approach is validated against static budget approaches with NSGA-II on noisy versions of the ZDT benchmark functions.

Introduction

Uncertainty handling is an important aspect of optimization since it concerns most, if not all, real-world applications. Optimizing uncertain objectives aims at taking into account modeling inaccuracies, measurement errors from sensors, or prediction errors, that will interfere with the beliefs of a decision maker about the environment. Therefore, optimization under uncertainty must include some mechanisms that ensure, one way or another, that the proposed solutions are effective, according to the user's point of view w.r.t. optimality. And whereas several definitions of such effectiveness can occur in the simplest case of a single objective, the complexity of optimizing multiple uncertain objectives increases drastically with the number of objectives.

The general framework of this work is that of multi-objective optimization in uncertain context. The degrees of freedom of the decision maker are the variables of the optimization problem, defined on the decision space, and observables are some responses of the system when setting these variables. However, the same setting will result in different responses every time it is used, and the output of the system thus defines a probability distribution, conditionally dependent on the decision variables. The goal of the optimization process is then to find the values of the decision variables that will optimize some statistics on this probability distribution. The choice of these statistics depends on the user's goal and preferences. The average or the median are common choices, though probably sometimes only because of the lack of efficient methods to handle other statistics. For instance, risk-adverse users will prefer to minimize the consequences of the worst outcomes, while risk-affine users will maximize their possible "profit" even if it comes at high risk, optimizing the value at risk for a given risk level.

Except when the type of noise is known -a totally unrealistic hypothesis-a common way to compute the desired statistics is to sample the fitness of each individual as many times as necessary to obtain a good estimation thereof, and the amount of computation per individual is user-defined, uniformly over the individuals and the generations. In the single-objective framework, an alternative has been proposed, using the idea of races [START_REF] Heidrich-Meisner | Hoeffding and Bernstein Races for Selecting Policies in Evolutionary Direct Policy Search[END_REF], minimizing the number of reevaluations while keeping a high confidence level on the results. But only limited attempts have been made in the multi-objective framework (see Section 2).

The approach proposed in this paper, Racing Selection Probability (RSP), is an attempt, in the multi-objective case, to dynamically decrease the number of sampling of all individuals by applying the principles of Hoeffding races directly on the estimation of the probability of being selected for an individual, using bounds on the behavior of that probability to decide as early as possible when to definitely select or discard an individual, for a given confidence level. Bounds on any statistics can be used, and straightforwardly embedded in any unmodified EMOA, thus allowing to handle any preference of the user. Furthermore, any type of noise can be handled that way.

The paper is organized the following way. Section 2 briefly surveys state-ofthe-art methods for uncertainty handling in Evolutionary Multi-Objective Optimization. Section 3 introduces the Hoeffding's inequality used in RSP. Section 4 presents experimental results on perturbed ZDT test functions (with different types of noise) where RSP is compared to the two basic noise-handling methods, the implicit and static averaging. More details and results are available in the corresponding Technical Report3 .

Uncertainty Handling in Multi-Objective Evolutionary Optimization

The context of this work is that of Multi-Objective Optimization with Uncertainty. On the space of decision variables X , several conflicting objectives f 1 , . . . , f k are defined (to be minimized, w.l.o.g.), and, as discussed above, the outcome of any given setting of the decision variables is a probability distribution f over the objective space F ⊂ R k that depends on the values of the variables and on some additional unknown external random variable ε, aka noise. In particular [START_REF] Basseur | A Preliminary Study on Handling Uncertainty in Indicator-Based Multiobjective Optimization[END_REF], there is no "true" value of the objectives to which some random noise is added. Formally [START_REF] Trautmann | Pareto-Dominance in Noisy Environments[END_REF][START_REF] Voß | New Uncertainty Handling Strategies in Multiobjective Evolutionary Optimization[END_REF], the Multi-Objective Noisy Optimization Problem (MNOP) can be written as

min x∈X ((f |x, ε) = (f 1 , . . . , f k |x, ε)) (1)
where f is a random variable taking values in F, and each f i is a real-valued random variable, a coordinate of f .

Even in the single objective case, the minimization of a random variable does not make much sense. So the user must complete the problem definition by providing some preferences through some statistics over that random variable (e.g., minimizing the mean, the median, the 5% percentile, the variance with constraint on the mean, . . . ). The situation is the same in the multi-objective case, except that there doesn't exist any total order on the samples of the random variable of interest. In the deterministic case, Pareto dominance has proved useful, and the notion of Pareto front is accepted as a way to describe interesting solutions of the multi-objective problem at hand. In particular, several multi-objective optimization algorithms have been proposed, among which Evolutionary Multi-Objective Algorithms (EMOAs) (see e.g., [START_REF] Zhou | Multiobjective evolutionary algorithms: A survey of the state of the art[END_REF]). And because uncertainty is ubiquitous in real-world problems, MNOPs have also been well studied, though not always with such a degree of generality.

Previous Work

A first approach is to port to multi-objective context the single-objective static averaging techniques, that re-evaluate every individual N times at each generation (also called implicit averaging if N = 1).

Several works consider the specific case of additive noise of known type: the random variable (f |x, ε) is of the form g(x) + ǫ for some function g(x) and some partly known noise ǫ. Depending on the form of ǫ, approximation of the probabilistic dominance (probability that an individual Pareto-dominates another one) can sometimes be computed at low computational cost. In [START_REF] Teich | Pareto-Front Exploration with Uncertain Objectives[END_REF], the noise is supposed bounded, and exact calculations are done for uniform noise; In [START_REF] Hughes | Evolutionary Multiobjective Ranking with Uncertainty and Noise[END_REF], the noise is assumed Gaussian with known variance (that can be computed off-line from static samples). This work is extended in [START_REF] Fieldsend | Multi-Objective Optimisation in the Presence of Uncertainty[END_REF] to the case of unknown (and non-uniform) variance. Later, [START_REF] Eskandari | Handling uncertainty in evolutionary multiobjective optimization: SPGA[END_REF] proposed another way to compute the probability with more general hypotheses, but going back to using a fixed number of samples (15 in experiments). In any case, it is clear that the hypothesis of a known type of noise is highly unrealistic in practice.

An approach that is specific to indicator-based algorithms is proposed in [START_REF] Basseur | A Preliminary Study on Handling Uncertainty in Indicator-Based Multiobjective Optimization[END_REF], that does not make any hypothesis about the noise and uses the general model of Equation 1: the indicator ε + is approximated using averages (over 5 samples), and is used within the environmental selection. However, the problem being solved there is the minimization of the expectation of the indicator at hand (w.r.t. some reference set), that cannot be adapted to the user's preferences.

Several works propose different approaches to probabilistic dominance for the general MNOP (equation 1). Pareto Dominance in Uncertain environments(PDU) [START_REF] Trautmann | Pareto-Dominance in Noisy Environments[END_REF] uses the convex hull of a fixed number of samples (10 in the paper) to estimate both the mean and its uncertainty. In [START_REF] Voß | New Uncertainty Handling Strategies in Multiobjective Evolutionary Optimization[END_REF], PDU evaluates the certainty of the mean using quartiles on each dimension, and some races are run for each objective, from [START_REF] Heidrich-Meisner | Hoeffding and Bernstein Races for Selecting Policies in Evolutionary Direct Policy Search[END_REF], with confidence 0.0001 and maximum race length 15. This latter work however assumes that the noise distribution is symmetrical, and Suzuki and collaborators propose another Pareto Dominance operator that does not need that hypothesis [START_REF] Boonma | A Confidence-based Dominance Operator in Evolutionary Algorithms for Noisy Multiobjective Optimization Problems[END_REF] using a CPU-expensive SVM construct over the samples; [START_REF] Phan | A Non-parametric Statistical Dominance Operator for Noisy Multiobjective Optimization[END_REF] improves the method using a non-parametric Mann-Whitney U-test. However, both works use a fixed number of samples (resp. 30 and 20) to estimate the dominance operator.

In [START_REF] Siegmund | Sequential Sampling in Noisy Multi-Objective Evolutionary Optimization[END_REF], six different resampling approaches are compared. All but one use some absolute criteria that only depend on some statistics on the previous samples and the individual at hand to decide on early stop of the resampling procedure and derive an estimation of the mean of the sample with known confidence. That mean is then used as the fitness in a standard EMOA. The last procedure (termed OCBA) is the closest to RSP proposed here, in that it makes the minimal global sampling allocation to estimate the confidence in a partition of the population into a non-dominated and a dominated sets. However, the calculation of the confidence assumes Gaussian noise on all objectives.

Discussion, and Rationale for RSP

Our goal is to design, within a given EMOA, an approach that will limit the number of resampling while preserving some confidence on the resulting Paretobased selection, for a wide range of statistics describing the user's preferences, and without any requirement on the type of noise. Most of the works listed above, however, use a fixed user-define resampling budget (except [START_REF] Siegmund | Sequential Sampling in Noisy Multi-Objective Evolutionary Optimization[END_REF] and [START_REF] Voß | New Uncertainty Handling Strategies in Multiobjective Evolutionary Optimization[END_REF]). Furthermore, either they derive estimations of the mean of a sample with some confidence interval -and this does not allow to derive confidence bounds on the comparison between those means (except in specific cases, e.g. Gaussian distributions); or they do derive probabilistic Pareto dominance, with known confidence, but omit the second component of Pareto-based selection, the diversity preserving mechanism (the case of indicator estimation [START_REF] Basseur | A Preliminary Study on Handling Uncertainty in Indicator-Based Multiobjective Optimization[END_REF] is different, but strictly limited to indicator-based EMOAs).

The idea of RSP borrows from [START_REF] Heidrich-Meisner | Hoeffding and Bernstein Races for Selecting Policies in Evolutionary Direct Policy Search[END_REF], like [START_REF] Voß | New Uncertainty Handling Strategies in Multiobjective Evolutionary Optimization[END_REF] cited above, is using Hoeffding races 4 to decrease the number of resampling while nevertheless guaranteeing some level of confidence on the statistic at hand. But contrary to the works above (including [START_REF] Voß | New Uncertainty Handling Strategies in Multiobjective Evolutionary Optimization[END_REF]), Racing Selection Probability, as its name claims, will perform the race on the probability of an individual to be selected by the selection mechanism of the chosen EMOA.

Racing Selection Probability

Let us assume some selection procedure in an existing MOEA (e.g., non-dominated sorting + crowding distance for NSGA-II [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF]) that aims at selecting µ individuals out of a population of size λ. The basic idea of RSP, inspired by [START_REF] Heidrich-Meisner | Hoeffding and Bernstein Races for Selecting Policies in Evolutionary Direct Policy Search[END_REF], is to estimate, for any individual i, with as few samples of fitnesses as possible thanks to Hoeffding bounds, the probability p sel i that i will be selected.

Hoeffding's inequality states that, for any random variable X with range width R, and for any confidence level 1 -δ, the absolute difference between the expectation and the empirical mean computed using t samples is upper-bounded by R log(2/δ)/2t.

Every time all λ individuals are resampled, the standard selection procedure of the EMOA at hand is applied to the current sample, determining the selected µ ones. This results in a new sample for probabilities p sel i . For any δ, lower and upper values for all p sel i can be computed at confidence level 1 -δ, thanks to Hoeffding's bound applied to p sel i . Any individual i whose lower bound for p sel i is larger than the upper bound of at least λ -µ other p sel k is definitely selected and leaves the race. Symmetrically, any individual i whose upper bound for p sel i is smaller than the lower bound of at least µ other p sel k is definitely discarded and leaves the race. Remain in the race the uncertain individuals, and only those are resampled again at next iteration. The race ends when either µ individuals are definitely selected, or λ-µ individuals are definitely discarded, or some maximum number of resampling T M ax have been done. In the latter case, selection is made without Hoeffding guarantee. Nevertheless, the proposed procedure allows to quickly select the most promising individuals with given confidence. Note that each selection step can also be done on some statistic for each individual given the past t samples. This will be illustrated in Section 4 where variants using the average or the median will be used, instead of the most recent sample. These variants will be termed RSP AV G and RSP M ED respectively, the variant that does not use any statistic being denoted by RSP I .

There are however some specificities to the multi-objective context. First, the selection step usually involves the whole population, be it indicator based, or using some diversity secondary criterion. Hence the individuals that have left the race should nevertheless be taken into account for the next selection steps -but without being themselves re-evaluated, of course. A bootstrap procedure is used here, to mimic an ever growing sample without any resampling. Note that bootstrap is also used to avoid resampling individuals that have not been modified by the variation operators from one generation to the next.

Finally, it might be beneficial to detect early that some race will not end before the maximum number of samples because of actual ties between individuals that remain in the uncertain set. Here, when the sum of absolute pairwise differences of the empirical mean of the p sel i becomes lower than a given threshold called Proximity Threshold, the race stops and the µ best individuals according to the current selection policy are returned.

Experimental Results

Experimental Conditions

Five methods have been experimentally compared: the implicit averaging, and two variants of the static sampling, whether the average or the median of the samples is used for the selection (see Section 2.1); and 3 variants of RSP, whether The testbench is based on the classical ZDT suite, used either as is (deterministic setting), or with known additional noise: Gaussian noise, that should favor the average estimator compared to the median estimator, the empirical average being the minimum-variance unbiased estimator of the expectation of a normal distribution with unknown mean and variance; Cauchy noise, that has an infinite mean, hence the mean estimator should be perturbed because of the outliers; and Gumbel noise, an asymmetrical distribution with finite moments that is used in extreme value theory to simulate rare events (its location parameter is chosen in order to center the median).

The goal of the experiments is to study the impact of the two parameters Sampling Budget and Confidence Level, and possibly their interaction, e.g., if the required confidence level is too high, all races will reach the maximum budget, and RSP amounts to static sampling. All parameter values (for the algorithms and the noise models) that have been used for these experiments are listed in Table 1). All runs were limited to 100k evaluations, except ZDT6 (500k), and 25 independent runs were run for each parameter setting. Due to space limitation, results on ZDT2 are not shown, but are quite similar to those of ZDT4.

Results

The performances are compared using the difference hypervolume indicator w.r.t. the real Pareto front, on the normalized objective space. The normalization is done with respect to the Nadir point, computed from the union of the exact Pareto front and every point generated by each algorithm for a given function and a given noise. Statistical significance is attested by p-values of the Wilcoxon signed-rank test. Pisa performance assessment tools (http://www. tik.ee.ethz.ch/sop/pisa) was used to compute the hypervolumes.

Each plot of the following figures summarizes the results obtained by all algorithms on one function with one type of noise (or no noise at all): each plot displays several boxplots, each boxplot represents the statistics of the 25 hypervolume values at the end of each of the 25 runs for the corresponding setting. Each plot is divided into six regions. First boxplot is that of the implicit averaging I. Next 2 regions give the results of the static sampling (resp. AV G and M ED), and display 6 boxplots each, corresponding to the 6 Sampling Budget values of Table 1. Next 3 regions give the results for RSP I , RSP AV G and RSP M ED resp. For each region, there are 6 subregions (the 6 values of Sampling Budget) with two boxplots each, one for each confidence level (25%, 95%).

Discussion

First of all, in the deterministic case, the results of implicit averaging assesses that the total budget of 100k samples is sufficient for NSGA-II to find a good approximation of the Pareto front. Furthermore, as expected, the performance of static resampling using an estimator degrades with the Sampling Budget, as more and more samples are wasted on the (sometimes useless) estimation of the statistic. In the same situation, RSP is able to detect the low (!) uncertainty and to stop the race early, at least when using a Confidence Level of 25%. A Confidence Level of 95% can sometimes, on the other hand, lead to a similar degradation than in the static setting. The anomalies in that respect for RSP on ZDT3 (discontinuous front) for small Sampling Budget (5 and 10) might come from races that stop too early with all selected individuals in the same component of the front. On the noisy instances, implicit averaging does not perform very well compared to the other uncertainty handling approaches. Surprisingly, even if the medians are higher, the spread of the performances is not greater than the other approaches, excepted for ZDT4-CA. It can be due to the fact that without any uncertainty handling approach, the probability that every individual of the population is good or bad is small and so, at the population level, the performance does not vary so much from one run to another.

Beside, implicit averaging is comparable to AVG in case of Cauchy noise, for all functions but ZDT4: choosing by default the mean (a common choice) can lead to poor results when the distribution of the noise is unknown. Using RSP seems to mitigate this effect, probably because it uses the probability of survival instead of the estimator of the mean.

Comparing, for each noisy function, the best configurations of RSP and static sampling leads to the following considerations: the results are statistically equivalent for all cases of noisy ZDT2 and ZDT4; RSP is significantly better (p-value < 10 -5 ) than static sampling in 5 cases (the 3 noisy ZDT6, and ZDT1 and 3 with Cauchy noise), is slightly worse (p-value in [0.01, 0.1]) in 2 cases (both ZDT1 and 3 with Gaussian noise), and both approaches are equivalent (p-value > 0.1) on the remaining 2 cases (both ZDT1 and 3 with Gumbel noise). On ZDT1 and 3 with Gaussian noise, static sampling with averaging performs best: this is most probably related to the fact that AVG is based on the minimum-variance unbiased estimator, while RSP AV G uses it indirectly to estimate the probability of survival.

Regarding the choice of estimator, racing seems to decrease the impact of the average vs median issue. Indeed, when using static sampling, average performs slightly better than median for Gaussian and Gumbel noises, whereas median is consistently and significantly better when facing Cauchy noise. On the opposite, all 3 variants of RSP perform in general similarly over all problems. In particular, the no-estimator, RSP I , performs as good as both others on most problems. This is good news, as it gives hope that the proposed racing approach might perform well with a lot of estimators, allowing the user to actually choose his favorite without having to care about the optimization algorithm in that respect.

Conclusion and Perspective

RSP is a general approach to uncertainty handling in existing EMOAs. It uses a (µ, λ) Hoeffding race at a given confidence level, inspired by [START_REF] Heidrich-Meisner | Hoeffding and Bernstein Races for Selecting Policies in Evolutionary Direct Policy Search[END_REF], though applied directly on the selection probabilities of the individuals in the population. It is agnostic w.r.t. the selection method, and hence can accomodate any user preference that could be carried by the algorithm selection.

First experimental results within NSGA-II on noisy versions of ZDT benchmarks, indicate that this path is worth following for future research: RSP performs significantly better than implicit averaging or static sampling in many situations, and never performs significantly worse. It is less sensitive to the Sampling Budget parameter, especially for small (on zero) levels of noise, and surprisingly almost insensitive to the choice of the estimator. On the other hand, it is very sensitive to the Confidence Level of the races. However, these partial conclusions should be sustained by deeper analyses and validated by more experiments, with different levels of non-homogeneous noise, and other test functions from real-world problems.

The main perspectives for further work are to couple RSP with other EMOAs such as SPEA-2, IBEA and HYPE, in order to study the interaction between racing and the indicator function. Also, RSP should also be compared to more sophisticated uncertainty handling methods (see Section 2). It is also mandatory to test other estimators within RSP, as well as different noise models and different noise intensities. On the more fundamental side, it should be possible to better understand the intricate relationship between estimating the selection probability and directly estimating the objective values.

A longer term research track is to come up with some adaptive procedure to dynamically tune the Sampling Budget and, maybe more importantly, the Confidence Level. Indeed, it is clear from the present results that this latter parameter has a strong effect on the performance of the algorithm and should be fixed carefully. In the case where its optimal value varies over the decision space, only adaptive tuning can perform well on most functions. To conclude, we feel that the use of RSP in EMOA is a promising avenue for taking into account the decision maker's preferences and increasing the reliability and robustness of the solutions in many real-world applications.

Fig. 1 .

 1 Fig. 1. Results for ZDT1 (4 top plots) and ZDT3 (4 bottom plots). See Section 4.2.

Fig. 2 .

 2 Fig. 2. Results for ZDT4 (4 top plots) and ZDT6 (4 bottom plots). See Section 4.2.

Table 1 .

 1 Parameters for (top to bottom) benchmarks and noise; static sampling; RSP the last sample, the average or the median of the previous samples are used in the selection (see Section 3). All RSP variants have been implemented within NSGA-II with standard SBX crossover and polynomial mutation. A common parameter of static sampling and RSP is the Sampling Budget, that will denote the fixed number of samples for each individual in the static case, and the maximum length of the races in RSP. RSP also requires a Confidence Level and the Proximity Threshold (see previous Section 3).

	ZDT Functions:	{1,2,3,4,6}	Number of runs:	25
	Deterministic(DE): Dirac Delta Function	Gaussian noise(GA): 0.25 * N (0, I)
	Cauchy Noise(CA): (0, 0.25)	Gumbel noise(GU): (2, 2 ln(ln(2))
	Population Size:	100	Nb Eval.:	{100k,500k}
	SBX Crossover:	pc = 1.0, η = 20	Polynomial Mut.:	pm = 1/|X |, η = 20
	Confidence Level: {0.25, 0.95}	Proximity Thres.:	0.5
	Sampling Budget: {5, 10, 15, 20, 30, 50}	Estimators: {None, AVeraGe, MEDian}

http://hal.inria.fr/hal-01002854/en

[START_REF] Heidrich-Meisner | Hoeffding and Bernstein Races for Selecting Policies in Evolutionary Direct Policy Search[END_REF] also advocates Bernstein races when the range of values is not know -which is not the case here. Hence Bernstein races will not be mentioned here.


Jan Margeta

Antonio Criminisi

Daniel C Lee
email: dlee@northwestern.edu

Nicholas Ayache
email: nicholas.ayache@inria.fr

Recognizing cardiac magnetic resonance acquisition planes

In this paper we propose a method for automatic wrangling of missing or noisy acquisition plane information of cardiac magnetic resonance images in order to simplify case filtering and image lookup in large collections of cardiac data. To recognize standard cardiac planes we use features based on image miniatures combined with a decision forest classifier. We show that augmenting the dataset with a set of nondestructive transformations can improve classification accuracy. Our approach compares favorably to the state of the art while requiring fewer manual annotations.

Introduction

Cardiac magnetic resonance image acquisitions are mostly performed along several oblique axes with respect to the structures of the heart instead of traditional body planes (coronal, axial and sagital) to better visualize and quantify different aspects of the heart. Knowledge of these planes is also essential for data collection organization, selecting adapted image processing algorithms, grouping of related slices into volumetric image stacks, filtering of cases based on completeness for a clinical study and to respond with the most relevant acquisition plane in content based image retrieval. This orientation information is often in some way encoded within the DICOM image tags: Series Description (0008,103E) and Protocol Name (0018,1030). The way to encode this view information is however not standardized, operator errors are present or such information is often completely missing. Searching through large databases to handpick desired views from cases in the collection is therefore tedious. Previous work in this field has concentrated mainly on real-time recognition of cardiac planes in echography acquisitions. These methods are based on appearance models, registration or part detection and require additional information to train view, part [START_REF] Park | Automatic Cardiac View Classification of Echocardiogram[END_REF] or landmark specific detectors [START_REF] Zhou | Automatic view classification for cardiac MRI[END_REF]. Therefore c 2014. The copyright of this document resides with its authors. It may be distributed unchanged freely in print or electronic forms. any new view will require these extra annotations to be made. Otey et al. [START_REF] Otey | Automatic view recognition for cardiac ultrasound images[END_REF] on the other hand addressed this problem by using view labels only from gradient based image features.

Cardiac planes of acquisition

Optimal cardiac planes depend on global positioning of the heart in the thorax. This is more vertical in young individuals and more diaphragmatic in elderly. Imaging in standard cardiac planes ensures efficient coverage of relevant cardiac territories and enables comparisons across modalities, thus enhancing patient care and cardiovascular research. A good overview of standard cardiac acquisition planes can be found in [START_REF] Taylor | Cardiovascular MR Imaging Planes and Segmentation[END_REF].

Left ventricular long axis acquisition planes. These are usually acquired as 2D or cine 2D+t stack. The 2-chamber, 3-chamber, and 4-chamber views (figs. Short axis acquisition planes. Short axis slices (fig. 1(e)) are oriented parallel to the mitral valve ring. These are acquired regularly spaced from the cardiac base to the apex of the heart, often as a cine 3D+t stack. These views are excellent for reproducible volumetric measurements or radial cardiac motion analysis but their use is limited in atrio-ventricular interplay or valvular disease study.

Method

We propose automatic cardiac view recognition pipeline (Fig. 2(a)) that learns to recognize the acquisition planes directly from cardiac magnetic resonance images by combining image miniatures with classification forests. These are trained in a supervised fashion to automatically select the relevant features. To augment the training dataset we jitter the images with small translations, rotations and scales. Apart from the view annotation of the training set, no other information is necessary to classify the views.

Slice normals from DICOM

For images where the DICOM orientation (0020,0037) tag is present we can use it to predict the cardiac image acquisition plane. Similarly to [START_REF] Zhou | Automatic view classification for cardiac MRI[END_REF] we extract image acquisition plane normal vectors as a cross-product of the two image orientation vectors specified in the tag and use this three-dimensional vector feature vector to describe each image (See Fig. 2(b)).

Image-based features

When the orientation tag is not present we rely on the image intensity information only. To capture differences between planes of acquisition, we extract a set of image based features.

To allow use of this for method for all of the above mentioned applications, we use 2D image slices individually rather than 3D or 3D + t volumes.

Preprocessing Our only preprocessing steps are image central square extraction, resampling to 128x128 pixels and linear rescaling of the intensities to range between 0 and 255.

Pooled image miniatures Radiological images are mostly acquired with the object of interest in the center. Therefore some rough alignment of structures can be expected (See fig. 3). Image intensity samples at fixed locations can therefore provide strong clues about the position of different tissues (e.g. symmetric dark lungs or bright cavity in the center).

It has been shown before [START_REF] Torralba | 80 Million Tiny Images: a Large Data Set for Nonparametric Object and Scene Recognition[END_REF] that significantly down-sampled miniatures can be used

2CH 4CH 3CH
LVOT SAXAp SAXA SAXM SAXB SAXSB for such image recognition. In this paper we subsample the cropped centers to two fixed sizes (20x20 and 40x40 pixels). We also create another set of pooled image miniatures from each cropped center by dividing the image into non-overlapping 4x4 tiles and computing intensity minimum and maximum in each of these tiles (fig. 2(c)). As in [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF], pooling can add invariance to small image translations and rotations (whose effect is within the tile size). We used directly pixel values from these miniatures sampled at random positions as features.

Each image can be seen as a data point in this high-dimensional feature space.

Decision forest classifier

Random decision forest [START_REF] Breiman | Random forests-random features[END_REF] is an ensemble classification method that consists of a set of binary decision trees. This method is computationally efficient and allows automatic selection of relevant features for the prediction. The tree structure is optimized by recursively partitioning the collection of data points X into the left or right branches such that points with different labels get separated while the same label points are grouped together. At each node of the tree a feature from a randomly drawn subset of all features is chosen such that impurity I in both branches is minimized. We weight samples from the under-represented classes more and we set sample weights inversely proportional to class probabilities and normalize them such that sum w le f t + w right to equal to one at each node.

I(X, θ ) = w le f t H(X le f t ) + w right H(X right ) (1)
H is weighted entropy and X le f t and X right are point subsets falling to either the left or the right branch, based on the tested feature value and threshold. Similarly, w le f t , w right are sums of sample weights at each branch and w c is sum of weights for a particular class c.

H(X) = -∑ c w c log(w c ) (2)
Only a random subset of features (i.e. single pixel values at 64 different locations of the miniatures) is tested at each node of each tree and a simple threshold on this value is used to divide the data points into the left or the right partition. This helps to make the trees in the forest different from each other which leads to better generalization.

When classifying a new image, features chosen at the training stage are extracted and the image is passed through the decisions of the forest (fixed in the training phase) to reach a set of leaves. Class distributions of the reached leaves are averaged across the forest and the most probable label is selected as the image view.

Augmenting the dataset with jittering

Our dataset contains in total 960 image slices from 100 cardiac patients (SAX:540, 4CH: 140, 2CH: 112, 3CH: 107, LVOT: 9). This might seem rather modest compared to recent image recognition work e.g. [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] where deep models are trained with 15 million images. The object of interest is in general placed in the image center, however deviations exist. To account for these differences we expand the training set with extra images by applying a set of transformations of the original training images. Similarly to [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] these are translations (all shifts in x and y between -10 and 10 pixels for a 5x5 grid), but we also add rotations (angles between -10 and 10 degrees with 20 steps) and scales (1-1.4 zoom factor with 8 steps while keeping the same image size) resulting in a total of 51894 training images. Note that this extra expense is mainly at the training time. The test time is almost unaffected except that now a deeper forest can be learnt.

Validation

Datasets We validated this method on a dataset of 100 patients from multi-center study on post myocardial infarction hearts DETERMINE [START_REF] Kadish | Rationale and design for the Defibrillators to Reduce Risk by Magnetic Resonance Imaging Evaluation (DE-TERMINE) trial[END_REF]. These consist of steady state free precession magnetic resonance acquisitions.

Results

We ran a 25-fold cross validation by dividing the dataset on patient identifier to prevent biasing our results. This means that images from the same patient (despite being a different view) and therefore the same acquisition could never appear in the training and testing set at the same time. Here we present results for the merged short axis view labels by using features extracted exclusively either from image-based features or from DICOM slice orientation information.

Image features only DICOM orientation No jitter

With jitter features only Accuracy 87.1 ± 1.0% 90.2 ± 1.6% 98.8 ± 5.4%

Increasing the number of trees from 10 to 160 improves the prediction accuracy, however at the expense of increased training and classification computational cost . Beyond 80 trees the forest reaches a plateau with no further benefit (See Fig. 4(a)).

Conclusion

In this paper we presented an automatic cardiac view recognition technique from magnetic resonance images based on image appearance and DICOM plane orientation information.

The DICOM orientation based features can be used for our population with confidence when this information is present. When this is not the case, the pure image information method can be used instead. It is achieving precision on par or better with pure image based results of [START_REF] Zhou | Automatic view classification for cardiac MRI[END_REF] while being much simpler and not requiring to train any landmark detectors.

By jittering the dataset with extra translations, rotations and scales we observed noticeable benefits to view recognition. Extra jittering steps might be useful to improve robustness to intensity differences between different acquisitions. If enough examples are collected, this method can be simply extended to other pathology specific views such as in congenital heart diseases and acquisition sequences different from SSFP.

We plan to use our method to organize collections of images, semantic image retrieval and parsing of medical literature.

Figure 1 :

 1 Figure 1: Examples of the main left ventricular long and short axis cardiac MR views.

  1(a) to 1(c)) are used to visualize different regions of the left atrium, mitral valve apparatus, and left ventricle. The 3-chamber and left ventricular outflow tract (fig. 1(d)) views provide visualization of the aortic root from two orthogonal planes. The 4-chamber view enables visualization of the tricuspid valve and right atrium.

  Discriminative pixels from image miniatures are chosen from a random pool as features for a classification forest. We jitter the training dataset to improve robustness to differences in the acquisitions without the need for extra data. (b) DICOM plane normals for different cardiac views.128x128 (4x4 tiles) 40x40 20x20 32x32 32x32 miniatures min & max pool Image tiles (c) Image feature maps as downsampled images and tile intensity minima and maxima.

Figure 2 :

 2 Figure 2: Our view recognition pipeline (above) and used features (below).

Figure 3 :

 3 Figure 3: Image miniature averages for each cardiac view after cropping the central square. Short axis view means are displayed from the apex SAXAp to above base slices SAXSB. Distinct patterns of the bright blood pool or dark lungs can be observed for each view.

  The effect of the number of trees.

  Confusion matrix from one fold using image features only.

Figure 4 :

 4 Figure 4: We selected the optimal parameters via cross validation. Despite the sample weighting the LVOT view recognition suffered from being underrepresented.

Acknowledgments. This work was partly supported by Microsoft Research through its PhD Scholarship Programme and ERC grant MedYMA. We used data and infrastructure made available through the Cardiac Atlas Project (www.cardiacatlas.org) [2][3]. This work uses scikit-learn toolkit [7].