D Lopes

S Clain
email: clain@math.uminho.pt

R M S Pereira

G J Machado

G Smirnov

I Vasilevskiy

Numerical simulation of breast reduction with a new knitting condition

Keywords: Breast surgery, numerical simulation, knitting conditions

Breast reduction is one of the most common procedures in breast surgery. The aim of this work is to develop a computational model allowing one to forecast the final breast geometry according to the incision marking parameters. This model can be used in surgery simulators that provide preoperative planning and training, allowing the study of the errors origin in breast reduction. From the mathematical point of view this is a problem of calculus of variations with unusual boundary conditions, known as knitting conditions. The breast tissue is considered as a hyperelastic material, discretized with three-dimensional finite elements for the body whereas the skin is modelled with two-dimensional finite elements on the curved surface. Although the model is of low precision, we show that it is sufficient for a satisfactory analysis of errors frequently done during breast reduction surgery and allows to understand how to avoid or correct them.

INTRODUCTION

Breast reduction surgery is a procedure used to remove excess fat, tissue, and skin from the breasts in order to achieve a breast size in proportion with the body and to alleviate the discomfort associated with overly weight breast. The surgeon first makes one or more cuts and then the excess tissue and skin are removed while the nipple and areola are repositioned. At last, the skin is closed with stitches to initiate the natural suturing process. To help surgeons to take a better informed decision, medical numerical simulation enables to carry out virtual surgery operations to visualize the post-operated breast and to analyse several scenarios. The two main areas of the breast modelling applications are aesthetic surgery [START_REF] Balaniuk | Cosmetic Breast Surgery Simulation[END_REF][START_REF] Cardoso | Computer Simulation of Breast Reduction Surgery[END_REF] and medical imaging analysis [START_REF] Azar | Methods for modeling and predicting mechanical deformations of the breast under external perturbations[END_REF].

Although the number of works on the breast soft tissue modelling has increased significantly during the two last decades [START_REF] Cardoso | Computer Simulation of Breast Reduction Surgery[END_REF][START_REF] Azar | Methods for modeling and predicting mechanical deformations of the breast under external perturbations[END_REF][START_REF] Ayache | Surgery Simulation and Soft Tissue Modeling[END_REF][START_REF] Perez-Del-Palomar | A finite element model to accurately predict real deformations of the breast[END_REF][START_REF] Delingette | Towards realistic soft tissue modeling in medical simulation[END_REF], the development of an adequate breast model still continues to be a challenging problem. The use of nonlinear models to simulate the skin, the muscles, and the tissues is well-developed in the biomechanics context [START_REF] Hernádez-Gascón | Computational framework to model and design surgical meshes for hernia repair[END_REF][START_REF] Fung | Biomechanics -Mechanical Properties of Living Tissues[END_REF][START_REF] Rajagopal | Development of a Three-Dimensional Finite Element Model of Breast Mechanics[END_REF] while the finite element method is a popular technique for the discretization [START_REF] Payan | Soft Tissue Biomechanical Modeling for Computer Assisted Surgery[END_REF][START_REF] Zhong | A real time finite element based tissue simulation method incorporating nonlinear elastic behavior[END_REF][START_REF] Lee | Biomechanical study on the edge shapes for penetrating keratoplasty[END_REF] using alternatively the weak formulation or the minimization framework [START_REF] Ball | Energy-Minimizing Configurations in Nonlinear Elasticity[END_REF][START_REF] Ciarlet | Three-Dimensional Elasticity[END_REF][START_REF] Picinbono | Non-linear anisotropic elasticity for real-time surgery simulation[END_REF]. The numerical simulation of this type of surgery gives rise to a very specific problem where one has to prescribe suturing conditions in order to provide realistic breast shape. Indeed, the junction of the two incision interfaces is complex: the interfaces are placed in contact while the border (skin) is sewn. From a numerical point of view, a technical difficulty derives from the non-matching of the left and right surface meshes which requires the development of a new numerical suturing algorithm. In [START_REF] Clain | Soft tissue modelling for analysis of errors in breast reduction surgery[END_REF] it was proposed a mathematical model and a numerical scheme to simulate the suturing after two typical incisions in the breast reduction context: a symmetric frontal cut and a rear cut. In this work we generalize the model in order to consider an asymmetric frontal cut.

From the mathematical point of view, the problem is formulated as a problem of calculus variations with unusual knitting conditions. The breast tissue is considered as a hyperelastic material and, although most of soft tissues are incompressible, we consider the breast as a compressible Neo-Hookean material as proposed in [START_REF] Cardoso | Computer Simulation of Breast Reduction Surgery[END_REF][START_REF] Cardoso | On the importance of Chassaignac's space in breast modelling[END_REF]. The breast tissue is modelled using three-dimensional elements while the skin is modelled using two-dimensional elements since its elastic properties are very different from those of the breast gland and fat tissue. Moreover, we include in the model the Chassaignac space since it is of great importance to simulate the breast submitted to the gravitational force.

The paper is organized as follows. We present in section 2 the mathematical formulation of the the suturing problem while in Section 3 we detail the numerical algorithms. The breast model in the gravitational field is given in Section 4 while its respective numerical scheme is designed in Section 5. Numerical tests are given in the sixth section where we present three typical configurations and in last section we present the conclusions of this work.

MODELLING

Suturing process model is based on two main ingredients: a nonlinear elasticity model of the breast and suturing conditions to perform the junction of the two interfaces. The last point is far from trivial since, in practice, the surgeon sutures the border of the interfaces between themselves while the inner parts of the interface only are in contact. Hence the boundary conditions should carefully be prescribed to correctly model the surgery.

Geometry

In order to mathematically model the suturing problem, we denote by Ω cut the domain which corresponds to the incised breast and prescribe the conditions to perform the suturing. Since the suturing process is very complex from a geometrical point of view, one has to carefully define the interfaces of the breast and its parametrization. In short, the breast is considered as a spherical cap where we perform two cuts. We plot in Figure 1 three technical views with the parameters given in Table I. We use the following convention: Γ represents a bi-dimensional surface while γ corresponds to a generic one-dimensional curve. To describe the suturing problem, we introduce the following set of boundaries we display in Figure 2 while Table I gives an overview of the notations. The breast is partially covered with a skin characterized by the boundary surface Γ F . The frontal

L = Γ F ∩ Γ L frontal left suturing line γ R = Γ F ∩ Γ R frontal right suturing line γ BL = Γ F ∩ Γ BL backward left suturing line γ BR = Γ F ∩ Γ BR backward right suturing line θ
angle of the suturing plane α L , α R left and right angles of the frontal incision incision is performed between the two interfaces denoted by Γ L and Γ R respectively sharing the line γ 0 = Γ L ∩ Γ R passing at point A, orthogonal to the trunk body. A second rear incision is performed between the breast and the trunk such that a portion of the breast remains in contact with the body via the interface Γ B while Γ BL and Γ BR are the left and right surfaces respectively deriving from rear incisions. Boundary lines γ L and γ R in red in Figure 2 correspond to the suturing line that the surgeon will join while γ BL and γ BR are the rear incision lines. Notice that one has

∂Γ α = γ α ∪ γ Bα ∪ γ 0 , α ∈ {L, R}.

The functional

To describe the breast, we use a nonlinear hyperelastic neo-Hookean compressible model as in [START_REF] Cardoso | On the importance of Chassaignac's space in breast modelling[END_REF] where we neglect the gravity force assuming that the suturing process is independent of the external field. The problem turns to be the minimization of the associated energy functional [START_REF] Ball | Energy-Minimizing Configurations in Nonlinear Elasticity[END_REF][START_REF] Ciarlet | Three-Dimensional Elasticity[END_REF][START_REF] Picinbono | Non-linear anisotropic elasticity for real-time surgery simulation[END_REF] where we determine the new position vector field f : p → f (p) which minimizes

J cut (f ) = Ωcut W br (∇f (p))dp + Γ F W sk (∇ || f (p))dS p (1)
subject to the constraints

f (p) = p, in Γ B , (2) b(f (p)) = 0, in Γ BL ∪ Γ BR , (3)
Π θ (f (p)) = 0 in Γ R ∪ Γ L , (4) f (p) -f (g(p)) = 0, in γ R , (5)
where Π θ is the equation of the plane where the suture of the left and right cuts will take place. Indeed, we suppose that for a small asymmetry, i.e. α L ≈ α R , the suturing interface may be approximated by the medium plane. Of course, for large asymmetries, this assumption is no longer valid and the problem of the location of the final interface turns to be more complicated.

The expressions for the volume strain-energy and the skin strain-energy densities are respectively represented by W br and W sk and given by

W br (F ) = µ br 2 tr(F F t ) -3 -2 ln (det(F )) + λ br 2 (det(F ) -1) 2 ,
where F = ∇f is the Jacobi matrix of f and (λ br , µ br ) are the Lamé parameters for the breast and

W sk (F || ) = µ sk 2 Ä tr(F || F t || ) -2 -2 ln det(F || ) ä + λ sk 2 det(F || ) -1 2 ,
where F || = ∇ || f is the Jacobi matrix of the superficial (skin) displacement f || and (λ sk , µ sk ) are the Lamé parameters for the skin. Notice that f (Ω cut ) will be the new breast configuration after the suturing that one has to determine. Condition [START_REF] Cardoso | Computer Simulation of Breast Reduction Surgery[END_REF] guarantees that the back interface remains in contact with the trunk while condition (3) states that the rear interfaces belongs to the trunk plane to model the suturing of the breast with the body. To achieve the suturing process with the frontal cut, two boundary conditions are prescribed. The first one (4) states that the junction occurs in a plane Π θ characterized by the angle θ we shall prescribe as the mean value of the left and right angles while the second condition (5) introduces the mapping g : γ R → γ L associated to the suturing. It provides the point-to-point correspondence between the left and right side of the cut, that will be given in section 6 for our specific application.

DISCRETIZATION OF THE SUTURING PROBLEM

We introduce a discrete version of problem (1)-( 5) based on a mesh with tetrahedrons of the breast and a piecewise continuous representation of the new position f .

Mesh and notations

We denote by T a mesh of Ω cut constituted of I non-overlapping tetrahedron cells τ i , i = 1, . . . , I, and N vertices P n = (P nx , P ny , P nz ) ∈ R 3 , n = 1, . . . , N . Moreover, T k , k = 1, . . . , K, represents the faces of the tetrahedrons of the mesh that belong to Γ F . Quantities |τ i | and |T k | represent the volume and the area of the cell and the triangle respectively. We also use a local indexation. We denote by P ij = (P ijx , P ijy , P ijz ) ∈ R 3 , j = 1, 2, 3, 4, the vertices of τ i and by P kj = (P kjx , P kjy , P kjz ) ∈ R 3 , j = 1, 2, 3, the vertices of T k . To discretize function f , we associate to each node P n an approximation f n ≈ f (P n ) and denote by f h the usual continuous piecewise linear function while vector Φ h = (f nx , f ny , f nz ) p collects the 3N components of the new positions. Vector Φ h corresponds to the new configuration but the entries are not the unknowns of the problem since some of them are characterized by the boundary or the suturing conditions. Therefore, we extract a sub-vector X h of Φ h which only contains the unknown values we shall use in the minimisation process while the boundary conditions define an operator

X h → Φ h = B(X h )
which provides the other entries to complete vector Φ h .

Boundary and suturing conditions

Condition (2) yields that f n = P n for any P n on Γ B while relation [START_REF] Azar | Methods for modeling and predicting mechanical deformations of the breast under external perturbations[END_REF] implies that for any P n on Γ BL ∪ Γ BR we set f ky = 0 to enforce the suturing of the rear cut on the trunk plane y = 0. To apply Constraint (4), for any P n ∈ Γ R ∪ Γ L , we use f ny and f nz as unknowns, i.e. f ny and f nz are entries of X h . Therefore, we compute the last component with the relation

f nx = tan(θ)f nz .
Constraint [START_REF] Perez-Del-Palomar | A finite element model to accurately predict real deformations of the breast[END_REF] is the main condition since it describes the surgical procedure. For a node P n ∈ γ L , f ny and f nz are entries of X h and we compute the last component with f nx = tan(θ)f nz thanks to condition [START_REF] Ayache | Surgery Simulation and Soft Tissue Modeling[END_REF]. Assume now that function g : γ R → γ L is given (defined in section 6), we proceed in the following way at the discrete level. For a node P n ∈ γ R , we determine the edge e = [P α P β ] ∈ γ L such that g(P n ) ∈ e. We then define the unique barycentric coordinates such that λ nα P α + λ nβ P β = g(P n ) (see Figure 3).

P α P β P n f n f β f α g(P n ) Figure 3
. Discretization of function g. For any point Pn ∈ γ R , we determine an approximation g(Pn) ∈ γ L between two nodes Pα, P β of an edge belonging to γ L . Barycentric coordinates then derive from the relative position of g(Pn).

To preserve the relative position of f n with respect to f α and f β in the final configuration, we set

f n = λ nα f α + λ nβ f β .
Note that the barycentric coordinates will be computed once as a preprocessing step and none of the coordinates of P n is an entry of X h . For any P n , the choice of the associated edge e and the barycentric coordinate is equivalent, at the discrete level, to completely define function g.

Discretization of the functional

We now introduce a discrete version of the energy functional setting

J h (X h ) = J h (B(X h )) = J h (Φ h ) = J 1 h + J 2 h
constituted of the volume energy J 1 h and the surface energy J 2 h we detail hereafter.

Triple integrals of the functional

For a new configuration characterized by the approximation p ∈ Ω cut → f h (p) ∈ R 3 and stored in vector Φ h , the internal energy on tetrahedron τ i is given by

W τi = |τ i | Å µ 2 tr(F i F t i ) -3 -2 ln(det(F i )) + λ 2 det(F i ) -1 2 ã ,
where F i is the 3 × 3 matrix solution of the linear system

f i2 -f i1 = F i (P i2 -P i1 ), f i3 -f i1 = F i (P i3 -P i1 ), f i4 -f i1 = F i (P i4 -P i1 ).
The total internal energy is then approximated by

J h 1 = τi∈Ωcut W τi . 3
(O) = O, f || (A) = A , f || (B) = B . The Jacobian matrix of f || is the constant matrix Jf || = A = ï a b c d ò .
To determine the matrix, one writes

ï a b c d ò ï OB 0 ò = ï O B 0 ò , ï a b c d ò ï OA cos(α) OA sin(α) ò = ï O A cos(α ) O A sin(α ) ò ,
where α = ∠(OA, OB) and α = ∠(O A , O B ). The first linear system gives c = 0 and a = O B OB . Substituting these expressions in the second linear system we obtain

d = O A sin(α ) OA sin(α) , c = O A cos(α ) -a OA cos(α) OA sin(α) .
The superficial energy on triangle T for the skin is then given by

W T = |T | Å µ sk 2 tr(Jf || Jf t || ) -3 -2 ln(det(Jf || )) + λ sk 2 det(Jf || ) -1 2 ã
and the whole superficial energy is approximated by

J 2 h = T ∈Γ F W T .

Minimization algorithm

For a given vector X h and the boundary conditions, we deduce vector Φ h = B(X h ), hence the continuous linear piecewise function f h we use to compute the discretize energy functional. In order to determine the minimizer Xh of the discrete functional J h (X h ), we employ the conjugate gradients method. To this end, one has to compute an approximation of the derivative of the discrete functional in directions = 1, . . . , #X h . Since the function is nonlinear and has complex structure, numerical derivatives are computed in a very simple way. For example, an approximation of derivative with respect to direction at point X is given by

∂ J h (X) ≈ J h (X + ) -J h (X) ,
where is a vector of zeros except value for the -th entry.

BREAST UNDER GRAVITATIONAL FIELD

After the suturing, the two faces Γ L and Γ R are in contact and the question arises of the breast shape in the gravitational field in function of the trunk position. Two new contributions are now considered for the energy function: the gravity force and the introduction of the Chassaignac space. The last item corresponds to a mobile zone located between the breast and the trunk which acts as a spring to maintain the breast close to the trunk. The initial geometry Ω gra we shall deal with is the breast after the suturing process where the two interfaces Γ L and Γ R are in contact. To determine the new shape, we assume that surfaces are stick together in such a way that the relative position of the nodes of Γ R (resp. γ R ) with respect to the nodes of Γ L (resp. γ L ) will be preserved when the gravity field is presented. Notice that no new meshing have been performed hence the two areas are still considered as distinct parts of the boundary with a non-matching mesh. Therefore a specific condition, named the stick condition, is added to maintain the two faces in contact. Figure 4 displays the breast geometry domain Ω grav deriving from the suturing process and we introduce two new ingredients for the modelling. A new curve γ D on the infra-mammary fold is defined where the breast is fixed. We also merge the interfaces Γ B , Γ BL and Γ BR into a unique interface Γ C where the Chassaignac condition is prescribed. The problem turns to the minimization of the associated energy functional where we determine the new position vector field f : p → f (p) which minimizes

J gra (f ) = Ωgra W br (∇f (p))dp + Γ F W sk (∇ || f (p))dS p - Ωgra ρ a g • f (p)dp + Γ C c f (p) -p dS p (6)
subject to the constraints

f (p) = p, in γ D , ( 7
) h(f (p)) = 0, in Γ C , ( 8
) f (p) -f (G(p)) = 0, in Γ R ∪ γ R , (9)
where a g is the gravity acceleration, c the Chassaignac coefficient and ρ the density. Condition [START_REF] Hernádez-Gascón | Computational framework to model and design surgical meshes for hernia repair[END_REF] corresponds to fix the breast displacement on curve γ D while condition (8) maintains the breast on the trunk. Condition ( 9) is the stick condition that preserves the junction between the left and right interfaces previously joined in the suturing process. Function G results from the suturing process which defined the mapping from Γ L (resp. γ L ) onto Γ R (resp. γ R ) and will be given in section 5.1.

DISCRETIZATION OF THE BREAST UNDER THE GRAVITATIONAL FIELD

We introduce a discrete version of problem ( 6)-( 9) based on a mesh of tetrahedrons of the breast and a continuous piecewise representation of the displacement f . Due to the new boundary and constraint conditions, we have to deal with a new subset of unknowns Y h from Φ h and a new operator C such that Y h → Φ h = C(Y h ) models all the constraints (7)-(9).

Boundary, Chassaignac and sticking conditions

Condition [START_REF] Hernádez-Gascón | Computational framework to model and design surgical meshes for hernia repair[END_REF] yields that f n = P n for any P n ∈ γ D while relation [START_REF] Fung | Biomechanics -Mechanical Properties of Living Tissues[END_REF] implies that for any P n ∈ Γ C , we set f ny = 0 to maintain interface Γ C on the trunk plane y = 0. Constraint [START_REF] Rajagopal | Development of a Three-Dimensional Finite Element Model of Breast Mechanics[END_REF] means that the two interfaces previously joined and stick will preserve their relative position and we define function G in the following way. At the initial stage, after reading the mesh, we compute the barycentric coordinates for points P n on γ R and Γ R corresponding to the suturing process:

• If P n ∈ γ R , we determine the edge e = [P α , P β ] ⊂ γ L such that P n also belongs to e since the two curves are joined. We then define the unique barycentric coordinates such that G(P n ) = λ na P α + λ nβ P β = P n . Hence the constraint on the new location f n associated to point P n writes f n = λ na f α + λ nβ f β where f α and f β are the new positions of P α and P β , respectively. Notice that f n does not belongs to the entries of Y h since we deduce it from the linear combination of P α , P β . • If P n ∈ Γ R , we determine the triangle f = [P i P j P k ] ⊂ Γ L such that P n also belongs to f since the two surfaces are sutured. We then define the unique barycentric coordinates such that G(P n ) = λ ni P i + λ nj P j + λ nk P k = P n . Hence the constraint on the new location f n associated to point P n writes f n = λ ni f i + λ nj f j + λ nk f k where f i , f j and f k are the new positions of P i , P j and P k , respectively. Notice once that f n does not belongs to the entries of Y h since we deduce it from the linear combination of P i , P j , P k .

Discretisation of the functional

We introduce a discrete version of the energy functional setting

J h (Y h ) = J h (C(Y h )) = J h (Φ h ) = J 1 h + J 2 h + J 3 h + J 4 h with J 3 h = - τi |τ i | 4 ρa g (f (P i1 ) + f (P i2 ) + f (P i3 ) + f (P i4 ))
and

J 4 h = Tj ⊂Γ C |T j | 3 c ( f (P i1 ) -P i1 + f (P i2 ) -P i2 + f (P i3 -)P i3 )) .
Minimization of functional J h is performed in the similar way as presented in section 3.4.

NUMERICAL SIMULATIONS

The original breast is modelled as a truncated sphere characterized by the its principal radius and the length of the truncated portion. Figure 5 displays the initial configuration of the breast, i.e. before the incisions. We would like to highlight that more complex geometries can easily be considered and the suturing method is not limited to that particular configuration.

The suturing mapping g

Function g is the fundamental information to describe the suturing process by defining the mapping from γ R onto γ L which physically represents the skin junction performed by the surgeon. To this end, we propose a method to provide such a mapping based on the principle of equivalent relative distance. More precisely, since A is a fixed point, for any P ∈ γ R , we compute the distance of arc ÃP denoted by | ÃP | and define the relative distance

r(P ) = | ÃP | |γ R |
where |γ R | is the length of the curve γ R . We then take the associated g(P ) such that | Ȧg(P )| = r(P )|γ L |. At the numerical level, for a node P n ∈ γ R , we seek the edge e = [v 1 , v 2 ] ⊂ γ L such that g(P n ) ∈ e and compute the barycentric coordinates by solving equation

g(P n ) = λ 1 v 1 + λ 2 v 2 , λ 1 + λ 2 = 1.
We store in a vector the correspondence between vertices and edges and the associated barycentric coefficients.

Symmetric suturing

In the first numerical test, we performed a symmetric frontal incision with α L = α R = 30 • with no rear cut. Figure 6 gives two views of the configuration before the suturing while we display in Figure 7 the final configuration after the suturing. We have chosen θ = 0 • to preserve the symmetry of the junction. We notice that the symmetry of the shape is preserved and we have check the convergence to a limit configuration using finer meshes.

We then evaluate the breast shape subject to the gravity field for a person in stand up position. Figure 8 displays the frontal and side view of the breast under the gravity field and we notice that the symmetry is also preserved and the junction obtained is a straight line. Computations have been carried out with finer meshes to check that the solution converges to a limit one. The presented solution involves a mesh of 1246 tetrahedra and is very closed to the one obtain with finer meshes (less that 1% of relative errors).

Asymmetric suturing

In some situations, the surgeon has to perform asymmetric incisions and we experiment other code to assess the its capacity to handle such a more complex case. Figure 9 presents the initial breast Computations are carried out with a 1376 tetrahedra mesh taking the angle θ 0 = (-15 + 45)/2 of the suturing plan as the (algebraic) mean values of the two angles. Figure 10 plots the shape of the breast after the suturing process and the asymmetry is noticeable. As in the former case, we obtained shape convergence using finer meshes which guarantee the stability of the method.

We perform the computation of the suturing breast under the gravity field and we report that the suture location is no longer a straight line as displayed in Figure 11. Convergence to a limit configuration has been checked using finer meshes.

Frontal and rear suturing

As a final test, we consider a more complex case where we perform both a front and a rear incision of the breast. Figure 12 shows the mesh corresponding to two incisions while Figure 13 prints the configuration after the numerical suturing process. We use a 1539 tetrahedra mesh presented in Figure 13 but have experimented several finer meshes to check the convergence. No significance differences (lower than 1%) have been reported with respect to the solution. also present in Figure 14 the picture corresponding to the deformation of the breast when applying the gravity field. Such a final snapshot will represent the expected shape of the breast after the operation.

CONCLUSION

The software developed for breast reduction modelling allows one to forecast the final breast geometry according to the incision marking parameters. The main point is to perform symmetric or non-symmetric suturing between two interfaces with non-matching meshes while preserving the final junction when we apply the gravitational field. Accuracy of the approximation with rough meshes is sufficient for a satisfactory analysis of errors frequently done during breast reduction surgery and allows to understand how to avoid or correct them. The next step will be the comparison of the simulations with real surgery cases to validate our approach from a practical point of view.

Figure 1 .

 1 Figure 1. Scheme of breast reduction surgery: the middle panel is the frontal view while the left and right panels report the parameters for the left and right parts of the breast.

γBRFigure 2 .

 2 Figure 2. Full breast (left) and half breast (right) in a 3D frontal view (top row) and rear view (bottom row).

Figure 4 .

 4 Figure 4. 3D rear views for the gravity problem: full breast (left) and half breast (right).

Figure 5 .

 5 Figure 5. Mesh of the pre-operative breast: frontal view (left) and lateral view (right).

Figure 6 .

 6 Figure 6. Mesh of the breast after a symmetric incision: frontal view (left) and lateral view (right).

Figure 7 .

 7 Figure 7. Shape of the breast after the suturing without gravity field: frontal view (left) and lateral view (right).

Figure 8 .

 8 Figure 8. The breast under the gravity field: frontal view (left) and lateral view (right).

Figure 9 .

 9 Figure 9. Mesh of the breast after an asymmetric incision: frontal view (left) and lateral view (right)..

Figure 10 .

 10 Figure 10. Mesh of the breast after suturing: frontal view (left) and lateral view (right).

Figure 11 .

 11 Figure 11. The breast under the gravity field: frontal view (left) and lateral view (right).

Figure 12 .

 12 Figure 12. Mesh of the breast after a symmetric front incision and a rear incision: frontal view (left) and lateral view (right).

Figure 13 .

 13 Figure 13. Breast shape after the suturing process: frontal view (left) and lateral view (right).

Figure 14 .

 14 Figure 14. The suturing breast in the gravity field: frontal view (left) and lateral view (right).

Table I .

 I List of symbols and notations.

	Γ B	back side of the breast that is fixed to the chest
	Γ F	surface associated to the skin of the breast
	Γ L , Γ R	left and right cuts performed during the frontal incision
	Γ BL , Γ BR	left and right cuts performed during the rear incision
	A	point on the front of the breast where the first incision starts
	γ 0 = Γ L ∩ Γ R	line between point A and the chest
	γ

  .3.2. Surface integrals of the functional The discrete piecewise linear function f h transforms a triangle T k with vertices OAB into a triangle T k with vertices O A B . Since the translation and the rotation do not change the stress due to the deformation, we assume that O B is collinear to OB and A belongs to the same plane than triangle OAB. Function f || is a two-dimensional function locally given by f ||

ACKNOWLEDGEMENTS The authors are grateful to Professor Augusta Cardoso and Professor Horácio Costa. This work was supported with the Portuguese Funds from the "Fundac ¸ão para a Ciência e a Tecnologia" through the Project EXPL/MAT-NAN/0606/2013. The authors were also supported by the Research Centre of Mathematics and the Research Centre of Physics of the University of Minho with the Portuguese Funds from the "Fundac ¸ão para a Ciência e a Tecnologia", through the Project PEstOE/MAT/UI0013/2014 and PEst-C/FIS/UI0607/2013.


domain: math spi

text:
A Erdozain
email: aralar.erdozain@inria.fr

V Péron

Impedance Transmission Conditions for the Electric Potential across a Highly Conductive Casing

Keywords: Impedance Conditions, Electric Potential, Borehole, Casing, Resistivity

We present Impedance Transmission Conditions (ITCs) for the electric potential in the framework of borehole through-casing resistivity measurements. Such ITCs substitute the part of the domain corresponding to a highly conductive casing. The naturally small thickness of the casing makes it ideal for exhibiting ITCs. We numerically observe the delivered order of accuracy.

Introduction

Borehole resistivity measurements are commonly used when trying to obtain a better characterization of the earth's subsurface. Often a metallic casing is employed to surround the well, which allows to protect the well and avoid possible collapses. The use of such casing highly complicates the analysis due to large contrast between the conductivities of the casing and the rock formations.

This work is motivated by realistic configurations [START_REF] Pardo | Sensitivity study of borehole-to-surface and crosswell electromagnetic measurements acquired with energized steel casing to water displacement in hydrocarbonbearing layers[END_REF] where the conductivity of the casing is σ c ≈ ε -3 when ε denotes the thickness of the casing. In this framework, our aim is to derive ITCs for the electromagnetic field across such a casing. As a first approach we derive ITCs for the electric potential.

We refer to [START_REF] Schmidt | Robust transmission conditions of high order for thin conducting sheets in two dimensions[END_REF] where the authors derive ITCs for eddy current models with a conductivity parameter of a thin sheet of the form σ ≈ ε -2 .

We first introduce the mathematical model. Then, we explicit two asymptotic models of order two and four. Finally, we numerically analyse the performance and order of accuracy of the ITCs.

The Mathematical Model

We consider a transmission problem for the static case of the electric potential, the governing equations read as follows The considered domain Ω ⊂ R 2 is formed by three subdomains Ω ε i , Ω ε c , Ω ε e , which are characterized by different conductivities σ i , σ e , σ c . These subdomains have rectangular shape and Ω ε c is a thin layer of uniform thickness ε, as shown in Figure 1. In (1), f i , f e , σ i , σ e and σ c are known data and u corresponds to the unknown, whose restrictions to the different subdomains are denoted as u|

                       σ i ∆u i = f i in Ω ε i σ e ∆u e = f e in Ω ε e σ c ∆u c = 0 in Ω ε c u i = u c on Γ ε i u c = u e on Γ ε e σ i ∂ n u i = σ c ∂ n u c on Γ ε i σ c ∂ n u c = σ e ∂ n u e on Γ ε e u = 0 on ∂Ω (1)
Ω ε i = u i , u| Ω ε e = u e , u| Ω ε c = u c .
In this framework, we address the issue of ITCs for u (as ε → 0) when the conductivity of the casing is σ c = αε -3 (α ∈ R).

Main Results

For developing ITCs, we perform a formal expansion of the solution u in power series of ε. This leads to a collection of problems which can be solved successively. Then truncating the series we build asymptotic models by obtaining ITCs between Γ ε i and Γ ε e .

Definition 1 Let u be the solution of problem (1). We say that u [k] satisfies an asymptotic model of order k + 1 when (for ε small enough)

||u -u [k] || L 2 (Ω ε i ∪Ω ε e ) ≤ Cε k+1 , C ∈ R.
We derive two asymptotic models of order two and four.

Order 2 model

σ i ∆u [1] i = f i in Ω ε i u [1] i = 0 on ∂Ω ε i σ e ∆u [1] e = f e in Ω ε e u [1] e = 0 on ∂Ω ε e Order 4 model                    σ i ∆u [3] i = f i in Ω ε i σ e ∆u [3] e = f e in Ω ε e u [3] = 0 σ∂ n u [3] = -α ε 2 ∆ Γ u [3] u [3] i = 0 on ∂Ω ε i \Γ ε i u [3] e = 0 on ∂Ω ε e \Γ ε e
Here, the jump and mean value of a function u across the domain Ω ε c are

[u] := u e | Γ ε e -u i | Γ ε i {u} := 1 2 u e | Γ ε e + u i | Γ ε i .
The validation of these models consists in proving estimates for u -u [k] (see Definition 1). Hereafter, we present numerical validations for each asymptotic model.

Numerical Results

We developed a finite element code to solve problem (1) and the two asymptotic models. We consider f i = 1, f e = 1 as right hand sides and we select σ i = 3, σ e = 5, σ c = ε -3 as the different conductivities.

We compute the L 2 error between the solution of problem [START_REF] Schmidt | Robust transmission conditions of high order for thin conducting sheets in two dimensions[END_REF] and the solution of each asymptotic model for different values of ε by using triangular elements and degree five Lagrange polynomials. We observe these results in Figure 2 and the corresponding slopes of the graphics in Table 1. Each numerical convergence rate converges to the formal order of accuracy. Table 1: Slopes corresponding to the curves of Figure 2 5 Perspectives Among the future perspectives, we would like to derive asymptotic models for different configurations including the non-static case of the electric potential and electromagnetic field. We plan to perform mathematical proofs to validate these asymptotic models.

Figure 1 :

 1 Figure 1: Domain of interest.

Figure 2 :

 2 Figure 2: L 2 error of the order 2 and order 4 model for different values of ε


A Shishlo

A Zhukov

Y C Chao

C Gong

F Jones

TRIUMF R Newhouse

P Chu

D Maxwell

T Pelaia Ii #

C K Allen

Y Zhang

R Fearn

L Fernandez

E Laface

ESS M Munoz

Sweden J Lund

Freed

P Gillette

P Laurent

G Normand

H Hale

CSNS Y Li

OPEN XAL STATUS REPORT 2015*

Keywords:

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Open XAL [ ] is an open source accelerator physics 1 software platform written in Java. The Open XAL project began in mid 2010 as a response to requests from the international accelerator physics community to adopt an open source accelerator physics platform based on XAL [ ] from the Spallation Neutron Source (SNS) at Oak 2 Ridge National Lab (ORNL) and establish a standard platform for accelerator physics software.

The goals of the project are to provide a common accelerator physics core to be developed in collaboration, provide a rapid development environment and to modernize the source code. Since the last status report [ ],

3 the development effort has shifted from porting code to verification, modernization and new functionality. The API has stabilized, and the project is production ready. This paper covers the collective contributions of the collaboration.

SOFTWARE HIERARCHY

The software is grouped among the core, extensions, plugins, services, applications and scripts. The core consists of packages that include the accelerator object graph, the online model, abstract channel interface and general supporting tools such as common math libraries, archiving, messaging, database access and concurrency. This core is intended to be shared in common across all Open XAL distributions, and it has no compile time dependency on other components.

Extensions and plugins are components which may optionally depend on other extensions and plugins and the core. A plugin differs from an extension in that the core depends on a plugin being implemented (but not the actual implementation) whereas the core has no dependency on an extension. Plugins consist of controls channel access support and database adaptors. For example, the core specifies generic database tools, but a database adaptor plugin provides support for a specific database driver. Extensions include packages of general use in applications, scripts and services but not in the core. Examples of extensions include the application framework, service framework, lattice generation, fitting, scanning and widgets such as plotting tools.

Applications, scripts and services are executables and may depend upon the core, plugins and extensions. Applications and scripts are launched by users and have a user interface whereas services are headless and run perpetually in the background independent of user interaction. Applications differ from scripts in that applications are written in Java and are thus fully compiled prior to runtime versus scripts that may be written in a Java variant of Ruby or Python scripting languages and are not compiled prior to runtime.

A top level site directory contains optional site specific resources and configurations that take precedence over resources and configurations of the same name in the _____________________________________________ * This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doepublic-access-plan). # pelaiata@ornl.gov main source tree and thus offer an opportunity for site specific customization.

COLLABORATION

The

Code Sharing

The code is placed in Git [ ] repositories [ ] hosted on 4 5 SourceForge and a ticket system is used to track tasks. The Git subtree feature is used extensively so subtrees of the code can be shared while allowing each site to choose which applications, scripts, extensions, plugins and services to include in their site specific branch.

SITE STATUSES

The operational status varies across the facilities and hence each facility has different needs for Open XAL, and the status at each site varies accordingly.

CSNS

MEBT commissioning for CSNS will use applications based on a mature XAL variant rather than Open XAL due to the confluence of commissioning with the nascent development of Open XAL. Several new applications have been developed to accommodate the Rapid Cycling Synchrotron (RCS), magnet database, lattice database, etc. These applications will be ported to Open XAL in the near future.

ESS

The Beam Physics group at ESS decided in November 2014 to use Open XAL for beam commissioning. The Open XAL model has been modified to handle field maps for cavities [ ] and was benchmarked against TraceWin 6 [ ]. A python integration environment based on JPype [ ] 7 8 has also been developed. The Open XAL solver extension was used for matching the initial parameters of the beam based on ESS specific criteria. Automated installation scripts were created for Mac and Linux operating systems. Optics files are generated as custom XML files called "LinacLego" and are convertible to standard Open XAL format using an extension.

FRIB

Open XAL is currently a backup option for high level accelerator physics applications at FRIB, and thus has an uncertain future here. Several Open XAL applications have been tested in the control room during the September 2014 commissioning of a new cryomodule successfully demonstrating use of Open XAL for linac commissioning and beam tuning. A model platform was developed which allows integration with IMPACT [ ] and 9 MADX [ ] as additional options for modeling. 10

GANIL (SPIRAL2 Project)

SPIRAL2 has adopted the OpenXAL September 2014 snapshot. Hibernate is used for database interaction, and the accelerator input XML files are generated from the database. The core has been modified slightly to support the concept of equipment in the accelerator object graph. Tools were developed to allow user interaction with process variables. Custom SPIRAL2 application adaptor and document subclasses have been developed as foundations for applications adding support for such things as displaying the accelerated beam label in the title bar, automatic combo-sequence selection, online/offline accelerator selection and the Nimbus look and feel. Several applications have been developed and others are being developed.

SNS

At SNS, the old XAL code has been frozen with the exception of critical bug fixes if necessary. All active applications (several dozen) have been ported to Open XAL, and it is now deployed as the default high level accelerator physics software platform. Several applications have been debugged and verified as production ready. XAL is still available as a fallback should critical bugs be encountered.

TRIUMF

At TRIUMF, a low energy empirical model has been successfully tested, and a magnet dithering and beam based alignment tool have been developed. The Experiment Automator application has been developed for automated data acquisition and contributed to the Open XAL community. An orbit correction package has been planned.

DEVELOPMENT HIGHLIGHTS

Several of the many developments since the last status report will be summarized here. Specifically, the ones discussed here are of common use across the sites.

Continuous Integration

ESS has setup a Jenkins continuous integration server for Open XAL providing analysis which has been used to improve the code.

Online Model

The online model has undergone major architectural improvements [ ] which we briefly highlight here. 11 Processing of simulation data has been separated from probe state population. This approach keeps analysis data separate from the trajectory and allows for richer simulation computations on the trajectory and its probe states.

Probes, probe states and trajectories have been modified to use Java Generics [ ] which allows for 12 compile time type checking and consistent typing across these classes.

The performance of fetching states from a trajectory has been improved significantly, and is especially noted when working with trajectories over a large number of states.

The lattice generator which maps model elements to device nodes has been rewritten making it simpler to maintain the code and support new device types.

The elapsed time mechanism has been refactored for RF Gaps leading to more sustainable code and removing the use of static variables to maintain state.

Site Customization

Site specific resources can be added and override existing resources of the same name if any. Site specific build properties can be specified and take precedence over the defaults.

Application Framework

The application framework now supports applications written entirely in a Java based scripting language such as JRuby [ ] or Jython [ ].

13 14

Services Framework

The

Tools

One old external library of open source widgets that had been included in XAL contained a very useful wheel switch which was the only control in use from that library. The code for the wheel switch was modernized, issues addressed and brought directly into the Open XAL source thus eliminating dependence on that external library.

Applications

Significant improvements have been made to the Virtual Accelerator application including, adding solenoid support, allowing set points to be edited directly and adding a data plot for beam size and positions. A new application that allows for experimental automation was developed and contributed to the master branch.

FUTURE PLANS

Here are a highlights of a few of the many tasks that lay ahead for Open XAL. Currently Open XAL depends on Java 7. Soon, we plan to embrace Java 8 as our new dependency so as to allow for new Java 8 constructs and packages in the code base.

Online Model

Plans are for each RF Gap to be able to maintain its own transit time factor functions. Also, the model needs to be benchmarked against other established machine codes.

JavaFX

With the adoption of Java 8, we can embrace JavaFX [ ] as a natural replacement for the Swing-based Bricks 18 [ ] GUI builder and framework. Support for JavaFX will 19 also be added to the application framework.

Unit Tests

Support for unit tests will be extended to other components beyond just the core.

  active collaboration consists of the China Spallation Neutron Source (CSNS) in Dongguan, China, European Spallation Source (ESS) in Lund, Sweden, Facility for Rare Ion Beams (FRIB) at Michigan State University in Lansing, MI, Grand Accélérateur National d'Ions Lourds (GANIL) in Caen, France, SNS in Oak Ridge, TN, and TRIUMF in Vancouver, Canada. Members from each facility contribute code, report issues and participate in workshops and monthly online meetings.

Proceedings of IPAC2015, Richmond, VA, USA

ACKNOWLEDGMENT

We are grateful for the leadership contributions of Sheng Peng of FRIB who organized online meetings and drove early momentum of the Open XAL collaboration.


Flavio Pressacco

Laura Ziani

Fibonacci numbers in a class of weighted majority games

Keywords: Fibonacci numbers, weighted majority games, homogeneous representation, minimal winning coalition, type weight vector, satellite games

come

Introduction and notations

Let Ω = {1, . . . , n} denote the set of non dummy players of a simple constant sum game in characteristic function form. A simple game is a mapping v : P(Ω) → {0, 1} such that v(∅) = 0 and v(Ω) = 1 and a coalition S ∈ P(Ω) is winning if its payoff v(S) = 1 and losing otherwise. A simple game is constant sum if v(S) + v(Ω/S) = 1 for any S. Moreover, S is minimal winning if, for any player i ∈ S, v(S/i) = 0.

A simple weighted majority game is described by the pair (w n ; q) where1 w n = {w 1 , . . . , w n } is the (weakly) ordered vector (w i ≤ w j for i < j) of individual weights and q the winning quota of the game. Thus w(S) = i∈S w i is the weight of S and w(S) ≥ q ⇔ v(S) = 1.

A simple constant sum homogeneous weighted majority game is described by its minimal homogeneous representation, that is the ordered vector h = (w n ; q) ∈ N n+1 with w 1 = 1, q = (1+w(Ω)) 2 and w(S) = q for any minimal winning coalition.

The vector w n of the minimal homogeneous representation induces a decomposition of Ω in equivalence classes K 1 , . . . , K j , . . . , K t . Each class groups all players of the same type, sharing the same individual weight and the corresponding (strongly) ordered type weight vector is w t = {w 1 , . . . , w j , . . . , w t } with w 1 = 1.

Coherently, h = (w t ; q) ∈ N t+1 , but to unequivocally describe a game it is necessary to introduce also the type cardinality vector k t = (k 1 , . . . , k j , . . . , k t ) with k j the number of players of type j in the game. Thus, a constant sum homogeneous weighted majority game is described by (k t ; w t ; q) ∈ N 2t+1 . Finally, we denote by s t = (s 1 , . . . , s j , . . . , s t ) with s j = |S ∧ K j | ≤ k j the "profile" of the coalition S, so that w(S) = t j=1 s j w j is an alternative formalization of the weight of S and w(Ω) = t j=1 k j w j = 2q -1.

Fibonacci games: main results

The Fibonacci sequence is defined by the well known finite difference equation:

f n = f n-1 + f n-2
holding for any n > 2 with initial conditions f 1 = f 2 = 1. Henceforth we exploit the "delayed" Fibonacci sequence g n = f n+1 for any n and denote, for any integer m, by g m ∈ N m the vector (g 1 , g 2 , . . . , g m ).

Definition 2.1. The set of Fibonacci games is the subset of constant sum homogeneous weighted majority games with h = g t+1 , that is (k t ; w; q) = (k t ; g t ; g t+1 ).

Remark 2.1. In any Fibonacci game with t types, t j=1 k j g j = 2g t+1 -1 constant for any k t : the winning quota is the same for all Fibonacci games with the same t.

Theorem 2.1. Both for t = 1 and t = 2 there is a unique Fibonacci game; their k t are respectively k 1 = (3), and k 2 = (3, 1); for any t > 2, a game is a Fibonacci game if and only if its k t = k(t, z) is given, for any value of the counter2 z = 1, . . . , ⌊(t -1)/2⌋ and with j 0 = t + 1 -2z, by:

k(t, z) = (k 1 , k ′ , k j 0 , k ′′ ) (2.1)
with k 1 = 2 + g t-1 -g j 0 , k j 0 = 2 and k ′ = 1 j 0 -2 , k ′′ = 1 2z-1 ,
or, putting z = ⌊(t + 1)/2⌋ in order to have a complete indexation of the games, by:

k(t, z) = k(t, ⌊(t + 1)/2⌋) = (k 1 , 1 t-1 ) (2.2)
with k 1 = 2 + g t-1 .
Corollary 2.1. For any t, there are exactly

Φ(t) = ⌊(t + 1)/2⌋ (2.3)
Fibonacci games with t types.

Let us denote by n(t, z) the total number of non dummy players in a Fibonacci game with t types and index z.

Result 2.1. For any t, z 1 < z 2 ⇔ n(t, z 1 ) < n(t, z 2 ).

Structures of k(t, z) vectors

According to Definition 2.1, given t there are no degrees of freedom in the choice of the sequence of type weights and of the winning quota. Hence, a Fibonacci game is unequivocally described by its k t vector which follows the rules described in Formulas (2.1) and (2.2). Table 1 We distinguish between two patterns of structures: for z = ⌊(t + 1)/2⌋, (the last structure of each row) all but the bottom3 type components are 1, while the bottom component plays a balancing role; for all the other z there is just one non bottom component equal to 2 in place j 0 = j(t, z) = t + 1 -2z, all the other non bottom components are still 1 and, again, the bottom component plays the balancing role. [START_REF] Rosenmüller | Weighted majority games and a the matrix of Homogeneity[END_REF] The behavior of n(t, z) The following properties hold for any feasible combination of (t, z):

Property 4.1. n(t, 1) = t + 2 (4.1) Property 4.2. n(t, z + 1) -n(t, z) = g t-2z (4.2) Property 4.3. n(t, 2) -n(t -1, ⌊t/2⌋) = 2 (4.3)
Here, we distinguish between two kinds of Fibonacci games: the first kind corresponds to z = 1 and Formula (4.1) implies that there is just one Fibonacci game (of this first kind) for any n > 2. Such games have been described by Isbell ([1], p. 185) while studying the class of coalitionally Parsimonious games (for details see also [START_REF] Pressacco | Bilateral symmetry and modified Pascal triangles in Parsimonious games[END_REF]).

For z > 1 (second kind), there is one Fibonacci game (of this second kind) for sparse values of n.

Summing up, there are values of n corresponding to two Fibonacci games (one of the first and one of the second kind), while for other values of n there is just one Fibonacci game of the first kind.

At first sight the values of n associated to the second kind of games seem to follow a chaotic rule. A regularity emerges once we consider the following table of the differences generated by Property 4.2 (see also Remark 6.1). Each row is a backward alternate sequence of ⌊(t -1)/2⌋) Fibonacci numbers ending with g 1 = 1 or g 2 = 2 for odd or, respectively, even values of t. All the columns are nothing but the whole sequence of g ∞ starting at the row t = 1 + 2z.

Proofs of Theorem 2.1

Proof. For t = 1:

q = g 2 = 2 w(Ω) = k 1 • g 1 = k 1 = 2q -1 = 3
Hence there are three players with weight g 1 = 1, k 1 = (3) and Φ(1) = 1.

For t = 2:

q = g 3 = 3 w(Ω) = 2 j=1 k j g j = 2q -1 = 5 Hence, k 1 • g 1 + k 2 • g 2 = k 1 • 1 + k 2 • 2 = 5.
There are two solutions k 2 = (3, 1) and k 2 = (1, 2) but the second is not feasible because the profile s 2 = (0, 2) identifies a minimal winning coalition S with w(S) = 4 > q, incompatible with the homogeneous character. On the other side, the first solution is feasible with profiles (3, 0) and (1, 1) of the minimal winning coalitions, so that Φ(2) = 1.

For t > 2, we show at first that the conditions (2.1) or (2.2) are necessary for Fibonacci games and then that they are also sufficient.

Proof of necessity

Property 5.1. In any Fibonacci game with t > 1, the coalition S whose profile s t = (0 t-2 , 1 2 ), i.e. made by a top and a last but top player, is minimal winning.

Proof.

w(S) = s t-1 g t-1 + s t g t = g t-1 + g t = g t+1 = q
Property 5.2. In any Fibonacci game with t > 1, k t = 1.

Proof. Otherwise the coalition S : s t = (0 t-1 , 2) would be minimal winning with w(S) = 2g t > g t + g t-1 = q a contradiction with the homogeneous character of the game.

Property 5.3. In any Fibonacci game with t > 2, k t-1 ≤ 2.

Proof. Otherwise the coalition S : s t = (0 t-2 , 3, 0) would be minimal winning as 2g t-1 < g t-1 + g t = q with w(S) = 3g t-1 > 2g t-1 + g t-2 = g t-1 + g t = q.

Property 5.4. In any Fibonacci game with t > 3, k t-2 = 1.

Proof. Suppose on the contrary, k t-2 > 1 and consider the coalition S : s t = (0 t-3 , 2, 0, 1). S would be minimal winning too as g t +g t-2 < g t +g t-1 = q with w(S) = g t +2g t-2 > g t +g t-2 +g t-3 = g t +g t-1 = q. A contradiction.

Properties 5.3 and 5.4 may be generalized to types (th) with h odd and respectively h even.

Property 5.5. In any Fibonacci game it is k t-h ≤ 2 for h odd < t -1.

Proof. Proof of Property 5.5 is based on the following lemma:

Lemma 5.1. In a Fibonacci game consider any player of type j 0 = (th) with h odd < (t -1); there exists a minimal winning coalition S in which this player is the weakest. The profile of S is given by: s t = (s 1 , . . . , s j , . . . , s t ) with s j = 0 for j < j 0 ; s j 0 = 1; s j = 1 for j > j 0 and j = (j 0 + 1) mod 2 and s j = 0 for j > j 0 and j = j 0 mod 2.

Proof. w(S) = g j 0 + j>j 0 and j=(j 0 +1) mod 2

g j = g t-1 + g t = g t+1 = q
and moreover 4 :

g t+2 - t j=1 g j = 2 (5.2)
Furthermore, Properties 5.6 and 5.7 imply that, for any j = 2, . . . , t, k t satisfies conditions (2.1) or (2.2). In the first case and putting j 0 = t + 1 -2z, we obtain:

k 1 = k 1 (t, z) = 2g t+1 -1 -( t j=2 g j + g j 0 ) = = g t+1 + (g t + g t-1 ) - t j=1 g j -g j 0 = = g t+2 + g t-1 - t j=1 g j -g j 0 = = 2 + g t-1 -g j 0 (5.3)
In particular, for z = 1, k 1 = 2. In the second case, we get immediately:

k 1 = k 1 (t, ⌊(t + 1)/2⌋) = 2 + g t-1 (5.4)
To resume, by Properties 5.6 and 5.7, and by Formulas (5.3) and (5.4), necessary conditions to be satisfied by the k t vector of a Fibonacci game with t > 2 are those described in Formulas (2.1) and (2.2).

Proof of sufficiency

The proof is based on an adaptation to our problem of the test for homogeneity of a weighted majority games developed as Basic Lemma in [START_REF] Rosenmüller | Weighted majority games and a the matrix of Homogeneity[END_REF] (Theorem 1.4) and [START_REF] Rosenmüller | Homogeneous games: recursive structure and Mathematics of Operations Research[END_REF] (pp. 312).

Henceforth we use the following definitions and notations:

• G 0 = (k 0 t ; g 0 t ; λ 0 = g t+1 ), "seed" game with k 0 t given by Formulas (2.1) or (2.2). In particular, G 0,z , seed game as a function of z = 1, . . . , ⌊(t + 1)/2⌋.

• G r = (k r t ; g r t ; λ r ), any game of a r generation 5 derived by the seed G 0 . Note that the dimension t of G r depends on r. Whenever useful, we utilize t(r). In particular, G r,z , game of a r generation derived by the seed G 0,z .

• Feasible (not feasible) game: a G r for which t j=2 k r j • g r j ≥ λ r (< λ r ).

• For any given type index j 0 > 1, an "intermediate" player in the game G r is a player of type j 0 , a "large" player is a player of any type j > j 0 , a "small" player is any player of type 1 < j < j 0 .

• No bottom coalition of a feasible game G r : a coalition S r whose profile s r has s r 1 = 0.

4 By induction: indeed it is immediate to check that (5.2) is true for t = 1 and, if true for t, is true also for t + 1.

5 Note that the winning quota λ r of G r may well be greater than t j=1 k r j g j r+1 2

, i.e. for r > 0, G r may not be constant sum.

• Dominant coalition6 of a feasible game G r a coalition S r such that:

1. S r is no bottom 2. there is an index j 0 (G r ) > 1 (or shortly j r 0 > 1) and a positive integer c r , such that s r = (s r j = 0 for j < j r 0 ; s r j = c r (1 ≤ c r ≤ k r j 0 ); s r j = k r j for j > j r 0 )

i.e. S r is made by all large players and by some (may be all) intermediate players, so that j r 0 is the type index corresponding to the weakest player in the dominant coalition S r . 3. w(S r ) = t(r) j=1 s r j g r j = λ r .

• Generation rules. Let G r a feasible game, S r its dominant coalition and j r 0 the corresponding type index. To any j r 0 ≤ j ≤ t(r), associate the game G r+1 (j) = (k r+1 , g r+1 , λ r+1 ) in which k r+1 = k r+1 (j):

for j > j r 0 and k r j 0c r > 0:

k r+1 j 0 = (k r 1 , k r 2 , . . . , k r j 0 -1 , k r j 0 -c r )
, the dimension of g r+1 is j 0 and λ r+1 = g j for j = j r 0 or j > j r 0 and k r j 0c r = 0:

k r+1 j 0 -1 = (k r 1 , k r 2 , . . . , k r j 0 -1
), the dimension of g r+1 is j 0 -1 and λ r+1 = g j Verbally the dominant coalition S r generates as many "satellite" games G r+1 (j) as the number of types j included in the coalition. The players of the satellite game associated to a type j are all players of type i < j not belonging to S r ; the winning quota is g j . The idea is that any satellite game looks for the dominant coalition of the satellite able to replace one player of type j in the generating game while preserving the homogeneous character of the "seed" game. In turn, the dominant coalition gives rise to other satellite games of next generation and so on. Now, our adapted version of Basic Lemma is:

Lemma 5.6. A seed weighted majority game G 0 is homogeneous if, in the set of all G r games obtained by the generation rules, do not exist feasible games lacking the dominant coalition.

To prove the Lemma 5.6 we show that any game of the set G r , (r = 0, 1, . . .) generated by a seed G 0 coherent with formulae (2.1) and (2.2) either is not feasible or admits a dominant coalition.

In the proof we will exploit the following relation7 concerning Fibonacci numbers:

g t-2 + 2g t-1 + h j=0 g t+j = g t+h+2 (5.5)
Proof. Let us distinguish three cases: Case α): z = ⌊(t + 1)/2⌋

Let us denote G 0,⌊(t+1)/2⌋ = G 0 = (k 0 t ; g 0 t ; g t+1 ) any seed game with k 0 t = (k 0 1 , 1 t-1 ). The profile of the dominant coalition S 0 is s 0 t = (0 t-2 , 1 2 ) with k 0 t-1c = 0. Hence the satellite games associated to type t -1 and t are G 1,a = G 1 (t -1) = (k 1 , 1 t-3 ; g t-2 ; g t-1 ) and G 1,b = G 1 (t) = (k 1 , 1 t-3 ; g t-2 ; g t ). By Formula (5.2), G 1,b is not feasible, while G 1,a is feasible for any t = t(0) > 4. Let shortly write G 1,a = G 1 ; it is immediate to check that G 1 mimics the structure of G 0 : only the dimension is different: it is t = t(0) for G 0 and t = t(1) = t(0) -2 for G 1 . The profile of the dominant coalition S 1 of G 1 is s 1 = (0 t(1)-2 , 1 2 ) = (0 t-4 , 1 2 ) and the procedure may be recursively repeated until the integer r for which also G 1,r is no more feasible. This is resumed by the following Property 5.8. Let G 0 be a seed game of type α with t = t(0) > 2. For any generation r = 0, 1, . . . , ⌊(t -3)/2⌋ there is a unique feasible G r = (k r t ; g r t ; λ r ). The k r t of such game is k 0 t truncated at dimension t(r) = t(0) -2r, which is coherently the dimension of g r t , while λ r = g t(r)+1 . Feasible G r games have dominant coalition S r whose s r t(0)-2r = (0 t(0)-2(r+1) , 1 2 ). Case β): z = 1 Let G 0,z = G 0,1 = (k 0 t ; g 0 t+1 ; g t+1 ) any seed game with k 0 t = (k 0 1 , 1 t-3 , 2, 1). The profile of the dominant coalition S 0 is still s 0 t = (0 t-2 , 1 2 ) with k 0 t-1 -c = 1. Hence, the satellite games associated to type t -1 and t are G 1,a = G 1 (t -1) = (k 1 , 1 t-3 ; g t-2 ; g t-1 ) and G 1,b = G 1 (t) = (k 1 , 1 t-2 ; g t-1 ; g t ). Both games share the properties of G 0 of case α), that is behave as seed games and generate for r = 2, . . . sequences G r,a and G r,b according to Property 5.8. Of course the dimension of G 1,a = t(1, a) = t(0)-2, and of

G 1,b = t(1, b) = t(0) -1. Case γ): 1 < z < ⌊(t -1)/2⌋
Let us write G 0,z = (k 0 1 , 1 t-1-2z , 2, 1 2z-1 ). The profile of the dominant coalition S 0,z is s 0,z t = (0 t-2 , 1 2 ) with k 0 t-1c = 0. Let us denote the two satellite games by G 1z,a = G 1z (t -1) = (k 1z,a t ; g 1z,a t ; λ 1z,a ) and

G 1z,b = G 1z (t) = (k 1z,b t ; g 1z,b t ; λ 1z,b ). It is k 1z,a t-2 = k 1z,b t-2 = (k 0 1 , 1 t-1-2z , 2, 1 2z-3
), λ 1z,a = g t-1 , λ 1z,b = g t . Note that G 1z,a mimics G 0,(z-1) ; hence recursively we go back to G 0,1 , i.e. the seed game of case β).

In turn, as by Formula (5.5), j s 1z,b j g j = g t , G 1z,b has dominant coalition with profile s 1z,b t-2 = (0 t-2z-1 , 1, 2, 1 2z-3 ).

Hence G 1z,b gives rise to 2z -1 satellite games G 2z,b (j) for j = t(0) -2z, . . . , t(0) -2. It is G 2z,b (j) = (k 0 1 , 1 t-2z-2 ; g t-2z-1 ; g j ). By Formula (5.2), for all j, except j * = t(0) -2z, the G 2z,b (j) are not feasible; the survival feasible G 2z,b (j * ) mimics the behavior of the G 0 seed of the case α).

Summing up, in all three cases all G r games either are not feasible, or if feasible admit a dominant coalition and the conditions of the Basic Lemma are verified.

Proofs of Properties of Section 4

By definition:

n(t, z) = k 1 (t, z) + t j=2 k j (t, z) (6.1)
Proof. Proof of Property 4.1

For z = ⌊(t + 1)/2⌋ it is t j=2 k j (t, z) = t -1, by Formula (2.2) and k 1 = 2 + g t-1 by Formula (5.4) so that: n(t, z) = t + 1 + g t-1 (6.2)

For z < ⌊(t + 1)/2⌋ it is t j=2 k j (t, z) = t, by Formula (2.1) and k 1 = 2 + g t-1g t+1-2z by Formula (5.3) so that: n(t, z) = t + 2 + g t-1g t+1-2z (6.3) and in particular n(t, 1) = t + 2

Proof. Proof of Property 4.2.

For (z + 1) < ⌊(t + 1)/2⌋ by Formula (6.3): n(t, z + 1)n(t, z) = (t + 2 + g t-1g t+1-2(z+1) ) -(t + 2 + g t-1g t+1-2z ) = g t-2z+1g t-2z-1 = g t-2z

(6.4)

For (z + 1) = ⌊(t + 1)/2⌋ by Formulas (6.2) and (6. n(t, z + 1)n(t, z) = (t + 1 + g t-1 ) -(t + 2 + g t-1g t+1-2z ) = g t+1-2z -1 (6.5) but z = ⌊(t -1)/2⌋ = (t/2) -1 for t even and (t -1)/2 for t odd so that:

t + 1 -2z = 3 for t even and t + 1 -2z = 2 for t odd which means also on one side t -2z = 2 for t even and t -2z = 1 for t odd and on the other side g t+1-2z = 3 for t even and g t+1-2z = 2 for t odd and finally g t+1-2z -1 = 2 for t even and g t+1-2z -1 = 1 for t odd so that in both cases still n(t, z + 1)n(t, z) = g t-2z (6.6) which completes the proof of Property 4.2. Note that this property implies the strict monotony of n(t, z) for any t and then, as a by product, gives a proof of Result 2.1.

Proof. Proof of Property 4.3.

Preliminarily note that putting t ′ = t -1 it is ⌊(t ′ + 1)/2⌋ = ⌊(t/2)⌋ so that by Formula (6.2): n(t ′ , ⌊(t ′ + 1)/2⌋) = t ′ + 1 + g t ′ -1 (6.7) equivalent to n(t -1, ⌊(t/2)⌋) = t + g t-2

Moreover, it is by Formula (6.3) for t > 4, or by Formula (6.2) for t = 3 and 48 . n(t, 2) = t + 2 + g t-1g t-3 = t + 2 + g t-2

and immediately: n(t, 2)n(t -1, ⌊t/2⌋) = 2 (6.8) Remark 6.1. Jointly Prop. 4.2 and 4.3 imply that, putting in lexicographic order the set Ξ of values of n(t, z), z > 1, a strictly monotone sequence is obtained. By Formulas (6.4), (6.6) and (6.8), all members of the sequence of first differences ∆ of the Ξ sequence are Fibonacci numbers (see Table 3).

  3):

Table 1 :

 1 helps to understand the structure of Fibonacci games for a set of small values of t. Structure of k(t, z) vectors for a set of small values of t.

	t	Φ(t)	k(t, 1)	n(t, 1)	k(t, 2)	n(t, 2)	k(t, 3)	n(t, 3)	k(t, 4)	n(t, 4)
	1	1	(3)	3
	2	1	(3, 1)	4
	3	2	(2, 2, 1)	5	(4, 1, 1)	6
	4	2	(2, 1, 2, 1)	6	(5, 1, 1, 1)	8
	5	3	(2, 1, 1, 2, 1)	7	(5, 2, 1, 1, 1)	10	(7, 1, 1, 1, 1)	11
	6	3	(2, 1, 1, 1, 2, 1)	8	(7, 1, 2, 1, 1, 1)	13	(10, 1, 1, 1, 1, 1)	15
	7	4	(2, 1, 1, 1, 1, 2, 1)	9	(10, 1, 1, 2, 1, 1, 1)	17	(13, 2, 1, 1, 1, 1, 1)	20	(15, 1, 1, 1, 1, 1, 1)	21
	8	4	(2, 1, 1, 1, 1, 1, 2, 1)	10	(15, 1, 1, 1, 2, 1, 1, 1)	23	(20, 1, 2, 1, 1, 1, 1, 1)	28	(23, 1, 1, 1, 1, 1, 1, 1)	30
	...	...	...	...	...	...	...	...	...	...

Table 2

 2 shows the behavior of the function n(t, z).

				z
	t	1	2	3	4	5
	1	3
	2	4
	3	5	6
	4	6	8
	5	7	10	11
	6	8	13	15
	7	9	17	20	21
	8	10	23	28	30
	9	11	32	40	43	44
	10	12	46	59	64	66
	...	...	...	...	...	...

Table 2 :

 2 List of the values n(t, z).

Table 3 :

 3 At row t and column z the entries correspond to g t-2z .

			z
	t	1	2	3	4
	3	1
	4	2
	5	3	1
	6	5	2
	7	8	3	1
	8	13	5	2
	9	21	8	3	1
	10	34	13	5	2
	...	...	... ... ...

Henceforth the subscript of the vector is its dimension, not to be confused with the subscript of a scalar which denotes a component of the vector.

As usual ⌊x⌋ denotes the floor(x). In particular, for t integer even (odd), ⌊(t + 1)/2⌋ = t/2 (or (t + 1)/2).

The first component of the sequence corresponds to the bottom type, i.e. to the weakest non dummy players.

It is the counterpart of the "largest min win" coalition (in backward lexicographic order) of a weighted majority game defined in[START_REF] Rosenmüller | Homogeneous games: recursive structure and Mathematics of Operations Research[END_REF], pp. 312.

By induction: it is true for j = 0, and if true for j it is true for j + 1.

Note that for t = 3 or 4 it is 1 + gt-2 = gt-1. Hence for such t: t + 2 + gt-2 = t + 1 + gt-1 so that both (6.2) and (6.3) describe n(t, ⌊(t + 1)/2⌋).

Now to prove Property 5.5 suppose on the contrary k j 0 > 2 and consider the coalition S ′ obtained from S by replacement of the player type j 0 + 1 with two additional players of type j 0 . S ′ would be minimal winning with w(S ′ ) > q, a contradiction.

Property 5.6. In any Fibonacci game it is k t-h = 1 for h even < t -1.

Proof. Proof of Property 5.6 is based on the previous Lemma 5.1 and on the following one:

Lemma 5.2. Suppose in a Fibonacci game there is a minimal winning coalition S and a j 0 = (t-h) > 2 such that s j = 0 for j = 1, . . . , j 0 -1, s j 0 > 0 (so that a player of type j 0 > 2 is the weakest in S), then k j 0 -1 = 1.

Proof. Suppose on the contrary k j 0 -1 > 1 and consider the coalition S ′ obtained from S by replacement of a player of type j 0 with two players of type j 0 -1. S

′ would be minimal winning with w(S ′ ) > q, a contradiction. This Lemma holds independently from the parity of h. Lemma 5.3. If a player of type j = th > 2 is the weakest in a minimal winning coalition S, then there exists a minimal winning coalition S ′ such that a player of type j -2 is the weakest in S ′ .

Proof. In S replace the player of type j with the unique (by Lemma 5.2) player of type j -1 and one player of type j -2. Also this Lemma holds independently from the parity of h.

Lemma 5.4. Suppose in a Fibonacci game k j 0 = 2 for j 0 = th with h odd < (t -1), then: a) k j 0 -1 = 1 (by Property 5.6) and b) there is a minimal winning coalition S ′ such that the player of type (j 0 -1) is its weakest player.

Proof. (of part b)) S

′ is obtained from the coalition S in Lemma 5.1 by replacement of the player of type t -(h -1) with one additional player of type th and the player of type t -(h + 1). Clearly S ′ is minimal winning with w(S ′ ) = q.

Lemma 5.5. If in a Fibonacci game k j 0 = 2 for j 0 = th with h odd < (t -1), then k j = 1 for any 1 < j < j 0 .

Proof. Lemma 5.4 and Lemma 5.3 imply that there exists a set of minimal winning coalitions such that each player of type j, 1 < j ≤ j 0 , is the weakest in a minimal winning coalition of the set. After that Lemma 5.5 comes as an immediate corollary of Lemma 5.2. Now, it is possible to give a more precise statement of Property 5.5 concerning the behaviour of k t-h for h odd < t -1. Indeed, Property 5.5 and Lemma 5.5 give immediately: Property 5.7. In a Fibonacci game there is at most one odd h < t -1 such that k j 0 = k t-h = 2 and k j = 1 for all j > 1 and j = j 0 . Now by Remark 2.1 it is:

k j g j (5.1)