Laurent Gilquin
email: laurent.gilquin@inria.fr

Clémentine Prieur

Elise Arnaud

Replication procedure for grouped Sobol' indices estimation in dependent uncertainty spaces

Keywords: sensitivity analysis, grouped Sobol' index, dependence, replicated designs, simplex

This paper deals with the estimation of grouped Sobol' indices in the case where the dependence inside each group is given by sets of linear ordered constraints. In the framework of independent inputs, the replication method allows to estimate first-order indices with only two designs. Through the use of orthogonal arrays of strength two the replication method can be applied to estimate closed second-order indices. We extend this methodology to estimate first-order and closed second-order grouped Sobol' indices under sets of linear constraints. The construction of the two designs required by the replication method is now based on the simplex geometric structure to handle the constraints within each set. We propose a space filling strategy to construct these designs.

Introduction

Sensitivity analysis studies how the uncertainty on an output of a mathematical model can be attributed to sources of uncertainty among the inputs. There are two main classes of sensitivity analyses called local and global sensitivity analysis. The former addresses sensitivity relatively to a nominal value of a given parameter. The latter examines sensitivity on the whole set of variation of the parameter. Global sensitivity analysis of complex and expensive mathematical models is a common practice to identify influent inputs and detect the potential interactions between them. Among the large number of available approaches, the variance-based method introduced by Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF] allows to calculate sensitivity indices called Sobol' indices. Each index gives an estimation of the influence of an individual input or a group of inputs. These indices give an estimation of how the output uncertainty can be apportioned to the uncertainty in the inputs. One can distinguish first-order indices that estimate the main effect from each input or group of inputs from higher-order indices that estimate the corresponding order of interactions between inputs. Closed k-th order indices estimate k-th order interactions in addition to the main effect of each of the k inputs. This estimation procedure requires a significant number of model runs, number that has a polynomial growth rate with respect to the input space dimension. This cost can be prohibitive for time consuming models and only a few number of runs is not enough to retrieve accurate informations about the model inputs. Saltelli [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF] proposed an improvement to reduce the number of runs but the total cost still depends linearly on the dimension of input space.

The notion of replicated designs to estimate first-order Sobol' indices probably goes back to McKay [START_REF] Mckay | Evaluating prediction uncertainty[END_REF] and appears later in Mara et al. [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF]. These last authors combine replicated designs with "pick-freeze" estimators [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF] to estimate first-order Sobol' indices. The procedure in Mara et al. [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF] has the major advantage of reducing drastically the estimation cost as the number of runs becomes independent of the input space dimension. This procedure has been further deeply studied and generalized in Tissot et al. [START_REF] Tissot | Estimating Sobol's indices combining Monte Carlo estimators and Latin hypercube sampling[END_REF] to the estimation of closed second-order indices. The generalization to closed second-order Sobol' indices relies on the replication of randomized orthogonal arrays (see Lemieux [6] or Owen [START_REF] Owen | Orthogonal arrays for computer experiments integration and visualization[END_REF]).

The motivation of our paper is to extend this methodology in presence of dependent inputs. Indeed, the case of correlated parameters has to be tackled with caution, as the calculation of single input indices does not provide anymore a proper information, that can be easily interpreted. One strategy is thus to define grouped indices for groups of correlated variables as proposed in Jacques et al. [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]. However, up to our knowledge, the problem of estimating those grouped indices, with a reasonable computing cost, has not been addressed in the literature. In this paper, we address this issue by proposing an approach based on replicated designs and randomized orthogonal arrays that enables to take into account dependency within inputs. We suppose that this dependency can be expressed through constraints. This approach can be used facing any set of constraints at the condition that one is able to provide points in the input space that verify the considered constraints. Guided by our application on a land-use and transport integrated model (LUTI) where some economical parameters are linked by order relations, this paper focus on the case of sets of linear ordered constraints. Thus we propose a sampling strategy based on the simplex geometric structure, that ensures a proper input space filling.

This paper is organized as follows. Section 2 recalls the notion of grouped Sobol' indices (see Jacques et al. [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]) to study the sensitivity of an output to groups of inputs where variables within a group are dependent but variables of different groups are independent. In Section 3 we review the use of replication method and randomized orthogonal array to estimate first and closed secondorder Sobol' indices. Section 4 is devoted to the sampling strategy adapted to our sets of linear ordered constraints and section 5 gives a summary of the whole procedure. Numerical experimentations are conducted in section 6. Classical case studies are addressed to compare our method to the ones in Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF] or Saltelli [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF].

Definition of Sobol' indices

Sobol' indices for independent inputs

Consider the following model defined from a black box perspective:

f : R d → R x = (x 1 , . . . , x d ) → y = f (x) ( 1
)
where y is the output of the model f , x the input vector. Let (Ω, A, P) be a probability space. We model the uncertainty on the inputs by a random vector X = (X 1 , . . . , X d ) whose components are independent. Let P X = P X 1 ⊗. . .⊗P X d denote the distribution of X. We assume that f ∈ L 2 (P X ). The model f can then be uniquely decomposed into summands of increasing dimensions (functional ANOVA decomposition [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF][START_REF] Hoeffding | A class of statistics with asymptotically normal distributions[END_REF]):

f (X) = f 0 + i f i (X i ) + i<j f ij (X i , X j ) + • • • + f 1...d (X 1 , . . . , X d ) , (2)
where each component of the decomposition verifies:

f i 1 ...is (x i 1 , . . . , x is )dP X i k (x i k ) = 0, ∀k ∈ {1, . . . , s}, ∀i 1 , . . . , i s ∈ {1, . . . , d}.
This implies that f 0 = E[Y ] and that the components are mutually orthogonal with respect to P X . Let I ⊆ {1, . . . , d}, each component is defined by:

f I (X I ) = E[Y |X I ] - J⊂I f J (X J ).
The functional decomposition can be used to measure the global sensitivity of the output Y to the input X i . By squaring and integrating (2), due to orthogonality we get:

V = i V i + i<j V ij + • • • + V 1,...,d
where

V I = Var[f I (X I )] = Var[E[Y |X I ]] - J⊂I V J and V = Var[Y ] .
Resulting from this decomposition, the Sobol' indices are defined by:

S I = V I V .
Let |I| denote the cardinal of I. This index measures the contribution to V of |I| th -order interaction between the X i , i ∈ I. Closed Sobol' indices are defined by:

S I = Var[E[Y |X I ]] V .
The closed Sobol' index S I measures the contribution of the X i , i ∈ I by themselves or in interaction with each other.

As an example, if there exist i = j ∈ {1, . . . , d} such that I = {i, j}, then

S ij = S ij + S i + S j .
At last, note that

Sobol' indices for dependent inputs

Consider again the model [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF]. We now suppose that X is a vector of d inputs among which one or multiple groups of variables are correlated. Variables belonging to the same group are dependent but variables of different groups are independent. This implies that one variable can not appear in two groups. Each group is denoted with a multidimensional variable. X can then be defined as follows:

X = (X 1 , . . . , X I , X I+1 , . . . , X I+k 1 X I+1 , . . . , X I+1+Σ i-1 , . . . , X I+Σ i X I+i , . . . , X I+k 1 +k 2 +•••+k m-1 +1 , . . . , X d X I+m ) where Σ i = i l=0 k l , k 0 = 0, (3) each X I+i , i ∈ {1, . . . , m} contains k i variables.
When inputs are dependent, equation [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF] given by the FANOVA decomposition of f does not hold as the components of the decomposition are not mutually orthogonal. The correlation of variables within X I+1 implies that a part of the sensitivity of one of these variables is explained by the other correlated variables. We can still compute the closed Sobol' index associated to each variable (or subset) of X I+1 but its value is difficult to interpret. An alternative is to define grouped Sobol' indices as in Jacques et al. [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]. These indices are defined like in the classical case of independent inputs, but refer to the multidimensional variables. We consider X defined as in [START_REF] Mckay | Evaluating prediction uncertainty[END_REF]. First-order grouped Sobol' indices are defined by:

S j = Var[E[Y | X j ]] V , ∀j ∈ {1, . . . , I + m} .
If X j is one-dimensional, then X j = X j and S j has the same expression as in the independent case. If X j is multidimensional (case j ∈ {I +1, . . . , I +m}), X j = X I+l , l ∈ {1, . . . , m}, the index is:

S j = S I+1+Σ l-1 ,...,I+Σ l = Var[E[Y |X I+1+Σ l-1 ,...,I+Σ l ]] V .
This index measures the impact of X j on the output Y . Higher order Sobol' indices can also be defined. In particular, second-order grouped Sobol' indices are given by:

S j,k = Var[E[Y | X j , X k ] -E[Y | X j ] -E[Y | X k ]] V .
This index measures the impact of the interaction between X j and X k on the output Y . The grouped Sobol' indices allow us to estimate the sensitivity of our model to the sets of linear ordered constrained inputs and the remaining unidimensional independent inputs.

standard estimation of Sobol' indices

We define by design a n × d matrix which consists in n evaluations of each X i , i = 1, . . . , d. The standard estimation of a closed Sobol' index S I proposed by Sobol' requires 2n evaluations of the model through two designs [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF]. Denoting by I c the complement of I ⊆ {1, . . . , d}, the second design of n evaluations is created from the first one by re-sampling the columns indexed by I c . This method is referred as the Sobol' method. Let Y and Y I be the vectors of evaluations associated to those two designs. The expression of S I can be rewritten as follows:

S I = E[Y Y I ] -E[Y ]E[Y I ] Var[Y ] .
An efficient estimator of S I with good asymptotic properties [START_REF] Janon | Asymptotic normality and efficiency of two Sobol' index estimators[END_REF][START_REF] Monod | Uncertainty and sensitivity analysis for crop models[END_REF] is:

S I = 1 n n j=1 Y j Y j I -1 n n j=1 Y j 1 n n j=1 Y j I 1 n n j=1 (Y j ) 2 -1 n n j=1 (Y j ) 2 . ( 4
)
Note that this method can be applied to independent inputs or independent groups of inputs. Its main drawback is the increasing number of needed evaluations. Estimating all first-order, respectively all closed second-order, Sobol' indices requires n(d + 1), respectively n( d 2 + 1), evaluations. For some models where a single evaluation takes a few minutes to one hour, this solution becomes rapidly intractable in case of large input space dimension.

Some improvements have been introduced by [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF] to reduce the number of evaluations but with a cost still depending of the dimension d of the input space. A solution to this issue lies in the use of the replication method as explained in the next section.

Replication method

The replication method has been described for first-order Sobol' indices estimation by Mara et al. [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF]. It has been combined with the use of randomized orthogonal array for closed second-order Sobol' indices estimation in Tissot and Prieur [START_REF] Tissot | Estimating Sobol's indices combining Monte Carlo estimators and Latin hypercube sampling[END_REF]. We propose hereafter a description of these methods and their generalization to independent groups of dependent variables.

Replication method for first-order Sobol' indices estimation

In their paper [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF], Mara et al. propose a comparison between a deterministic sensitivity analysis and a stochastic sensitivity analysis. The deterministic approach is achieved using the improved FAST method. The stochastic methodology is performed with a modification of the Sobol' method using the concept of replication. This improved Sobol' method consists in two replicated Monte Carlo designs in the case of d independent inputs. As described below, this method can be extended for groups of independent inputs.

Let us consider p independent one-dimensional or multidimensional input variables X 1 , . . . , X I , X I+1 , . . . , X p as in [START_REF] Mckay | Evaluating prediction uncertainty[END_REF]. Thus p = I + m.

Let

D = {X j = (X j 1 , . . . , X j I , X j I+1 , . . . , X j p ), 1 ≤ j ≤ n} be a design where X j i denotes the j-th sample of X i . Each X j is therefore a point in the d-dimensional input space. We say that a design D ′ is replicated from D if it is obtained through a column-wise permutation of D. More precisely it means that:

D ′ = {X ′j = (X π 1 (j) 1 , . . . , X π I (j) I , X π I+1 (j) I+1 , . . . , X πp(j) p ), 1 ≤ j ≤ n}
where the vectors π i i ∈ {1, . . . , p} are p independent random variables uniformly distributed in the set of permutations Π n and π i (j) = π j i denotes the j-th component of π i .

Let us denote by f the model and by Y D and Y D ′ the two vectors of model outputs associated to D and D ′ :

Y D = {Y j D = f (X j 1 , . . . , X j I , X j I+1 , . . . , X j p ), 1 ≤ j ≤ n} , Y D ′ = {Y j D ′ = f (X π 1 (j) 1 , . . . , X π I (j) I , X π I+1 (j) I+1 , . . . , X πp(j) p ), 1 ≤ j ≤ n} .
To estimate the index S k , k ∈ {1, . . . , p}, the values of Y D are rearranged with the corresponding permutation π k . As a result, it looks like Y D ′ has been constructed by varying all groups of inputs except the k-th. This is the same concept of "freezing" as the one introduced by Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF](see Appendix B for a detailed example).

Then using (4) with Y D rearranged and Y D ′ in place of Y and Y I we obtain S k (see Appendix B for an illustration of the method).

Replication method for closed second-order Sobol' indices estimation

In the case of closed second-order Sobol' indices estimation, the replication method does not consist anymore into a column-wise set of permutations. Indeed, to estimate closed second-order Sobol' indices, a way to "freeze" each subset of two variables has to be found. We need to find a structure allowing to get the same 2-sets of values in each subset of two columns of each replicated design. Such a structure has been introduced by Kishen [START_REF] Kishen | On latin and hyper-graeco cubes and hypercubes[END_REF] and further extended by Rao [START_REF] Rao | Hypercubes of strength "d" leading to confounded designs in factorial experiments[END_REF] and is known as orthogonal array. It is defined as follows:

Definition 1. A t -(q, d, λ) orthogonal array (t ≤ d) is a λq t × d matrix
whose entries are chosen from a q-set of N such that in every subset of t columns of the array, every t-subset of points of this q-set appears in exactly λ rows.

From this definition by setting t = 2, we can construct a structure consisting of points in {1, . . . , q} λq 2 where each 2-set of columns have the same 2-set of points λ times.

We present now a result from the method of differences introduced by Bose and Bush [START_REF] Bose | Orthogonal arrays of strength two and three[END_REF] to construct a 2 -(q, p, 1) linear orthogonal array. We begin with some useful definitions and a theorem resulting from the method of differences: Definition 2. A linear orthogonal array is a t -(q, d, λ) orthogonal array where λ = 1 and q is a prime number. The set of rows of a linear orthogonal array is a subspace of (Z/qZ) d . Definition 3. Let G be the Galois field of order q, GF (q), a λq × d matrix D(λq, d, q) is a difference matrix (or difference scheme) of strength 2 if for every subset of 2 columns of the matrix, the vector of row-wise difference of the subset contains every elements of GF (q) λ times.

Theorem 1. Let M be a D(λq, d, q) difference matrix of strength 2, B the vector of elements of GF (q), ⊕ the Kronecker sum and T the transpose operator then (M ⊕ B) T is a 2 -(q, d, λ) orthogonal array. Now, consider q a prime number and G = GF (q) the Galois field associated, the multiplication table of G is a D(q, q, q) difference matrix of strength 2 with λ = 1. From Theorem 1 we can construct a 2 -(q, p, 1) linear orthogonal array: A = (M ⊕ B) T where M is a sub-matrix of p columns of the multiplication table of G where p still denotes the number of groups of independent inputs. What is left is to connect this 2 -(q, p, 1) orthogonal array to the replication method. We consider a q × p design:

D = {X j = (X j 1 , . . . , X j I , X j I+1 , . . . , X j p ), 1 ≤ j ≤ q}
where X 1 , . . . , X I , X I+1 , . . . , X p are one-dimensional or multidimensional independent variables and q is a prime number. Denote by A the previous 2 -(q, p, 1) orthogonal array. We can establish a direct bijection between the indices {1, . . . , q} of rows of D and the set {1, . . . , q} from which A is constructed. Each column of A will serve as new indices to build the respective column of the new design. Resulting from the construction procedure, in each column of A, each entry of {1, . . . , q} appears q times. This is due to the fact that each entry is involved in q different 2-sets. The number of rows of A is equal to q 2 meaning that using the previous bijection we can construct a design D OA of q 2 rows where each column has the same entries of the corresponding column of D but re-indexed from the elements of the corresponding column of A:

D OA = {X j = (X A j 1 1 , . . . , X A j I I , X A j I+1 I+1 , . . . , X A j p p ), 1 ≤ j ≤ q 2 }
. This gives us the first of the two needed replicated designs. To replicate D OA we consider once again random vectors of permutations π i in Π q . We also have a direct bijection between the elements {1, . . . , q} of each π i and the set {1, . . . , q} from which A is constructed. This allows the following construction of the second needed replicated design:

D ′ OA = {X ′j = (X π 1 (A j 1 ) 1 , . . . , X π I (A j I ) I , X π I+1 (A j I+1 ) I+1 , . . . , X πp(A j p ) p ), 1 ≤ j ≤ q 2 }
where each π i is re-indexed from the elements of the corresponding column of A. In fact, the matrix whose columns are the vectors π 1 (A 1 ), . . . , π p (A p ) is also an orthogonal array (see Appendix A for details).

We obtain two designs of length q 2 . Denote then Y D OA and Y D ′ OA the two vectors of model outputs associated to D OA and D ′ OA . To get the same "freezing" concept as in the last subsection, but this time for each subset of two coordinates (k, l) ∈ {1, . . . , p} 2 , we apply the following transformation:

∀j ∈ {1, . . . , q 2 } : Y A j k +A j l * q-q = Y D OA j , Y ′ π k (A j k )+π l (A j l ) * q-q = Y D ′ OA j .
As a result, Y and Y ′ have been constructed by varying all groups of inputs except the k-th and l-th. This allows us to estimate the index S kl with the formula ( 4) by replacing Y I with Y ′ (see Appendix C for a detailed example).

Remark 1. The construction of a 2 -(q, d, 1) linear orthogonal array is possible only if q ≥ d -1. The total number of evaluation required by the replication method will be 2q 2 .

Remark 2. For the sake of clarity, we have treated the case of first-order and closed second-order indices estimation separately. Actually, given that a 1 -(q, d, 1) orthogonal array is a matrix of permutations we could have written a general definition for both cases.

Sampling strategy under linear ordered constraints

In the previous section a method based on replicated designs to estimate Sobol' indices for independent groups of inputs was proposed. The fact that it handles groups of inputs is a way to tackle the problem of dependence within parameters. Indeed, each set of dependent inputs is treated as a group. The first-and second-order indices are calculated for the group, and not independently for each parameter (which would provide values that can not be interpreted easily). This method can be used for any kind of dependency, under the condition that one is able to properly sample from the input space; samples that are further needed to feed the designs. When inputs are independent, getting samples uniformly distributed in [0, 1] is sufficient. When dependency is expressed through constraints, one has to create samples that satisfy the constraints, which is not an easy task in general. As mentioned in the introduction, the work presented in this paper is motivated by a variancebased sensitivity analysis we want to conduct on a land use and transport integrated model, called Tranus [START_REF] Barra | Mathematical description of TRANUS[END_REF]. Such a model creates a numerical representation of the transport network on a territory, and calculates the spatial repartition of the economic activities. It is based on economic laws whose parameters can be linked though ordered linear constraints. We thus focus on such constraints and propose an adapted sampling strategy based on the simplex geometric structure. This strategy is further enhanced with a space filling methodology.

Sampling strategy based on the simplex to handle constraints

We consider the same vector of inputs X as defined in [START_REF] Mckay | Evaluating prediction uncertainty[END_REF]. Each unidimensional variable is uniformly distributed in [0, 1]. Each multidimensional variable X I+i , i ∈ {1, . . . , m} contains k i variables valued in [0, 1] and is subject to the following linear ordered constraints:

X I+1+Σ i-1 ≤ X I+2+Σ i-1 ≤ • • • ≤ X I+Σ i . (5)
In case of absence of constraints, the joint distribution of each multidimensional variable X I+i is a uniform distribution on the unit k i -hypercube. The introduction of constraints ( 5) transforms the support of the joint uniform distribution in a k i -simplex. The general definition of a simplex is the following:

Definition 4. A simplex, or k-simplex, is the generalization of a triangle (2-simplex) or tetrahedron (3-simplex) to any dimension k. Let ∆ k be the simplex with vertices P 0 , . . . , P k in R k :

∆ k = {M = k j=0 w j P j | w j ≥ 0, k j=0 w j = 1}
To sample from a uniform distribution over a simplex we refer to the following theorem of L. Devroye [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF]:

Theorem 2 (Devroye,1986). Let ∆ k be a simplex in R k with vertices P 0 , . . . , P k and U 1 , . . . , U k be independent variables uniformly distributed in [0, 1]. We define W 0 , . . . , W k the associated differences as follows:

∀j ∈ {1, . . . , k}, W j = V j+1 -V j , V 0 = 0, V k+1 = 1
where V j denote the j-th ordered statistics constructed from the set (U i ) i : ∀j ∈ {1, . . . , k}, V j = min i=1,...,k

U i \ V j-1 .
Then the random variable Z = k j=0 W j P j is uniformly distributed in the simplex ∆ k .

To each set of constraints ( 5) is associated a k i -simplex, defined as follows:

∆ * k i = {M = k i j=0 w j e * j | w j ≥ 0, k i j=0 w j = 1}
with vertices the ordered vectors of the standard basis of R k i :

e * 0 = (0, 0, . . . , 0, 0), e * 1 = (0, 0, . . . , 0, 1), . . . , e * k i = (1, 1, . . . , 1, 1)

Since the vertices of ∆ * k i satisfy ordered linear constraints (5), then every point sampled in ∆ * k i will also satisfy those constraints. Sampling in ∆ * k i using the strategy of Theorem 2 is equivalent to sampling uniformly under a set of k i linear ordered constraints [START_REF] Tissot | Estimating Sobol's indices combining Monte Carlo estimators and Latin hypercube sampling[END_REF]. This method enables to sample each of the m multidimensional variables contained in the vector of inputs X. For the remaining independent one-dimensional variables we use uniform random numbers. Algorithm 1 summarizes the whole sampling strategy to get n samples:

Space filling strategy

In the classical case of independent variables, one should prefer to use a Latin Hypercube rather than sampling each variable from a uniform distribution. Such a procedure ensures that the input space will be properly explored. To ensure the same property in our case of constrained input space, the proposed sampling method has to be improved as nothing prevents all the points to be located in a small area of the space. A way to do so is to first subdivide each k i -simplex into sub-simplices of order k i and then sample from the latter. Algorithm 2 summarizes our subdivision method (see Appendix D for a detailed example). For the sake of clarity, in the following, k i has been replaced by k. Figure 1 shows two examples of subdivisions for k = 2 and k = 3. Let us remark that from Theorem 2, the uniform distribution in the unit ordered k-simplex in R k corresponds to the joint distribution of the k ordered statistics. The l-th ordered statistic of a uniform distribution is a Beta variable : X (l) ∼ B(l, k + 1l).

This remark can serve as a goodness of fit measure for our space filling sampling improvement of the uniform distribution over the unit ordered ksimplex. Indeed, since the marginal distributions are known, a Kolmogorov Smirnov test (ks.test) can be performed to compare our space filling strategy to a standard uniform sampling in the simplex. For each l-th ordered statistic, we perform the test for k = 3 and different sizes n of the sample. The test procedure is repeated r = 100 times. We take the average value of the r p-values obtained. Table 1 gives the averaged p-value obtained for each (n, l)-set. The general observation is that for each l-th ordered statistic the average p-value obtained for the space filling sampling is better than the one obtained for the standard sampling, whatever the number of samples is. This clearly shows the advantage of using the space filling strategy.

Algorithm 2 Space filling sampling of n points for a set of k linear ordered constraints 1. create the k 2 vertices of the unit k-hypercube and assign a number to each vertice according to the order of creation 2. identify the k + 1 vertices whose coordinates satisfy the set of constraints 3. store those vertices into a (k + 1) × k matrix M coor , each column corresponding to one of the k coordinates 4. create the k! vectors of permutations of Π k 5. apply each permutation to the columns of M coor , this creates a total of k! matrices 6. for each one of the k! matrices, match the vertices with those of the unit k-hypercube and retrieve the vector of matching numbers, this creates a total of k! vectors of matching numbers that we note ω j , j ∈ {1, . . . , k!} 7. choose α = floor(n

1 k ) (respectively α = n 1 k if n 1 k ∈ N) the number of levels 8. subdivide the unit k-hypercube into (α + 1) k (respectively α k if n 1 k ∈ N) sub-hypercubes.
For each sub-hypercube, assign a number to each vertices following the same order of numeration as in 1. 9. for each sub-hypercube:

• for each ω j , conserve the vertices whose indices correspond to ω j , this gives k! sets of vertices

• among those sets keep only those whose coordinates satisfy the constraints 10. randomly discard (α + 1) kn of all the sets created to keep only n sets 11. sample one point using each set in theorem 2 (each set are the vertices of a simplex) to get a total of n points

Summary of the whole procedure

Consider the vector of inputs

X = (X 1 , . . . , X I , X I+1 , . . . , X I+k 1 X I+1 , . . . , X I+1+Σ i-1 , . . . , X I+Σ i X I+i , . . . , X I+k 1 +k 2 +•••+k m-1 +1 , . . . , X d X I+m ) where Σ i = i l=0 k l , k 0 = 0 .
Each multidimensional variable X I+i , i ∈ {1, . . . , m} contains k i variables and is subject to the following linear ordered constraints:

X I+1+Σ i-1 ≤ X I+2+Σ i-1 ≤ • • • ≤ X I+Σ i .
Let d be the dimension of input space, K = m i=1 k i be the number of constrained inputs, m be the number of sets of constrained variables then the number of grouped Sobol' indices to estimate is p = d -K + m for first-order indices and p 2 for closed second-order indices. Let (A j i ) i=1..p,j=1..n be an orthogonal array of strength t. t = 1 respectively 2 to estimate first-respectively second-order indices. Let n = q t ∈ N * be the number of samples of the input space with q a prime number if t = 2. Denote (D OA j i ) i=1,...,d, j=1,...,n and (D ′ OA j i ) i=1,...,d, j=1,...,n the two replicated designs.

1. For the independent inputs X i i ∈ {1, . . . , I}, construct the replicated designs sampling from Latin Hypercube [START_REF] Tissot | Estimating Sobol's indices combining Monte Carlo estimators and Latin hypercube sampling[END_REF] to get space filling uniform distribution.

∀i ∈ {1, . . . , I} :

D OA j i = A j i -U A j i i q , D ′ OA j i = π i (A j i ) -U π i (A j i ) i q ,
where the π i and the U j i are independent random variables uniformly distributed in Π n and [0, 1].

2. For the constrained inputs X I+i i ∈ {1, . . . , m}, construct the replicated designs sampling from algorithm 2 to get a space filling uniform distribution in the corresponding ordered k i -simplices.

∀i ∈ {1, . . . , m} :

D OA j I+1+Σ i-1 ,I+2+Σ i-1 ,...,I+Σ i = k i l=0 W A j i l P A j i l , D ′ OA j I+1+Σ i-1 ,I+2+Σ i-1 ,...,I+Σ i = k i l=0 W π i (A j i ) l P π i (A j i ) l , where W A j i 0 , . . . , W A j i k i (respectively W π i (A j i ) 0 , . . . , W π i (A j i ) k i
) are the associated differences constructed from U

A j i I+1+Σ i-1 , . . . , U A j i I+Σ i (respectively U π i (A j i ) I+1+Σ i-1 , . . . , U π i (A j i )
I+Σ i ) and P 0 , . . . , P k i are the columns of the n×(k i +1) matrix P storing the sets of vertices resulting from Algorithm 2.

Denote then Y D OA and Y D ′

OA the two vectors of model outputs associated to D OA and D ′ OA we compute the Sobol' indices with the methodology presented in section 3.

Remark 3.

As we stated at the beginning of section 4, we can apply the replication method to any kind of distributions under any type of correlation. The only requirement is to be able to generate a design or in other words to be able to sample each multidimensional variable (see section 6 for a multivariate Gaussian example).

Application to test functions

In this section the replication method is applied to standard test functions with sets of linear ordered constraints and to a case of a multivariate normal distribution. We compare the results to those given by the Saltelli method and the standard method of Sobol' that we both generalized to the case of groups of inputs (a first generalization of the standard method of Sobol' can be found in [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]). For both the standard method and our replication method we used the efficient estimator with formula given in [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF]. The size N of the design of experiments (DoE) is different for each method. Let n denotes the numbers of points in the d-dimensional input space:

• for the replication method we get N = 2n (each of the two replicated designs has n points) for first-order indices and N = 2n for closed second-order Sobol' indices where √ n has to be a prime number

• for Saltelli's method we get N = n(p + 2) for first-order indices and N = n(2p + 2) for closed second-order indices where p is the number of independent one-dimensional or multidimensional input variables

• for the standard method we get N = n(p + 1) for first-order indices and N = n( p 2 + 1) for closed second-order indices From the user point of view, the total number of evaluations required is one of the most important aspect of the method. Thus, the comparison are made for similar computational time. For each method, we chose n to get the same value N .

Sobol's g-function

We consider the g-function introduced by Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF]. The g-function is defined as follows:

f (X 1 , . . . , X d ) = d i=1 g i (X i ), where g i (X i ) = |4X i -2| + a i 1 + a i , a i ≥ 0.
Each value a i determines the relative importance of the X i . When the value of a i gets closer to zero the variable X i becomes more influent. We choose here d = 4 and (a 1 = 0, a 2 = 1, a 3 = 3, a 4 = 6). We consider the following sets of characteristics:

• X 3 ≤ X 4 with (X 3 , X 4 ) uniformly distributed in ∆ * 2
• X 1 and X 2 are independent and both uniformly distributed in [0, 1].

Since f has an analytical expression, Sobol' indices can be calculated by symbolic integrals to evaluate their theoretical value. Each of the three methods are used to compute the three first-order Sobol' indices and the three closed second-order Sobol' indices. The estimation procedure is repeated r = 100 times to get a sample of estimated indices; and this for various values of N .

To evaluate the precision of each method, are shown for each index:

• boxplots drawn depending on increasing values of N ;

• evolutions of the error relative to the true value, depending on N . The error calculated is a variation of the SAE (Sum of Absolute Error) indicator used in [START_REF] Mara | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF], and is defined as:

SAE N = r k=1 | Ŝk,N -S|,
where Ŝk,N is the k-th estimation of the index, and S the theoretical value of the index.

Results on first-order indices are shown on figure 2, while results on closed second-order Sobol' indices are shown on figure 3. In both cases, the first row corresponds to boxplots and the second row to error curves. Note that for all illustrations, black is associated to our replication method, blue to Saltelli approach and green to Sobol' method (also referred as standard method). The theoretical values of indices are indicated by red horizontal broken lines.

For the calculation of first-order indices, the observation is that the replication method gives better results, and particularly for low values of N . This observation is explained by the independence between the computational cost of our approach and the input space dimension. Thus, this observation is emphasized when we increase the number p of independent groups of inputs. The Saltelli's method gives better results when it comes to estimate indices with low value rather than indices with high value. At the opposite both the standard method and the replication method are better at estimating indices with high value but this could be due to the choice of the estimator rather than from the method itself (see Owen [START_REF] Owen | Better estimation of small sobol' sensitivity indices[END_REF] for a discussion on the estimation of small indices).

The replication method performs also better compared to the other two approaches for the calculation of closed second-order indices. The Saltelli method is once again better at estimating indices with low value and the other two methods are better at estimating indices with high value. Despite the difference, for low value indices our method still stays competitive.

In order to test our method facing a high number of inputs, we consider the Sobol' g-function with a dimension of input space now equals to 50. The following characteristics are considered:

• (X 1 , X 50 ), (X 2 , X 49 ), (X 3 , X 48 ), (X 4 , X 47 ), (X 5 , X 46 ) uniformly distributed in ∆ * 2 ,
• (X 6 , X 10 , X 40 ), (X 7 , X 11 , X 39 ), (X 8 , X 12 , X 38 ), (X 9 , X 13 , X 37 ), (X 19 , X 30 , X 31 ) uniformly distributed in ∆ * 3 ,

• (X 14 , X 20 , X 21 , X 36 ), (X 15 , X 22 , X 23 , X 35 ), (X 16 , X 24 , X 25 , X 34 ), (X 17 , X 26 , X 27 , X 33 ), (X 18 , X 28 , X 29 , X 32 ) uniformly distributed in ∆ * 4 ,

• (X 41 , X 42 , X 43 , X 44 , X 45 ) uniformly distributed in ∆ * 5 .

The values of the a i are the following: (a 1 = ... = a 10 = 0, a 11 = ... = a 20 = 1, a 21 = ... = a 30 = 2, a 31 = ... = a 40 = 4, a 41 = ... = a 50 = 6).

Since there is a substantial number of indices, we choose not to draw all boxplots but instead to compute the two errors. The average of the p SAE N calculated for each first-order index is referred as ASAE 1 N . The average of the p 2 SAE N calculated for each closed second-order index is referred as ASAE 2

N . Table 2 gives the values of ASAE 1 N for each method and each value of N . Figure 4 right represents the evolution curves of the ASAE 2 N . Since there is 120 closed second-order indices to estimate, the difference in that case between the results of the three methods is even more noticeable. Figure 4 left represents the evolution curves of the ASAE 1 N given N for the three methods. It clearly illustrates the benefit of the independence of the replication method's total cost to the dimension of the input space. Even for N = 10 4 the replication method is still in average two times better than the other two methods.

f (X 1 , . . . , X d ) = d i=1 (-1) i i j=1 X j .
We consider the following characteristics:

• X 3 ≤ X 4 with (X 3 , X 4 ) uniformly distributed in ∆ * 2
• X 1 and X 2 are independent and both uniformly distributed in [0, 1].

We use the same estimation procedure as in the last subsection and we draw the same graphical results (Figure 5 and Figure 6). Conclusions are similar to the ones for the g-function of Sobol'. It shows the efficiency of our approach with respect to more standard approaches such as Saltelli or Sobol'.

Correlated multivariate normal distribution

As stated in Remark 3, this example of a correlated multivariate normal distribution is to shows that the replication method can be used in a general case of correlated inputs at the condition that we know how to sample each multidimensional input. Furthermore, the case of multivariate normal distribution with correlation terms is a common example of dependence for the inputs.

Let Σ be a positive-definite matrix, non necessarily diagonal. As in the previous subsection we write X = (X 1 , . . . , X I , X I+1 , . . . , X I+m ) where the multidimensional variables X I+j , j ∈ {1, . . . , m} are independent random gaussian vectors with mean vector µ and covariance matrix Σ. Based on section 2.2 we can estimate the Sobol' indices for the groups of correlated variables X I+i .

To sample such a distribution we first use a Cholesky decomposition of the covariance matrix Σ = LL T . Then, let Z = (Z 0 , . . . , Z d ) be a vector whose components are independent and follow a standard normal distribution. We obtain the true sampling from the formula: X = µ + LZ. We obtain the two replicated designs of experiments D OA and D ′ OA with the formula from section 3 applied to our previous sample of X. As an example we consider the following function inspired from the function presented in [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF]: The first two variables X 1 ,X 2 are correlated and we chose the following values for the terms in Σ:

f (X 1 , X 2 , X 3 , X 4 ) = g 1 (X 1 , X 2 ) + g 2 (X 2 ) + g 3 (X 3 ) + g 4 (X 4 ) where        g 1 (X 1 , X 2 ) = (2(X 1 -µ 1 ) + 1)(3(X 2 -µ 2 ) + 2) g 2 (X 2 ) = 2(X 2 -µ 2 ) 2 + X 2 -µ 2 + 3 g 3 (X 3 ) = 1 + 2(X 3 -µ 3 ) + 2(X 3 -µ 3 ) 2 + 3(X 3 -µ 3 ) 3 g 4 (X 4 ) = 1 + 4(X 4 -µ 4 ) and X = (X 1 , X 2 , X 3 , X 4 ) follows a centered multivariate normal distribution N (µ, Σ) with parameters: Σ =     σ 2 1 ρ 12 σ 1 σ 2 0 0 ρ 12 σ 1 σ 2 σ 2 2 0 0 0 0 σ 2 3 0 0 0 0 σ 2 4     µ =     0 0 0 0     .
σ 1 = σ 2 = 0.3, σ 3 = 0.5, σ 4 = 0.6, ρ 12 = 0.2
We use the same estimation procedure as in the last two subsections and we draw the results only for first-order Sobol' indices since the sum of the true values of those indices is equal to 0.999 which means that interactions are negligible.

Looking at the results in Figure 7 our replication method still gives good results even in this case where there is no space filling strategy involved. The results given by the Saltelli method are in accordance with the results from the previous two subsections.

Conclusion

In this paper we proposed a methodology to estimate first-order or closed second-order grouped Sobol' indices using replicated designs to handle groups of dependent inputs. In our application on the LUTI model of Tranus this dependency translates into groups of linear ordered inputs. Thus, through an algorithm (Algorithm 2) we have presented a space filling strategy based on the simplex structure to sample points satisfying those constraints and a detailed procedure to estimate the associated Sobol' indices. Compared to the standard method of Sobol' [START_REF] Sobol | Sensitivity indices for nonlinear mathematical models[END_REF] or to Saltelli's approach [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF], our new estimation procedure has a better precision for a given cost, especially in the case of high number of inputs. An application on the Tranus LUTI model will be the subject of a future work via a grid deployment. given that (x, x ′ ) ∈ {1, . . . , q} × {1, . . . , q} we obtain x = x ′ .

Then (y ′y) * q = 0 implies y = y ′ .

Since A j 1 and A j 2 ∈ {1, . . . , q} (resp. π 1 (A j 1 ) and π 2 (A j 2 )), this means that every element of {1, . . . , q 2 } is obtained from a unique 2-set of (A 1 A 2 ) (resp. (π 1 (A 1 ) π 2 (A 2 )) thus each value of the vectors Y and Y ′ are calculated from the same set of the first two coordinates. We get the same 'freezing' result for every other set of two coordinates.

Algorithm 1

 1 Sampling n points from d inputs under m disjoint k i -sets of linear ordered constraints 1. for each of the m multidimensional variables X I+i • create the unit ordered k i -simplex • sample n points with Theorem 2 2. for each of the d -m i=1 k i one-dimensional variables sample n random uniform numbers.

Figure 1 :

 1 Figure 1: examples of subdivision of the unit ordered k-simplex

Figure 2 :

 2 Figure 2: Estimation of first-order Sobol' indices given by the three methods for r = 100 repetitions. The color black is for the replication method, the blue for the Saltelli method and the green for the standard method. At the top: boxplot representation for different values of N . At the bottom: curve of the SAE N for different values of N .

Figure 3 :

 3 Figure 3: Estimation of closed second-order Sobol' indices given by the three methods for r = 100 repetitions. The color black is for the replication method, the blue for the Saltelli method and the green for the standard method. At the top: boxplot representation for different values of N . At the bottom: curve of the SAE N for different values of N . Sobol' index S {1}{2} Sobol' index S {1}{3,4} Sobol' index S {2}{3,4}

Figure 4 :

 4 Figure 4: Curves of ASAE 1 N (left figure) and ASAE 2 N (right figure) for each method depending on N

Figure 5 :

 5 Figure 5: Estimation of first-order Sobol' indices given by the three methods for r = 100 repetitions. The color black is for the replication method, the blue for the Saltelli method and the green for the standard method. At the top: boxplot representation for different values of N . At the bottom: curve of the SAE N for different values of N .

Figure 6 :

 6 Figure 6: Estimation of closed second-order Sobol' indices given by the three methods for r = 100 repetitions. The color black is for the replication method, the blue for the Saltelli method and the green for the standard method. At the top: boxplot representation for different values of N . At the bottom: curve of the SAE N for different values of N .

Figure 7 :

 7 Figure 7: Estimation of first-order Sobol' indices given by the three methods for r = 100 repetitions. The color black is for the replication method, the blue for the Saltelli method and the green for the standard method. At the top: boxplot representation for different values of N . At the bottom: curve of the SAE N for different values of N .

X 1 X 2 X 3 X 1 X 2 X 3

 123123 Appendix B. Example of the replication method for first-order indices D = X =   index S 1 we re-sample Y D with π 1 : , 22, 32, 42) f (13, 23, 33, 43) f (11, 21, 31, 41)   looking at both Y Dnew and Y D ′ , X 1 has been "frozen". Appendix C. Example of the replication method for closed secondorder indices D = X =   2 -(3, 4, 1) orthogonal array A = (M ⊕ B) T : add another column corresponding to the vector of elements we construct Y with this new indexation:                           Y [9] = Y D OA [1] = f (13, 23, 33, 43) Y [3] = Y D OA [2] = f (13, 21, 32, 41) Y [6] = Y D OA [3] = f (13, 22, 31, 42) Y [1] = Y D OA [4] = f (11, 21, 31, 43) Y [4] = Y D OA [5] = f (11, 22, 33, 41) Y [7] = Y D OA [6] = f (11, 23, 32, 42) Y [5] = Y D OA [7] = f (12, 22, 32, 43) Y [8] = Y D OA [8] = f (12, 23, 31, 41) Y [2] = Y D OA [9] = f (12, 21, 33, 42) construct Y ′ with this new index:                           Y ′ [4] = Y D ′ OA [1] = f (11, 22, 32, 43) Y ′ [1] = Y D ′ OA [2] = f (11, 21, 31, 41) Y ′ [7] = Y D ′ OA [3] = f (11, 23, 33, 42) Y ′ [2] = Y D ′ OA [4] = f (12, 21, 33, 43) Y ′ [8] = Y D ′ OA [5] = f (12, 23, 32, 41) Y ′ [5] = Y D ′ OA [6] = f (12, 22, 31, 42) Y ′ [9] = Y D ′ OA [7] = f (13, 23, 31, 43) Y ′ [6] = Y D ′ OA [8] = f (13, 22, 33, 41) Y ′ [3] = Y D ′ OA [9] = f[START_REF] Rao | Hypercubes of strength "d" leading to confounded designs in factorial experiments[END_REF] 21, 32, 42) retrieve the vectors of numbers of matching vertices: w 2 = 1 2 4 matrix of subdivision of the 2-hypercube into 3 2 = 9 sub-hypercubes (one sub-hypercube per row of the matrix): each sub-hypercube, for each w j conserve the vertices whose index corresponds to w j : over the list of sub-hypercubes, to get all the sets of vertices, then conserve only those satisfying the constraints:

Table 1 :

 1 Averaged p-values given by the ks.test to compare our space-filling sampling to the standard sampling on the unit ordered 3-simplex. The test is repeat r = 100 times for each l-th ordered statistic and for different value of number n of points sampled.

			standard sampling	space filling sampling
	n	l	1	2	3	1	2	3
	100		0.53	0.56	0.50	0.80	0.82	0.76
	250		0.54	0.53	0.47	0.82	0.82	0.75
	500		0.51	0.48	0.53	0.93	0.98	0.93
	1250	0.49	0.51	0.52	0.93	0.96	0.93
	2500	0.46	0.47	0.46	0.93	0.97	0.95
	3750	0.46	0.52	0.50	0.96	0.98	0.96
	5000	0.46	0.49	0.49	0.92	0.94	0.93

Table 2 :

 2 Values of ASAE 1 N for each method and for different size N of the design of experiments

	size	replicated	Saltelli	standard
	N	method	method	method
	200	4.80	54.39	16.62
	500	3.40	39.27	7.40
	1000	2.35	16.59	5.80
	2500	1.61	4.18	3.35
	5000	1.18	3.41	2.97
	7500	1.01	2.12	2.54
	10 4	0.93	2.02	1.78

I⊂{1,...,d},I =∅ S I = 1, allowing a direct interpretation of the value of each index.

Sobol' index S 1Sobol' index S 2 Sobol' index S 3,4

Sobol' index S 3 Sobol' index S 4 Sobol' index S 1,2

Appendix A. Proofs Proposition 1. Let A be a 2-(q, p, 1) linear orthogonal array of elements in {1, . . . , q}. Let π 1 (A 1 ), . . . , π p (A p ) ∈ Π q , p independent vectors of permutations. Then the matrix A π whose columns are the vectors π 1 (A 1 ), . . . , π p (A p ) is also a 2-(q, p, 1) orthogonal array.

Proof. From Definition 1, we need to demonstrate that for every subset of 2 columns of A π , every 2-set appears exactly one time. This would mean that every subset of 2 columns of A π is identical to the space {1, . . . , q}×{1, . . . , q}. Let us consider only A 1 , A 2 and π 1 , π 2 , the demonstration will be identical for the other subset of 2 columns. Consider the following application:

f is a bijection since π 1 , π 2 are two bijections from {1, . . . , q} to {1, . . . , q}. Now if x are elements of A 1 and y are elements of A 2 this means that the subset (π 1 (A 1 ) π 2 (A 2 )) have the same 2-set of (A 1 A 2 ). Applying this to every other subset of 2 columns leads to the result.

Proposition 2. We want to prove that the following transformation allows us to freeze each set of two coordinates (k, l) ∈ {1, . . . , p} 2 allowing us to estimate each closed second-order index S kl : ∀j ∈ {1, . . . , q 2 } :

Proof. To simplify the proof consider k = 1, l = 2, the method will be the same for the other set of two coordinates. Consider the following application:

f is a bijection. To prove it, we only have to demonstrate that f is injective. Consider (x, y) and (x ′ , y ′ ) in {1, . . . , q} × {1, . . . , q} such that:

creation of the vertices of the unit 2-hypercube and numeration:

vertices satisfying the constraints and creation of M coor :

vectors of permutations and creation of the 2! matrices by column-wise permutations of M coor :


Jana Trgalova

Ulrich Kortenkamp

Ana Paula Jahn

Paul Libbrecht

Christian Mercat

Tomas Recio

Sophie Soury-Lavergne

I2GEO.NET -A PLATFORM FOR SHARING DYNAMIC GEOMETRY RESOURCES ALL OVER EUROPE

Keywords: Intergeo, I2Geo platform, dynamic geometry

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

DGS are well known computer tools to support mathematics teaching and learning by means of personal explorations and experience. Despite their availability and the recommendations, or even requirements, in European countries curricula to use them, their integration in schools is still unsatisfactory [START_REF] Hendriks | Status quo report on DGS usage[END_REF]. The Intergeo project [START_REF] Kortenkamp | Interoperable Interactive Geometry For Europe -First Tech. and Educ. Results and Future Challenges of the Intergeo Project[END_REF] tackles main obstacles to DGS integration: (1) difficulties to find suitable resources due to the lack of metadata describing accurately their content, (2) the impossibility to exploit with a given DGS resources created with another one, (3) the lack of quality guarantees of available resources. We briefly expose and illustrate solutions proposed to overcome these obstacles.

SHARING DG RESOURCES

An accurate resource content description is necessary for easy finding and sharing resources. The challenge in Intergeo consisted in defining international metadata allowing a cross-curricular resource search regardless of language and cultural differences. This was achieved by defining a mathematical topics and competencies ontology making it possible to find resources even if they are written in different languages. Moreover, different sophisticated search tools (such as a search associated to key words from the different national curricula) have been developed.

INTEROPERABILITY OF DGS

One of the obstacles in the use of DG is the issue of user lock-in with respect to a particular software product. Learning new software is time consuming and usually unrealistic for teachers. The project defines a common interoperable file format that allows describing any construction created with a DGS and that can be used to interchange content between DGSs. The standards are supervised by a non-profit association, Intergeo A.s.b.l.

QUALITY ASSESSMENT PROCESS

Quality assessment of DG resources in the Intergeo project [START_REF] Trgalová | Quality process for dynamic geometry resources: the Intergeo project[END_REF] aims at promoting access to best quality resources, as well as at ensuring their continuous improvement. The main tool supporting DG resource quality assessment on the platform is a standardized questionnaire organized around nine dimensions of a resource related to its mathematical, didactical, pedagogical, technical and ergonomic aspects. Quantitative evaluation of a resource along these aspects in terms of a 4-level range of agreement can be complemented by qualitative comments, which are crucial for the resource improvement.

I2GEO PLATFORM AND TEACHER PROFESSIONAL DEVELOPMENT

Around 3000 resources are available on the platform and more than 500 evaluations have been performed to this day. The analysis of the resource reviews and of the interviews with reviewers provides evidence that involvement into a quality assessment process is an interesting means of life-long teacher professional development by reflexive analysis and the diffusion of best pedagogical practices.

Figure 1 .

 1 Figure 1. Main page of the platform displaying a list of resources.

Figure 2 .

 2 Figure 2. I2Geo online quality questionnaire.

Intergeo was co-funded by the European Union within the eContentPlus programme,

[2007][2008][2009][2010]. See http://i2geo.net


R Flamary
email: remi.flamary@unice.fr

C Aime

Optimization of starshades: focal plane versus pupil plane

Keywords: Astronomical instrumentation, methods and techniques -Techniques, high angular resolution -Methods, analytical -Methods, numerical

Aims. We search for the best possible transmission for an external occulter coronagraph that is dedicated to the direct observation of terrestrial exoplanets. We show that better observation conditions are obtained when the flux in the focal plane is minimized in the zone in which the exoplanet is observed, instead of the total flux received by the telescope. Methods. We describe the transmission of the occulter as a sum of basis functions. For each element of the basis, we numerically computed the Fresnel diffraction at the aperture of the telescope and the complex amplitude at its focus. The basis functions are circular disks that are linearly apodized over a few centimeters (truncated cones). We complemented the numerical calculation of the Fresnel diffraction for these functions by a comparison with pure circular discs (cylinder) for which an analytical expression, based on a decomposition in Lommel series, is available. The technique of deriving the optimal transmission for a given spectral bandwidth is a classical regularized quadratic minimization of intensities, but linear optimizations can be used as well.

Results. Minimizing the integrated intensity on the aperture of the telescope or for selected regions of the focal plane leads to slightly different transmissions for the occulter. For the focal plane optimization, the resulting residual intensity is concentrated behind the geometrical image of the occulter, in a blind region for the observation of an exoplanet, and the level of background residual starlight becomes very low outside this image. Finally, we provide a tolerance analysis for the alignment of the occulter to the telescope which also favors the focal plane optimization.This means that telescope offsets of a few decimeters do not strongly reduce the efficiency of the occulter.

Introduction

Almost twenty years ago Mayor & Queloz (1995) made the first detection of an exoplanet, and more than a thousand are now recognized 1 . External occulters will be formidable instruments for performing detailed spectral analysis of exo-Earths that have previously been detected by other methods.

The first use of an external occulter dates back to Evans (1948), who used an occulting disk set in front of a very small coronagraph of Lyot (1939). Very soon solar astronomers have identified diffraction effects of a circular occulter and considered various kinds of apodized, shaped, or multiple occulters for space missions, as well described in the review paper of [START_REF] Koutchmy | [END_REF]. Because of implementation difficulties, toothed or multiple discs are preferred over radially graded ones, which have been described in Purcell & Koomen (1962) and Newkirk Jr & Bohlin (1965). Next to the successful Solar and Heliospheric Observatory mission [START_REF] Brueckner | [END_REF]), future solar coronagraphy envisages formation-flying spacecrafts with an occulter of diameter 1.5 m at about 150 m in front of a small telescope (Vives et al. (2009)).

The parameters are completely different for exoplanets,. In rounded numbers, the occulting angle is about 0.2 arcsec against 2000 arcsec for the Sun, and a large 4m diameter telescope is needed with an occulter of 50 m diameter at a distance of 80 000 km. The positive aspect of exoplanet experiments is that diffrac-1 see for example http://exoplanet.eu/ tion phenomena are less difficult to simulate numerically than in the solar case. This is due to two reasons. There are fewer than 20 Fresnel zones compared with 8000 in the solar case, and the main source of diffracted light is a point source instead of a huge 1/2 • extended object. Nevertheless, the problem remains difficult because of the high expected dynamic range.

Progression of ideas for exoplanets first followed developments for the Sun, with little interaction, but now the theoretical calculations and instrumental developments have resulted in much more developments than in the solar case. In his review of space astronomy [START_REF] Spitzer | [END_REF] mentioned the technique of apodization envisaged by solar astronomers. Elaborated petalshaped occulters can already be found in Marchal (1985), and Copi & Starkman (2000) proposed a circular apodizing transmission on a square mask. Studies of shaped and apodized occulters have been reported in many publications, such as Cash (2006), Arenberg et al. (2007), Vanderbei et al. (2007), [START_REF] Soummer | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF], Cash (2011) and Wasylkiwskyj & Shiri (2011), to cite just a few. Although partially transparent petal shaped occulters are envisaged by Shiri & Wasylkiwskyj (2013), the most advanced technological developments are for petaled star shades. The two-step procedure leading to an optimal shaped occulter is well described in [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF]. These authors explained that they first seek for the optimal variable transmission of an apodized occulter and then chose a sufficient number of petals for the shaped occulter to best approach the theoretical result given by the variable transmission. We here focus only on the A&A proofs: manuscript no. Starshade first part of this approach, that is, on the search for an ideal transmission, leaving the final design of a shaped occulter to a future work.

Our purpose is to show that the occulter can be optimized advantageously from the telescope focal plane, which will minimize the level of light that is diffracted by the star in a zone of interest of the focal plane for the exoplanet. The resulting illumination on the telescope aperture appears to be apodized from center to edge. The integrated flux on the telescope pupil is no longer minimized. The observation is nevertheless improved because most of the residual flux is trapped behind the geometrical image of the mask, in a blind area for the observation of the exoplanet. Differences in transmission with an occulter that minimizes the flux over the aperture are small, but the gain in the residual background light at the level of the exoplanet range from 1.9 to 50 depending on the spectral bandwidth, which might correspond to a gain in integration time of a factor 3.6 to 2500 for certain observing conditions, in photon counting mode.

The paper is organized as follows: The fundamental relations that describe the complex amplitudes diffracted by the occulter in the aperture and focal planes are given in Sect. 2. The procedure for optimizing the transmission of a circularly symmetric apodized occulter to minimize the flux on the telescope aperture or for selected regions of the focal plane are described in Sect. 3. Sect. 4 contains results on the optimization problem and a discussion. Additional information is given in two appendices based on an analytic approach of the Fresnel diffraction of a pure circular disc and of a linearly apodized disc.

Expressions of complex amplitudes and intensities in the aperture and focal planes

We denote with 2Ω the overall diameter of the occulter. For an occulting mask, it is convenient to write its transmission as t(r) = 1 -f (r), with the attenuation function f (r) constrained by 0 ≤ f (r) ≤ 1 for |r| ≤ Ω, and f (r) = 0 for |r| > Ω. For a wave of unit amplitude, the Fresnel diffraction at the telescope aperture at a distance z from the occulter (see Fig. 1) can be written as

ψ(r) = 1 - τ z (r) iλz Ω 0 2πξ f (ξ)τ z (ξ)J 0 (2π ξr λz )dξ, (1)
where τ z (r) = exp(iπr 2 /λz) is a quadratic phase-term corresponding to a diverging lens, and J 0 (r) is the Bessel function of the first kind. This relation is equivalent to Eq. 4 of Vanderbei et al. (2007), neglecting here the term of propagation of a plane wave. In general, this integral does not admit an analytical solution, except for f (r) = 1, as shown by Aime (2013) and outlined here in Appendix A. Therefore, Eq.1, which appears as the Hankel transform of f (r)τ z (r), must be evaluated numerically. Computations of these transforms are delicate, as discussed by Lemoine (1994) for example, who recommended a sampling of the function based on zeroes of the Bessel function J 0 . In practice, after several tests described in Appendix A and B and below, we directly used the function NIntegrate of Mathematica, which performs very well for this kind of calculation, thanks to recent improvements of numerical integration. This wavefront ψ(r) arrives onto the aperture of the telescope of transmission P(r), and we denote with ϕ(r) = P(r) ψ(r) the complex amplitude of the wave going through the telescope aperture. For the sake of simplicity, we assume in the following a perfectly circular telescope of diameter 2R, centered on the optical axis of the system, so that P(r) = Π(r/R), defining here for convenience that the box distribution Π(r) equals 1 for |r| < 1 and 0 elsewhere. Note that a more realistic telescope with a central obscuration can be easily inserted into the calculations. Using the circular symmetry of the problem, the complex amplitude of the wave in the focal plane can be written as a Hankel transform of ϕ(r) of the form:

φ(r) = τ F (r) iλF R 0 2πξϕ(ξ)J 0 (2π ξr λF )dξ = τ F (r) iλF φ( r λF ), ( 2
)
where F is the focal length of the telescope. Note that here the phase term τ F (r) is compensated for by the lens phase function of the telescope, meaning that a simple 2D Fourier transform or Hankel transform denoted by φ substitutes the more complex Fresnel propagation of Eq. 1.

It is moreover convenient to calibrate the focal plane in terms of angular units θ = r/F on the sky, and the intensity in the focal plane of the telescope can be written

Φ(θ) = 1 λ 2 | φ( θ λ )| 2 , (3)
instead of just the Airy pattern

Φ 0 (θ) = 1 λ 2 | φ0 ( θ λ )| 2 = π 2 R 4 λ 2 J 1 (2πθR/λ) πθR/λ 2 (4)
for the direct observation without external occulter with a perfect telescope of circular aperture of diameter 2R.

Limiting the analysis to only geometrical aspects, the occulter prohibits the planet observation for an angular radius smaller than Ω/z to the star, and the light coming from the planet passes over the occulter for an angle θ 0 (Ω + R)/z, a variant of the geometric inner working angle Ω/z of [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF], taking into account the telescope size. Note that all the wave and intensity functions defined in this section are illustrated in Figure 1.

Procedure used for optimizing f (r) for integrated intensities

The basic idea of an external occulter coronagraph is to avoid letting the stellar light enter the telescope while leaving the light coming from the exoplanet unchanged . In a first approximation, if the angular position of the planet is larger than θ 0 , we can expect that the flux coming from the planet is almost unaffected by the occulter, and a fair measurement of the efficiency of the external occulter can be the normalized residual integrated intensity over the telescope aperture, for example the quantity Γ defined by

Γ = 2π ∆λ λ M λ m R 0 r Ψ(r)drdλ, (5)
where Ψ(r) = |ψ(r)| 2 is the wavelength-dependent intensity and ∆λ is the spectral bandwidth of the experiment, simply taken equal to λ M -λ m here. In this analysis, the same weight is assigned to each wavelength, independently of the brightness of the source, the quantum efficiency of the detector, or the presence of acquisition filters [START_REF] Soummer | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF]). It would be easy to perform a differently weighted average for a more realistic analysis, if necessary. Our numerical computation has been made from λ m = 380 nm to λ M = 750 nm and with R = 2 m.

To find an optimal shape for the external occulter, a natural approach would be to minimize Γ with respect to f . Note that this approach has been used in Wasylkiwskyj & Shiri (2011), but this is not the only way to optimize the shape of an occulter. For instance, Vanderbei et al. (2007) and [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF] proposed to minimize the upper bound c on the real and imaginary part of ψ(r), ∀r ∈ [0, R]. This approach will ensure that the intensity in the aperture plan will be below 2c 2 and also minimize the flux. We emphasize that while we decided here to minimize the flux, our approach can be readily adapted to the minimization of the upper bound instead of the whole flux Γ.

Because of the conservation of energy, the flux is conserved from aperture to focal planes. Nevertheless, we are in fact interested in the level of background light produced by the star in the area where we observe the planet. A quantity representative of the residual flux may be given by the residual light γ over a zone of the focal plane, such as

γ = 2π ∆λ λ M λ m θ 2 θ 1 θ Φ(θ)dθdλ, (6)
where θ 1 and θ 2 define the region of interest of the observation (Fig. 1), corresponding for example to the habitable zone around the star, or the possible excursion of a known planet over a star, with probably θ 1 close to θ 0 . We therefore have γ ≤ Γ, the equality holding in the limit θ 1 = 0 and θ 2 = ∞.

The two measures of residual light Γ and γ are based on different points of view. If one manages to completely shut off the light in the aperture, then the residual light in the observation area will also shut off, but the reverse is false. In other words, minimizing Γ or γ with respect to f will lead to different optimal solutions. We believe that minimizing γ is better when the objective is to observe an exoplanet in a known area of the focal plane.

Decomposition of f (r) on a basis of functions f k (r)

As already presented in Jacquinot & Roizen-Dossier (1964) for the systematic search for apodizing properties, a common way to optimize the shape of a function is to force this function to be a weighted sum of basis functions. We here aim to optimize the weights α k such that

f (r) = K k=1 α k f k (r) and 0 ≤ f (r) ≤ 1. ( 7
)
Jacquinot & Roizen-Dossier (1964) have shown that, if the expansion contains an infinite number K of terms, "the absolute optimal function according to the criterion chosen" is obtained regardless of the bases used. This result was obtained in the different context of Fraunhofer diffraction, but given the linearity of the equations, the result still holds for our purpose. In practice, however, the number K of elements of the bases must be limited and the choice of basis functions becomes important.

1 0 0 f 1 (r) r f 2 (r) f 3 (r) 1 0 0 f(r) r
After various tests, we used trapezoidal functions for the f k (r), equal to 1 for r ≤ Ω m + (k -1)∆ and equal to 0 for r ≥ Ω m + k∆. The functions linearly decrease from 1 to 0 between Ω m + (k -1)∆ and Ω m + k∆. These function, illustrated in Fig. 2, are very similar to the binary disk for small ∆. It is the same for their Fresnel diffraction pattern, as shown in Appendix B, but this small difference is sufficient to lead to much better numerical results. Note that these basis functions lead to a transmission that is piece-wise linear as opposed to the use of binary disk that leads to piece-wise constant transmission function. For the numerical analysis, we set ∆ equal to 5 cm and used K = 300 functions spreading between Ω m = 10 m to Ω M = Ω = 25 m. An example of such functions is given in Fig. 2. The reason for the existence of an innermost opaque section (i.e. t(r) = 0 or f (r) = 1 for r ≤ Ω m ) is of technical origin linked with the spacecraft, as indicated for example by Vanderbei et al. (2007).

Note that expressing function f (r) as a sum of basis function leads to a more general model potentially at the cost of more complex numerical computation. For instance, one might choose to use a more elaborated basis of functions such as prolate functions. This would potentially lead to very good results because those functions are known to provide efficient apodizers, and might also be excellent occulters. Nevertheless, the physical constraints on f (r) as discussed in the following would be more difficult to express. In this particular case, it would be advantageous to compute these constraint on a finite sampling of the function, as in Vanderbei et al. (2007).

A&A proofs: manuscript no. Starshade 3.2. Constraints on f (r) and f k (r)

The final attenuation function f (r) has to follow a number of physical constraints. On the one hand, some of those constraints can be easily enforced through a wise choice of the basis functions. For instance the smallest radius of apodization is set by Ω m , while the largest radius is selected by Ω = Ω m + K∆. On the other hand, the attenuation function f (r) has to be equal to 1 for r < Ω m , and since all f k (r) functions are equal to 1 for r ≤ Ω m , we have to enforce the linear constraint K k=1 α k = 1. This constraint can also be expressed as a linear equality of the form 1 ⊤ α = 1, where 1 is a vector of ones.

Another more complicated constraint is the fact that 0 ≤ f (r) ≤ 1, ∀r. These block constraints, thanks to the simple structure of the f k functions, can be expressed as a linear inequality constraints Aα ≥ b with

A = C -C C =                     1 1 . . . 1 1 0 1 . . . 1 1 0 . . . . . . . . . 1 0 0 . . . 1 1 0 0 . . . 0 1                     , and b =                             0 . . . 0 -1 . . . -1                            
.

Another constraint that has been used in previous works is the smoothness of the transmission function. This is typically enforced using derivatives of the function, but in our case smoothness can be readily enforced using a classical quadratic regularization term α ⊤ α. Indeed, minimizing the euclidean norm of the α vector under the constraint 1 ⊤ α = 1 tends to promote similar values on all α k , leading to a smooth transmission function (linear apodization for an infinite regularization). This kind of regularization is equivalent to the constraint on the curvature of f described in [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF].

An additional constraint that has been enforced in previous literature is to force f to be monotonic, which in this case is a decreasing function w.r.t. the radius r. This constraint, also enforced in [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF], has an extremely simple form in our formulation, in effect it is the positivity constraint α ≥ 0.

To have a grasp of the performance loss due to this additional constraint, we performed an optimization with and without this constraint in the numerical experiments.

Note that in this work, we made the design choice of using linear apodization functions as basis functions for the transmission. These choices lead to the above-mentioned physical constraints on the weights α. An equivalent formulation has been proposed in Vanderbei et al. (2007) and [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF], where the basis function are trapezoidal functions overlapping such that the final transmission is also piece-wise linear. Copi & Starkman (2000), Cash (2011) and Wasylkiwskyj & Shiri (2011), on the other hand, used offset high-order polynomial or hyper-Gaussian functions.

Below, we present the optimization problems derived from minimizing the intensity Γ in the telescope aperture and minimizing the intensity γ in a region of the focal plane. The optimal functions are denoted as f Γ (r) and f γ (r) for the optimal attenuation w.r.t. the aperture plane and w.r.t. a selected region of the focal plane, respectively.

Minimizing the integrated starlight in the telescope

aperture plane : f Γ (r)

Substituting the decomposition of f (r) of Eq.7 into Eq.1, we obtain

ψ(r) = K k=1 α k 1 -τ z (r) Ω 0 f k (r)τ z (ξ)J 0 (2π ξr λz )dξ = K k=1 α k ψ k (r), (8)
and the corresponding intensity becomes

Ψ(r) = |ψ(r)| 2 = K k=1 α k ψ k (r) 2 = K,K k,l=1 α k α l ψ k (r)ψ l (r) † = α ⊤ K ψ (r)α, (9)
with

K ψ (r) = ψψ † with ψ ⊤ = [ψ 1 (r), . . . , ψ K (r)] a rank one ma- trix of general coefficient K l,k (r) = ψ k (r)ψ * l (r)
, where ψ k (r) is the Fresnel diffraction of a wave of unit amplitude for the basis function occulter f k (r). The flux Γ can then be expressed as

Γ = 2π ∆λ λ M λ m R 0 r α ⊤ K ψ (r)αdrdλ = α ⊤ K Γ α, (10)
where

K Γ = 2π ∆λ λ M λ m R 0 r K ψ (r)
drdλ is a K×K matrix integrating all the surface of the aperture and all the wavelength in (λ m , λ M ). The optimization problem is finally defined as

min α α ⊤ (K Γ + µI)α (11) s.t. 1 ⊤ α = 1 and Aα ≥ b,
where I is the identity matrix and µ is a positive regularization parameter. As discussed in Wasylkiwskyj & Shiri (2011), regularization is important in this case because the matrix K Γ is poorly conditioned and would lead to difficulties when solving the problem. Moreover, as discussed in the previous section, the µ parameter weight the quadratic regularization on the α vector, promoting smoothness in the transmission function. This smoothness has also been shown to be important in practice in [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF], where an equivalent smoothness constraint was inserted in the optimization problem. Also note that in practice one can easily express the smoothness constraint of [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF] in our optimization problem instead of using a regularization term. Since with the trapezoidal basis, the secondorder derivatives are exactly computed using the finite difference between adjacent α components, the corresponding constraints can be obtained by adding lines to the A matrix that contain finite difference operators, and lines to b that contain the constraint parameter σ of [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF].

The optimization problem defined Eq. ( 11) is a constrained quadratic optimization [START_REF] Boyd | Convex optimization[END_REF]. The linear constraints forbid the use of the closed-form solution as proposed in Wasylkiwskyj & Shiri (2011), but there exist several efficient optimization strategies in the literature, among which the active set approach presented in Vanderbei & Shanno (1999).

Minimizing the starlight in a region of the telescope focal plane: f γ (r)

To express, the intensity in the focal plane, we used the same decomposition of the f (r) as a sum of functions f k (r)

as above, the complex amplitude in the focal plane can be written as a sum of elementary function φ k (r), Fourier transforms of P(r)ψ k (r) according to Eqs 7 and 2. The intensity in the focal plane of Eq. 3 becomes

Φ(θ) = 1 λ 2 K k=1 α k φ k ( θ λ ) 2 = α ⊤ K φ (θ)α, (12)
and the intergrated density

γ = 2π ∆λ λ M λ m θ 2 θ 1 θ α ⊤ K φ (θ)αdθdλ = α ⊤ K γ α. (13)
The resulting optimization problem is exactly the same that of Eq. 11, but with a different metric K γ .

The optimal functions f Γ (r) and f γ (r) are different because the solutions α are different.Because of conservation of energy between aperture and focal planes, these functions become identical in the special case θ 1 = 0, θ 2 = ∞, where K Γ = K γ . Our optimization in the focal plane was performed for θ 1 = 0.1 arcsec and θ 2 = 0.5 arcsec.

Results and discussion

Numerical computations

In the following, we denote Ψ Γ (r) and Φ Γ (θ) the intensities obtained in the aperture and focal plane of the telescope using an external occulter of transmission 1 -f Γ (r) minimizing the integrated intensity Γ and described in Sect. 3.3. Similarly, we denote Ψ γ (r) and Φ γ (θ) the results obtained using an occulter transmission 1 -f γ (r), as described in Sect. 3.4.

Note that while obtaining the matrices K Γ and K γ is a computationally intensive problem, we can use heavily parallel computing to estimate as a first step the ψ k (r) and φ k (θ) functions through numerical integration. This was performed using Mathematica, with K computational jobs by ≈ 100 parallel processes on eight core Intel processors. In the numerical experiments, both ψ k (r) and φ k (θ) were regularly sampled on 3001 samples from respectively r ∈ (0, R+1) and θ ∈ (0, 0.644) arcsec. Note that ψ k (r) was computed outside of the telescope so that we can used it for the tolerance analysis in section 4.4. Moreover, the behavior of the different optimization strategies outside of the telescope aperture is also interesting. The wavelength λ has also been sampled regularly with 21 samples in the visible light interval (380, 750) nm. After computing ψ k (r) and φ k (θ), a numeric 2D integration can be performed to obtain the matrices K Γ and K γ . The regularization parameter µ in the optimization problem was chosen adaptively with µ = µ 0 max k,l |K k,l | to be less dependent on the metric matrix K. µ 0 was chosen by hand to allow good attenuation performance with a limited dynamic in the solution function f . For the monochromatic study, we set µ 0 = 10 -8 and for the more complex study on a wide bandwidth, less regularization was necessary with µ 0 = 10 -10 because of a more complex problem and a better conditioning of the K matrix.

Finally, the code for all the numerical experiments will be freely available on the authors website to promote reproducible research 2 .

Optimization for monochromatic light

For a single monochromatic wave, here λ = 562 nm, we have represented the apodizing curves f Γ (r) and f γ (r) in Fig. 3 and the resulting intensities Ψ Γ (r) and Ψ γ (r) (top curves), and Φ Γ (θ) and Φ γ (θ) (bottom curves) in Fig. 4. We emphasize that these curves were obtained assuming a perfect circular telescope of radius R. The f -curves are only slightly different and are really similar to a linear apodization. Moreover, in the monochromatic case, enforcing a monotonic f leads to the same results for the thickness of the trace, which will not be the case anymore for a wide spectral bandwidth, as we show below.

The results for the Ψ and Φ-curves are quite different. The residual starlight Ψ Γ (r) is very well attenuated across the telescope pupil, as found by Vanderbei et al. (2007). In contrast, Ψ γ (r) presents a center-to-limb decrease which suggests an apodized structure. It is interesting to note the behavior of the curves beyond the region of optimization. There, the relative position of the curves is reversed very rapidly, the focal plane optimization becoming lower than the pupil plane optimization. This is an aspect that deserves further studies.

The intensity in the focal plane is also very interesting. While f Γ tends to minimize the intensity at all angles θ, f γ only focuses on the region of interest, and a strong residual light is observed in the blind area behind the occulter.

The resulting flux Γ, γ for all the different optimizations is given in the upper part of Table 1. Interestingly, f Γ leads to a flux Γ in the aperture plane 10 5 times fainter than f γ . Nevertheless, in the focal plane the flux γ in the observation area is ≈ 50 times fainter when using f γ . This clearly shows the importance of using focal plane optimization for optimal apodization.

Optimization for a wide spectral bandwidth

Next we investigated the two optimization strategies on a wide bandwidth. For this we compute a regular sampling in the interval λ ∈ (380,750) nm and computing the average flux across λ. Note that both a sampling of 20 and 100 values were investigated, and while we used in the experiments a sampling of 100 for a better precision, the final results were very similar. In prac- Angle θ in arcsec

Φ(r) Focal plane Φ Γ (θ) Φ γ (θ) Φ Γ (θ) ≥ 0 Φ γ (θ) ≥ 0
Airy pattern Apod. area Fig. 4. Top curves: intensity in the telescope aperture plane for an optimization that minimizes the integrated intensity Γ in the aperture plane (blue curve) or γ for optimizing the observation between 0.1 and 0.5 arcsec in the focal plane (red curve), for λ = 562 nm. Note the apodization profile for the γ curve and the inversion of relative position of the curves outside the region of optimization. Bottom curves: focal plane intensity obtained for the two optimizations. An airy pattern corresponding to a 10 -12 fainter exoplanet in the observation area is reported in black. The apodization area defined by the smallest and largest radius of the occulter is also reported in green. Note the inversion of behavior of curves with regard to the overall occulter profile.

tice the variation of the wavefront is extremely smooth w.r.t. λ and the selected sampling in sufficient.

For multiple wavelengths, λ ∈ (380,750) nm, the results are very different. First the optimal functions are much more similar, as illustrated in Figure 5. The functions also tend to oscillate and the monotonic variants are this time different from the free variant of the optimal functions.

The resulting intensities in the aperture and focal plan are plotted in Figure 6. In the aperture plane, the intensity is still apodized for f Γ , but the flux is much more attenuated than in the monochromatic case. In the focal plane, the intensities retain the same tendencies, but the gain due to the focal plane optimization is lower. Interestingly, enforcing a monotonic function through α ≥ 0 leads to a clear loss in terms of attenuation for both f Γ and f γ .

In terms of quantitative performances, Table 1 (lower part) shows that the flux in the observation area of the focal plane is still more attenuated with f γ with a gain of 1.9 in the chromatic case (1.6 with positivity constraint). While the gain in performance is smaller than in the monochromatic case, it is far from negligible when observing faint objects around a star. Interestingly, in the chromatic case, the α ≥ 0 constraint leads to a loss of 2 in attenuation (for both f Γ and f γ ), proving that relaxing the constraint can lead to additional gain as long as the apodization is physically realizable.

An additional visualization comparing the performances of f Γ versus f γ in the focal plane is given in Fig. 7. The intensity Φ γ (θ) is plotted as a function of Φ Γ (θ) for θ ∈ (0, 0.7) arcsec. The red line allows a comparison between the methods; if a curve is below this line, it means that Φ γ (θ) < Φ Γ (θ) and vice versa. We clearly see that Φ γ (θ) > Φ Γ (θ) for θ < θ 1 in the blind area behind the occulter and that Φ γ (θ) < Φ Γ (θ) ∀θ ∈ (θ 1 , θ 2 ) in the

Tolerance analysis to positioning errors of the telescope

A study of the tolerance sensitivity to positioning errors for a 4m circular telescope out of its on-axis position was carried out up to an offset of 1 m.

Because of this offset, the resulting pupil plane and focal plane intensities are no longer circularly symmetric, and the computation of the final focal plane image can no longer use the Hankel transformation of Eq.2. The focal plane image must Angle θ in arcsec

Φ(r) Focal plane Φ Γ (θ) Φ γ (θ) Φ Γ (θ) ≥ 0 Φ γ (θ) ≥ 0
Airy pattern Apod. area Fig. 6. Same representation as in Fig. 4, but for an optimization made for the whole spectral bandwidth between 380 nm and 750 nm. The behavior of the inversion of the curves is similar to the monochromatic case, although it is somewhat reduced by the effect of the chromatism.

Table 1. Resulting light flux in the aperture (Γ) of focal plane (γ) for all the optimization schemes. The upper part of the table corresponds to an optimization on a unique wavelength λ, whereas the lower part of the table corresponds to an optimization on a large spectrum.

λ Flux f Γ (r) f Γ (r) ≥ 0 f γ (r)
f γ (r) ≥ 0 Flux ratio Flux ratio ≥ 0 λ=562 nm Γ 9.07e-13 9.55e-13 5.06e-08 5.52e-08 1.79e-05 1.73e-05 γ 7.66e-13 8.06e-13 1.38e-14 1.44e-14 5.57e+01 5.61e+01

λ ∈(380,750) nm Γ 2.91e-13 5.88e-13 2.63e-12 7.01e-12 1.11e-01 8.38e-02 γ 7.08e-14 1.42e-13 3.74e-14 8.99e-14 1.89e+00 1.58e+00 be computed taking the modulus squared of the two-dimensional Fourier transform of the complex amplitude of the wave on the telescope aperture.

In practice, the study was made using a discrete Fast Fourier Transform of an array of 1024 × 1024 points corresponding to an overall zone of 20 × 20 meters, inside which the telescope of diameter 2m was defined by 1 or 0 transmitting pixels over about 205 points in diameter. As a result, altogether 33 081 points were set equal to 1. This percentage of zero-padding is enough to obtain a satisfactorily sampled focal plane image.

Positioning errors were simulated by moving the complex amplitude obtained in the shadow of the occulter over the telescope aperture in one direction, by steps of 1 pixel or about 2 cm. The two-dimensional array was filled using an interpolation of the one-dimensional complex amplitude of the wavefront computed with Eq. 1.

The focal plane image was obtained taking the modulus squared of the Fourier transform of the wave on the telescope aperture. The effect of wavelength was taken into account afterward. The diffraction pattern increases in size with the wave-length. For that, a two-dimensional interpolation of the result was used to resample the image. Moreover, a scaling factor in intensity inversely proportional to the square of the wavelength was applied, as in Eq.4. As a result, the integrated flux in the focal plane becomes wavelength independent and equal to that crossing through the telescope aperture. As before, for the sake of simplicity, we assumed that the received flux is constant over the spectral bandwidth.

Figures 8 and9 present results obtained for transverse displacements of the telescope up to 1 m from its nominal position. We recall that the occulter shapes optimize the measurement when the telescope is at its nominal on-axis position, either for the integrated intensity on the telescope aperture, or for a region in the focal plane between 0.1 and 0.5 arcsec. In both cases the computation was made for the whole bandwidth of λ ∈ (380,750) nm.

Three-dimensional representations of the normalized intensities in the telescope aperture were used to enhance differences of intensities in the results of the two optimizations, and the apodization-like pattern obtained for the focal plane optimiza- tion. Both curves were clipped to make the central parts of the images well visible. As expected, there is much less light for the aperture plane minimization when the telescope is at its nominal position. Since the optimization was calculated for a 4m diameter telescope, there is a rise in the intensity at the edges for this large offset. This increases becomes somewhat surprisingly in favors of the focal plane optimization for a large off-set, as already discussed in Fig. 6.

In Fig. 8 focal plane intensities are also represented in color levels, with a common color scale for the telescope on-axis and off-axis of 50 cm and 1m. The two circles of radius θ 1 = 0.1 arcsec and θ 2 = 0.5 arcsec define the angular area inside which the residual intensity is integrated. The focal plane optimization is found to always be better than the aperture plane optimization. From comparing these curves, it is clear that the residual diffracted light remains concentrated longer inside the first circle for the focal plane optimization, which is especially clear for the 50 cm offset images. This behavior is confirmed by Fig. 9 which represents the integrated intensity in this region with a logarithmic scale.

Conclusion

We addressed the problem of optimizing of the shape of a circularly symmetric apodized occulter. To this end, we proposed a general expression for the attenuation function, as a weighted sum of basis functions. We expressed the optimization problem as the minimization of the flux, first in the pupil plane, as is commonly done in the literature, and then in the focal plane of a 4.m diameter telescope, for a region between 0.1 to 0.5 arcsec in which the exoplanet is supposed to be observed.

Numerical experiments show that the focal plane optimization leads to a better attenuation of the starlight than the pupil plane optimization. One interesting result is that while the focal plane optimization leads to a strong flux in the pupil plane, most of this light is concentrated in the area behind the occulter in the focal plane, which is not a problem for a perfect telescope. Finally, the robustness of our approach to lateral misalignment of the occulter was investigated and showed that the focal plane optimization also yields a better robustness.

As already indicated, the optimization study was made for a 4m telescope assumed to be at its on-axis position. A simple procedure to make the experience less sensitive to defect position would be to perform the optimization for a larger aperture than is actually used. A better procedure would be to estimate the statistics of pointing errors and include them in the optimization procedure. Then a better tolerance to pointing errors would be obtained, probably at the expenses of a tolerable loss in the starlight rejection. Nevertheless, our analysis showed that telescope offsets of a few decimeters will not strongly reduce the efficiency of the occulter. This conclusion was obtained for the condition that the telescope has a perfectly circular aperture and another study is required to model a telescope with central obscuration. Moreover, it would be interesting to study the effect of a a more realistic telescope with small defaults of phase and amplitude. Deviations from this perfection could lead to somewhat different results, but the exigency of quality is limited to a maximum of 0.5 arcsec, or about twenty Airy rings for a 4m diameter telescope. Finally, with the knowledge of the exact PSF of the telescope, one can readily adapt our optimization approach. A comparison of pupil vs focal optimization for the maximum-magnitude minimization as in [START_REF] Kasdin | Proceedings of the 2013 IEEE Aerospace Conference[END_REF] would also be interesting.

Fig. 1 .

 1 Fig. 1. Schematic diagram of an external occulter coronagraph with the principal notations used in the paper. The observation, occulter, aperture, and focal telescope planes are shown.

Fig. 2 .

 2 Fig. 2. Left: illustration of the basis function used to represent the attenuation function f (r) with an example of the combination of three of these functions (right). For clarity of representation, the values of Ω m and ∆ are not to scale of those used for the f k (r) basis functions.

2Fig. 3 .

 3 Fig. 3. Radial cut of the attenuation of the functions f (r) for aperture and focal optimization at a wavelength of 562 nm. The dashed curves denoted as ≥ 0 correspond to an optimization with α ≥ 0 constraint.

Fig. 5 .

 5 Fig. 5. Radial cut of the attenuation of the functions f (r) for aperture and focal optimization for the whole spectral bandwidth between 380 nm and 750 nm.

  Fig. 7. Parametric plot of the observed intensities in the focal plane: xaxis, aperture plane optimization, y-axis, focal plane optimization. The curve is drawn for the wide spectral bandwidth without positivity constraints. From top-right to bottom-left the curve follows the increase of θ, clearly showing the concentration of light behind the occulter for the focal plane optimization. The line y = x clearly shows that the γ-curves always give better results than the Γ-curves, especially in the observation area (reported in green).

Fig. 8 .

 8 Fig. 8. Top curves: pupil plane intensities for an offset of 1 m for λ ∈ (380,750) nm. Density curves below represent the focal plane intensities (log scale) for telescope positions: on-axis (upper curves) and offsets of 50 cm (intermediate curves) and 1 m (lower curves). Left: aperture optimization, right: focal plane optimization.

Fig. 9 .

 9 Fig. 9. Average residual light level γ for the focal plane (a) and pupil plane (b) optimization, as a function of the telescope offset in centimeters for λ ∈ (380,750) nm. The average is computed over the region from 0.1 to 0.5 arcsec delimited in Fig. 8 by two circles.

  Fig. B.1. Differences of moduli of Fresnel diffractions for raw and linearly apodized functions for Ω values of 10 m (thick curve) and 10 m (thin curve).

Acknowledgements. We want to thank the referee Jeremy Kasdin for his very constructive comments, and in particular for his suggestion to study the effect of telescope misalignment. This work was granted access to the HPC and visualization resources of Centre de Calcul Interactif hosted by Université Nice Sophia Antipolis (http://calculs.unice.fr/en).

Appendix A: Quality test of the numerical computation using an analytical expression for a sharp-edged circular occulter Aime (2013) has shown, using an approach similar to [START_REF] Born | Principles of optics: electromagnetic theory of propagation, interference and diffraction of light[END_REF], that an analytic expression using Lommel series could be derived for the Fresnel diffraction of a sharp circular occulter of the form f (r) = f • (r, Ω) = Π(r/Ω), defining here Π(r) as the function equal to 1 for |r| ≤ 1 and 0 otherwise. At the wavelength λ, the complex amplitude ψ • (r, Ω) of the wave at the distance z for (i) inside the geometrical umbra (r < Ω) and (ii) outside (r > Ω), can be written using two series:

At r = Ω, both series converge to the same value given by

These series are alternative series for the real and imaginary terms, and according to the Leibniz estimate, an upper bound of the error for the sum limited to n = N terms is given by the absolute value of the N + 1 term. The convergence is easy when r/Ω is smaller than 1, for the first sum for r < Ω inside the geometrical umbra, or the inverse outside. The main error occurs near the transition zone r = Ω. In Fig.

A.1 we compare the analytical and numerical results for an example of a circular occulter of 50 m diameter set at 80 000 km, and for a wavelength of 550 nm. The Lommel series were computed for 100 terms, and there is no difference between the two methods. Results were drawn for the real, imaginary, modulus, and unwrapped phase of the Fresnel diffraction. Note that this is just an example, in particular, the real and imaginary part of the diffraction pattern are highly sensitive to a small variation of any of the parameters. The modulus of the intensity pattern clearly shows the visible Poisson spot at the center of the pattern. The unwrapped phase presents a strong phase variation in the center of the umbra and a quieter structure outside the geometric umbra, that is for a region utilized by the telescope for the planet observation. The consistency between the two analytic and numerical results is recognized here as a proof of the quality of the numerical calculation.

Appendix B: Numerical calculation of linearly apodized basis functions f k (r)

As indicated in the body of the paper, we tried several forms for the set of apodized basis functions f k (r), and decided to use the simple trapezoidal functions described in the body of the paper because they affect the final solution the least. We did not succeed in deriving an analytic form for the Fresnel diffraction of these functions, but comforted by the results obtained on the Fresnel diffraction of Π(r/Ω), we used the direct numerical calculation of Mathematica.

We recall that these functions are linearly apodized disks of radii varying from 10 m to 25 m in steps of 5 cm, the apodization