domain: spi math sde

text:
R S Mouradi

Y Audouin

C Goeury

N Claude

P Tassi

K El Kadi Abderrezzak

Sensitivity analysis and uncertainty quantification in 2D morphodynamic models using a newly implemented API for TELEMAC2D/SISYPHE

Keywords: Uncertainty quantification, Morphodynamic modelling, API, Monte-Carlo like algorithms, Sensitivity analysis

In this study, the modules TELEMAC-2D and SISYPHE of the Telemac-Mascaret Modelling System (TMS) have been used in combination with the OpenTURNS library (www.openturns.org) to perform an uncertainty quantification analysis of two-dimensional morphodynamic problems. Open-TURNS is a scientific library usable as a Python module dedicated to uncertainties treatment.

A recently implemented API (Application Program Interface) allowed the communication between OpenTURNS and TELEMAC-2D/SISYPHE, and therefore the efficient implementation of Monte-Carlo like algorithms. In this problem, each uncertain sedimentological parameter, e.g. inlet mean diameter, Shields parameter, etc. has been associated to a statistical distribution, defined with OpenTURNS. A number of TELEMAC-2D/SISYPHE simulations has been proposed regarding the predefined random entries in order to guarantee the convergence of the Monte Carlo-like algorithms. This work allowed the implementation of uncertainty quantification analysis of computationally intensive morphodynamic simulations in the TMS. Thanks to the access to computer resources and optimized software, we were able to perform the uncertainty quantification analysis with a large set of variables, and therefore push the study further with the correlations effects analysis.

I. INTRODUCTION

Morphodynamic simulations have been increasingly used in the last few decades to model the bed evolution in rivers, coasts and estuaries. In this context, most of the equations are empirical and the parameters involved in the calculations are generally poorly defined in literature. The impact of the uncertainties related to those parameters remains unknown.

In order to quantify the impact of inputs uncertainties on simulations results, an uncertainty study is conducted. Ranking the variables in terms of influence allows to orientate the investigations when performing measurements or calibrating parameters for the simulations. In this study, the uncertainty quantification is applied to SISYPHE [START_REF] Tassi | SISYPHE USER'S MANUAL -Release 6[END_REF], a sediments transport module, coupled with TELEMAC2D [START_REF] Ata | TELEMAC MODELLING SYSTEM -Release 7[END_REF] for hydrody-namics, that integrate the TELEMAC-MASCARET modelling system.

The Monte Carlo method is used to propagate the uncertainties through SISYPHE. This approach requires random generation of several configurations of inputs, using their probability distributions. Successive deterministic model simulations are then submitted to generate a set of responses that correspond to the set of inputs.

In order to have total control over the simulation's parameters and therefore conduct an efficient uncertainty study, an API (Application Program Interface) is developed for SISYPHE. This interface, when coupled with TELEMAC2D's available API, makes running hundreds of cases simultaneously possible through a cluster, taking total benefit from the available processors. Running an optimal number of cases guarantees the statistics convergence.

Finally, The pre-processing of uncertain data, as well as the post-processing of the results, are done using the OpenTURNS uncertainty library [START_REF]Airbus-EDF-IMACS-Phimeca[END_REF]. This paper is organized as follows: a description of the context and general goals, followed by the present study objectives are given. Section III deals with the API's implementation and coupling with TELEMAC2D. Section V discusses the uncertainty quantification steps and alternates theory and results for each of these. In this section, a sensitivity analysis followed by an uncertainty propagation are investigated. Correlations between variables are also studied using copulas and an ANCOVA (ANalysis of COVAriance) sensitivity analysis is applied. In the last section, conclusions and perspectives are drawn.

II. CONTEXT AND GOALS

This study is set out in the context of EDF numerical tools development. EDF's R&D National Laboratory for Hydraulics and Environment department (LNHE) uses the TELEMAC-MASCARET system to simulate complex hydroenvironmental phenomenons (such as dam breaks and flooding risks) in order to anticipate the risks related to electrical production.

TELEMAC-MASCARET results are therefore expected to produce highly reliable results. However, a great number of parameters used in these studies, specifically in morphodynamic simulations, can be set by the user and are uncertain. In order to determine the uncertain parameters impact on the system's result, an uncertainty quantification study is necessary.

In this context, an uncertainty quantification study is conducted in the morphodynamic simulation module SISYPHE.

To make this study possible, an API (Application Program Interface) is implemented and coupled to TELEMAC2D, to guarantee the inter-operability with SISYPHE.

III. API'S DEVELOPMENT

The work here focused on the implementation of SISYPHE's API and its coupling to TELEMAC2D's already available API via FORTRAN modules. The API's main goal is to have control on a simulation while running a case. For example, it must allow the user to stop the simulation at any time step, retrieve some variables values and change them if necessary. SISYPHE's API is contained in the source folder of the TMS, as shown in Figure 1, and can therefore call all SISYPHE's subroutines to conduct morphodynamic simulations.

Fig. 1: Extract of TELEMAC-MASCARET sources folders

In order to make this possible, a FORTRAN structure called instance was developed in the API. It contains a list of variables declared as pointers (memory addresses [START_REF] Delouis | Cours Fortran 95 -IDRIS[END_REF]) that are pointing to SISYPHE's variables. This gives direct access to the physical memory of variables, and allows therefore to retrieve their values, and modify them. Furthermore, modifications have been made in SISYPHE's main subroutines to make morphodynamic cases execution possible time step by time step. Finally, parallel runs have also been treated.

In addition to this, to make running coupled cases via the API possible, a communication interface is developed in FOR-TRAN. This interface contains communication subroutines that send TELEMAC2D's variables to SISYPHE and vice-versa. It also contains subroutines that manage coupled cases and take into consideration the coupling period.

A number of modifications in TELEMAC2D and SISYPHE sources were necessary. These modifications, along with the API and the coupling developments, were validated using three different compilers (NAG, IFORT, GFORTRAN), on classical SISYPHE cases and coupled TELEMAC2D-SISYPHE cases, available in the system.

IV. UNCERTAINTY TREATMENT LIBRARY OPENTURNS

OpenTURNS (Open source initiative to Treat Uncertainties, Risks'N Statistics) [START_REF]Airbus-EDF-IMACS-Phimeca[END_REF] is an open source C++ Library for uncertainty treatment used through python scripts. It is codeveloped since 2005 by EADS IW, EDF R&D and PHIMECA Engineering. Various statistical methods are implemented in this library and allow to follow the uncertainty study steps [START_REF] Pasanisi | Guide méthodologique pour le Traitement des Incertitudes[END_REF] represented in Figure 2.

                 ∂h ∂t + ∇ • (hu) = 0 ∂u ∂t + u ∂u ∂x + v ∂u ∂y = -g ∂Z s ∂x + F x + 1 h ∇ • (hν t ∇u) ∂v ∂t + u ∂v ∂x + v ∂v ∂y = -g ∂Z s ∂y + F y + 1 h ∇ • (hν t ∇v) (1)
Where u and v are the velocity components, h the water depth, Z s the free surface elevation, F x and F y are source terms (for example the friction) and ν t the eddy viscosity.

2) Morphodynamic model: The main goal of a morphodynamic simulation is to predict the bed evolution of a given domain considering the flow conditions. Two types of sediment transport exist [START_REF] Garcia | Sedimentation engineering : processes, measurments, modeling, and practice ASCE[END_REF]:

• Bed-load : sediments are transported without losing contact with the bed. Their speed is therefore lower than the flow velocity.

• Suspension : corresponds to sediment transport in the absence of direct contact with the bottom. Sediments are here transported with a speed equal to the velocity of the flow.

In this study, only bed-load uncertain parameters are investigated. The main equation for this phenomenon is Exner's equation 2, that calculates the evolution of the bottom elevation

Z f . (1 -λ) ∂Z f ∂t + ∇ • - → q b = 0 ( 2
)
Where λ is the porosity of the sediment, ∇• the divergence operator and -→ q b = q b × (cos α t -→ e x + sin α t -→ e y ), α t being the angle between the flow direction and the sediment transport direction and q b the bed-load transport rate per width.

The transport rate q b can be calculated using several empirical formulas. In this study, we focus on Meyer-Peter and Müller formula [START_REF] Meyer-Peter | Formulas for Bed-Load Transport[END_REF], as written below:

q b g(s -1)d 3 = α M P M (θ -θ c ) 3/2 (3)
• α M P M being a transport coefficient specified by the user;

• d the sediments mean diameter;

• s = ρ s /ρ f the relative density, where ρ s is the sediments density and ρ f the fluid density;

• θ and θ c resp. the Shields number and its critical value indicating the movement's threshold.

The Shields number θ is the dimensionless shear stress τ calculated as in equation 4. The threshold θ c is given as a parameter of the simulation by the user.

θ = τ g(ρ s -ρ)d (4)
3) Summary of the uncertain parameters: The porosity λ, the mean diameter d, the Shields critical parameter θ c and five other parameters, described in [START_REF] Tassi | SISYPHE USER'S MANUAL -Release 6[END_REF], are considered as uncertain for the current study. A summary of these variables and their definitions are given in Table I Section V-B.

4) Study cases:

The uncertainty quantification is conducted for two distinct cases.

Experimental case -Channel : The channel case is a numerical reproduction of Ashida and Michiue's experiment [START_REF] Wu | Depth-averaged two-dimensional numerical modeling of unsteady flow and uniform sediment transport in open channels[END_REF] in which the erosion downstream of a dam is studied. The case, illustrated in Figure 3, is modelled in 2D with boundary conditions on flow discharge (Red boundary in Figure 3 -Q = 0.0314 m 2 • s -1 ) and imposed water depth (Blue boundary in Figure 3 -H = 0.06 m), as well as walls on the remaining boundaries. The channel has a 20 m length, 0.8 m width and 0.2 m slope. An experimental duration of 10 hours is simulated using a time step ∆t = 0.1 s within 1 minute of sequential TELEMAC2D/SISYPHE calculations on an Intel-Xeon(R) 3.40GHz core.

B. Uncertainty sources quantification

The goal here is to define variation intervals and probability density functions (PDF) for each uncertain parameters. We remind that, for a random variable X defined on an interval [a, b], a PDF f (x) is defined as follows :

P(X ∈ [a, b]) = b a f (x)dx (5)
For all the variables from 2 to 8 defined in Table I, variation intervals were found in literature, without further information about their probabilities. Consequently, using the principle of maximum entropy [START_REF] Sudret | Contributions to structural reliability and stochastic spectral methods[END_REF], uniform probability density functions are chosen for those variables.

For the experimental channel case, only one mean diameter measure is available and is subject to measurements errors. The errors interval is considered as a variation interval and a uniform PDF is applied for the mean diameter on its measurement interval. For the real bifurcation case, several samples of sediments are extracted in different sections of the river. Different values of the mean diameter are possible. In order to take all of them into consideration, a PDF that corresponds to the sample is selected and validated via the QQ-Plot method [START_REF]Airbus-EDF-IMACS-Phimeca[END_REF].

The summarized variables and their probability density functions are presented in Table I. Mean :

E[Y ] = µ Y = 1 N N i=1 M (x i ) (6)
Variance : V ar(Y ) = 1 N -1 N i=1 [M (x i ) -µ Y ] 2 (7)
Standard deviation :

σ Y = V ar(Y ) (8)
These statistical moments are useful for both the uncertainty sensitivity analysis and uncertainty propagation.

The convergence order of the Monte-Carlo sampling method is given by the Central Limit Theorem [START_REF] Goeury | Analyse de sensibilité et propagation d'incertitude dans les modèles numériques : Application à l'hydraulique[END_REF] as

O 1 √ N .
The estimated statistics are also random quantities and are impacted by estimation uncertainties. Confidence intervals on estimators should therefore be calculated. The non-parametric "bootstrap" method provides information about the statistics uncertainties given few hypothesis [START_REF] Kadiri-Ottmani | Traitement statistique d'un échantillon[END_REF]. Let x = (x 1 , ..., x N ) denote a sample of N independent and identically distributed realizations according to a probability density function f (x). The statistical moment θ = T (F ) (mean, variance, etc.), is estimated by θ = T ( F ), where F is the empirical cumulative density function that gives equal probability 1 N to each observed value x i defined by :

F (x) = 1 d N i=1 1 xi≤x (9)
The idea of the non-parametric bootstrap is to simulate data from the empirical cumulative density function. Given that F is build upon equal probability for the observations (x 1 , ..., x N ), a sample of same size N from F would simply be a selection from (x 1 , ..., x N ) with repeated values. A number of B samples are generated following this strategy, and estimators properties can therefore be deduced as shown in Figure 5. The confidence interval I θ is then estimated as follows :

I θ = [0.025 -q θ b , 0.975 -q θ b ] (10)
Where α -q X is the α-quantile of a variable X defined as :

P(X ≤ q X (α)) = α ∀α ∈ [0, 1] (11)
2) Analysis of variance: The main goal of a sensitivity analysis is to rank the uncertain parameters according to their influence. In order to do so, a definition of ranking indices is necessary. The indices used here are called Sobol Indices [START_REF] Ioss | Revue sur l'analyse de sensibilité globale de modèles numériques[END_REF].

The definition of Sobol Indices is a result of the ANOVA (Analysis Of VAriance) variance decomposition. In fact, given a set of V independent uncertain parameters X = (X 1 , ..., X V ), the variance of a response Y = M (X) can be calculated, using the total variance theorem [START_REF] Caniou | Global sensitivity analysis for nested and multiscale modelling[END_REF], as follows [START_REF] Ioss | Revue sur l'analyse de sensibilité globale de modèles numériques[END_REF]:

V ar[Y ] = V i=1 V i (Y ) + i<j V ij (Y ) + ... + V 12..V (Y ) (12) Where V i (Y ) = V ar[E(Y |X i )] and V ij (Y ) = V ar[E(Y |X i X j )] -V i (Y ) -V j (Y ), etc. E(Y |X i ) is Y 's conditional expectation with the condition that X i remains constant.
One can therefore calculate first order sensitivity indices that estimate the influence of a variable X i without its interactions with other variables:

S i = V i (Y ) V ar[Y ] = V ar(E[Y |X i ]) V ar[Y ] (13)
And total indices that estimate the global influence of a variable (including interactions):

S T i = S i + j =i S ij + j =i,k =i,j<k S ijk +... = 1- V -i (Y ) V ar[Y ] (14)
V -i (Y ) being conditional expectations variances that do not involve X i .

There are two ways to estimate the Sobol indices defined in equations 13 and 14 :

SALTELLI method [7]:

Step 1: Two independent samples A and B of size N are generated for the V uncertain variables. For example, sample A can be written as follows:

A =       x A,1 1 x A,1 2 ... x A,1 V x A,2 1 x A,2 2 ... x A,2 V . . . . . . . . . . . . x A,N 1 x A,N 2 ... x A,N V       (15)
A new sample C is created using columns of B except for column i that is replaced with data from A:

C =       x B,1 1 ... x A,1 i ... x B,1 V x B,2 1 ... x A,2 i ... x B,2 V . . . . . . . . . x B,N 1 ... x A,N i ... x B,N V       (16)
Simulations using the samples A, B and C result with :

     y A k = M ((x A,k 1 , ..., x A,k V )) k = {1, N } y B k = M ((x B,k 1 , ..., x B,k V )) k = {1, N } y C k = M ((x C,k 1 , ..., x C,k V )) k = {1, N } (17)
Which can be used to estimate Sobol indices as follows:

S i = 1 N N k=1 y A k y C k -(µ Y A ) 2 σ 2 Y A (18) S T i = 1 - 1 N N k=1 y B k y C k -(µ Y B ) 2 σ 2 Y B (19)
Overall, for a given sample size N , (V + 2) N simulations are necessary to estimate the first order and total Sobol indices for each variable X i .

Polynomial chaos method (PCE):

The models response can be approached by an analytical function :

M (X) = M 0 + V i=1 M i (X i ) + 1≤i<j≤V M i,j (X i , X j ) + ... + M 1,...,V (X 1 , ..., X V ) (20)
M i (X i ) represents the variable X i 's contribution to the result of the simulation. The variance can therefore be calculated as follows [START_REF] Caniou | Global sensitivity analysis for nested and multiscale modelling[END_REF]:

V ar[Y ] = u⊆{1,...,V } [V ar[M u (X u )]+ v⊆{1,...,V },v∩u=∅ Cov [M v (X v ), M u (X u )]] (21)
In this particular case of independent variables, the covariance term Cov [M v (X v ), M u (X u )] vanishes, and the decomposition becomes then equal to ANOVA. Sobol indices can therefore be estimated as:

S i = V ar [M i (X i )] V ar[Y ] (22)
S T i = u⊆{1,...,V },i∈u V ar [M u (X u )] V ar[Y ] (23)
The contributions M i (X i ) can be calculated by estimating the models response using a polynomial chaos expansion, which can, in a simplified way, be written as:

M (X) = |α|≤P a α Ψ α (X) (24)
Where {Ψ α , α ∈ N V } is a multivariate polynomial basis and a α adequate coefficients for the estimation of the model's response, that can be determined using projection methods [START_REF] Caniou | Global sensitivity analysis for nested and multiscale modelling[END_REF].

The X i -univariate polynomials shares are the exact contribution of X i to the polynomial expansion, and are therefore an estimation of M i (X i ).

3) Results:

Channel's results: For a sample of size N = 30000, a number of (2+8)×30000 simulations are necessary to estimate Sobol indices via a SALTELLI method. Up to 500 cores are used to launch the simulations, for a calculation time of 20 hours overall. The observed variances of Z f in m 2 and some interest points (1, 86 and 150) are given in Figure 6. The skin rugosity coefficient α S is by far the most influencing variable, followed by θ c the Shields critical parameter. The influence of other variables can be considered negligible.

The SALTELLI method and the PCE method for the estimation of Sobol indices are compared in the interest points for the same sample size N = 30000, as shown in Figure 8. Both methods give equal estimations of Sobol indices for N = 30000. Given that the PCE method needs a total of N runs to estimate Sobol indices, which is 10 times less than the cost of SALTELLI method, it will be used from now on.

Bifurcation's results: For the bifurcation application, a sample of size N = 12800 is used . An overall 19 hours of simulation using 800 cores is necessary. The observed variances are given in Figure 9. For this case, the observed variances seem to result only from the diameters uncertainty. All the other uncertain variables seem to be none influencing.

D. Uncertainty propagation

In this section, the impact of the parameters uncertainty on the model's response will be investigated using statistics defined in section V-C1, like the mean and the variance.

1) Channel's results:

The sensitivity analysis has shown a great difference between the observed results for an experimental and a real case. In fact, for a real case, the mean diameter of sediments seems to be the only influencing parameter, while it is nearly meaningless for an experimental case. In order to investigate the influence of the sediment diameter, a new set of N = 30000 uncertainty runs is launched, considering as an only uncertain parameter the mean diameter (other parameters are constant). The observed statistics (mean in m and variance in m 2 ) are given in Figure 11. Uncertainties seem to be more important in the upstream direction, which corresponds to a higher erosion upstream and therefore bigger chances of bottom variations. This could also be confirmed through comparing the statistical mean of the final bottom elevations to a deterministic result with the statistical mean of sediment diameters on the channel's center line, as shown in Figure 12. The statistical mean of final bottom elevation is closer to the final deterministic result in the downstream direction. To quantify the uncertainty propagation, a statistical number called variation coefficient is defined in equation 25.

c v (Z f ) = V ar(Z f ) µ(Z f ) = σ(Z f ) µ(Z f ) (25)
Variation coefficients c v are evaluated in interest points with high and low variances (points 86 and 150 in Figure 11), and their convergence investigated with the Bootstrap method, as shown in Figures 13 and14 In comparison with a variation coefficient of c v = 0.07 for the diameter's sample, the uncertainty seems to be highly propagated in some interest points as point 86. This in fact confirms that even if the diameter seems to be none influencing while conducting the sensitivity analysis, its uncertainty is still propagated. In other interest points like 150, the variation coefficient of the output is lower than that of the input. In both cases (low and high c v ) the results have converged at about N = 16000. The confidence interval (INF and SUP in Figures 13 and14) for c v = 0.36 is ±0.005 and for c v = 0.06 is ±0.001. The convergence of the statistical estimations and the narrow confidence intervals indicate that our calculations can be trusted for the chosen size N = 30000.

2) Bifurcation's results: For the Bifurcation case, there is no need to investigate the diameter as an only uncertain parameter, as it is influencing the result at about 100%. The same calculations as in sensitivity analysis have been exploited to estimate the statistics in Figure 15. The highest variances seem to be concentrated in specific areas (for example, interest points 317 and 313 in Figure 15). This could be interpreted as the result of higher movements in this areas, which is demonstrated in Figure 16, where the initial and final bottom elevations have been shown for point 317 with a diameter of d = 0.149mm. This diameter corresponds to the minimum value of the studied sample. For the maximum value, no bed evolution has been observed. This difference of movement is explained by the threshold defined by the Shields parameter, which most certainly doesn't allow the maximum diameters sediments to be transported. For d = 0.149mm, a sedimentation area around point 317 is observed, which is confirmed by the bottom elevation values, going from Z f = -1.68m in the beginning of simulation to Z f = 6.27m at the end.

Nevertheless, interpreting the variances using its minimummaximum scale can give a false idea about where variances shall appear. In fact, one can interpret that there are no other variances areas than close to points 317 and 313, which is not true and can be proved with rescaling the variances as shown in Figure 17.

Furthermore, the variation coefficients c v have also been estimated for this case in low and high variance points and their convergence studied. The same conclusions can be drawn. For a variation coefficient of the diameter being c v = 0.42, uncertainty is more or less propagated in different points. For example, in point 317 (Figure 18), c v = 1.879 which corresponds to a high uncertainty propagation. For point 55 (Figure 19), c v = 0.137 which is far less than the variation coefficient of the diameter's sample and corresponds to an uncertainty that is not propagated.

E. Correlations impact

Several physical relationships can exist between the considered uncertain parameters. These correlations between the variables can impact the uncertainty study, given that the parameters are no longer sampled independently. In order to model the correlations, copulas are introduced in section V-E1. For the sensitivity study, the ANOVA method can no longer be used, because of the independent parameters hypothesis it implies. A new method called ANCOVA (ANalysis of COVAriance) is introduced in section V-E2. The uncertainty study is conducted with the correlations consideration on the channel's case. In fact, given that the only influencing parameter in the Bifurcation is the diameter, it has been concluded that the correlation study for this case would be useless.

1) Copulas for correlations modelling:

A copula is a function that defines a dependency structure between a set of variables [START_REF] Chastaing | Indices de Sobol généralisés pour variables dépendantes[END_REF]. In fact, it links the multivariate probability density function of random set of variables (X 1 , ..., X V ) to their univariate probability density functions.

A copula is a V-dimensional function C defined on [0, 1] V that verifies:

• ∀u ∈ [0, 1] V ∀i ∈ [|1 : V |], if u i = 0 then C(u) = 0 • ∀i ∈ [|1 : V |] and u i ∈ [0, 1], C(1, ..., 1, u i , 1, ..., 1) = u i • ∀u, v ∈ [0, 1] V verifying ∀i ∈ [|1 : V |] u i ≤ v i then V C ([u, v]) ≥ 0 Where V C ([u, v]) is the C-volume of the space [u 1 , v 1 ] ⊗ . . . ⊗ [u V , v V ]
defined as follows :

V C ([u, v]) = ∆ v V u V ...∆ v1 u1 C(w) ∆ vi
ui being the i th finite differential: ∆ vi ui C(w) = C(w 1 , ...w i , v i , w i+1 , ..., w V ) -C(w 1 , ...w i , u i , w i+1 , ..., w V )

The Sklar theorem [START_REF] Chastaing | Indices de Sobol généralisés pour variables dépendantes[END_REF] allows to define a relation between the multivariate PDF f X of the set X = (X 1 , ..., X V ) and the univariate probability density functions f i of X i as follows :

f X (x 1 , ..., x V ) = c(F 1 (x 1 ), ..., F V (x V )) × V i=1 f i (x i ) (26)
Where F i are the univariate cumulative distribution functions of X i associated to the probability density functions f i . On the other hand, c is the probability density function of the copula C defined as follows:

∀u ∈ [0, 1] V c(u 1 , ..., u V ) = ∂ V C ∂u 1 ...∂u V (u 1 , ..., u V ) (27)
In this study, a classical Gaussian copula is used [START_REF]Airbus-EDF-IMACS-Phimeca[END_REF]. It requires the calculation of a correlation Matrix using Spearman 23rd TELEMAC-MASCARET User Club Paris, France, 11-13 October, 2016 indices [START_REF]Airbus-EDF-IMACS-Phimeca[END_REF] from the relationships between variables. The relationships modelled here are the following:

• The empirical relation between the Meyer-Peter and Müller coefficient and the Shields number, that was introduced by Wiberg and Smith [START_REF] Wiberg | Model for calculating bed load transport of sediment[END_REF] as follows:

α M P M = 9.64 × θ 0.166 (28)
Which implies a correlation between α M P M , α S and d.

• The relationship between the Shields critical parameter and the dimensionless sediments diameter [START_REF] Tassi | SISYPHE USER'S MANUAL -Release 6[END_REF].

• The modified Komura porosity formula [START_REF] Wu | Formulas for sediment porosity and settling velocity[END_REF]:

λ = 0.13 + 0.21 (d + 0.002) 0.21 (29)
• The following relationship between the deviation parameter and the sediments diameter [START_REF] Van Ormondt | FThe effects of bed slope and wave skewness on sediment transport and morphology[END_REF]:

β 2 = 9 d H 0.3 (30)
2) Analysis of covariance: For dependent variables, it is possible to calculate the variance with the ANCOVA decomposition as follows:

V ar[Y ] = u⊆{1,...,V } [V ar[M u (X u )] + v⊆{1,...,V },v∩u=∅ Cov [M v (X v ), M u (X u )]] (31)
New sensitivity indices can be defined as:

                   S U i = V ar [M i (X i )] V ar[Y ] S C i = Cov M i (X i ), v⊆{1,...,V },v∩{i}=∅ M v (X v ) V ar[Y ] S i = S U i + S C i = Cov [M i (X i ), Y ] V ar[Y ] (32)
Where S i is the total influence of the variable X i , S U i the uncorrelated part of influence and S C i the correlated part. ANCOVA indices can be negative because of the covariance term. In order to interpret the signification of negative indices, their absolute values are compared. In fact, if |S C i | has a high value, this means that S U i is close to S i , which signifies that correlations of the variable X i have weak influence on the result. Inversely, if it has a high value, this means that correlations of X i have great impact on the simulation's result.

Lastly, as show in section V-C2, the terms M i (X i ) can be estimated via the polynomial chaos expansion. However, in order to guarantee the orthogonality of the polynomial chaos basis, it is necessary to estimate the coefficients of the expansion using uncorrelated entries X. The M i terms are estimated afterwards using the correlated values of the entries.

3) Results: The ANCOVA indices are compared to the Sobol (ANOVA) ones in order to quantify the impact of correlations on the sensitivity analysis, as shown in Figure 20. In this case, the global order of influence doesn't change. The gap between the variables in terms of influence is reduced, for example between the variables α M P M and θ C . Furthermore, some of the variables that seemed to be of weak influence at first are now considerably influencing, for example α M P M . Other variables for which no correlation was considered didn't change their sensitivity indices (eg. α C ). In order to estimate the part of the correlations influence in the total ANCOVA influence, the comparative Figure 21 is drawn. Other variables that were initially influencing also have part of their influence that is due to the correlations (θ C and α S ).

Furthermore, the variances with independent uncertain parameters and correlated uncertain parameters are compared as shown in Figure 22.

It is clearly noticed that variances considerably increase when adding correlations in this case. It can be explained by observing the new variation interval of α M P M for example.

VI. CONCLUSION AND PERSPECTIVES

In this study, an uncertainty quantification of a morphodynamic problem has been proposed.

In a sensitivity analysis step, differences were observed between a real case and an experimental case. In fact, for the real bifurcation case, the diameter was the only influencing parameter, which is not the same for the channel.

In order to analyse the influence of sediments diameter on the model's response, an uncertainty propagation study was conducted, considering as an only uncertain parameter the sediments diameter. This study has shown that the diameter has highly propagated uncertainties where there is movement. In fact, high variances were observed in maximum erosion points for the channel, and in a deposition zone for the bifurcation.

Finally, correlations were added and increased the variances significantly. An ANCOVA method was implemented in order to conduct a sensitivity analysis. The gap between the variables' influences decreased and variables that seemed first none-influencing (sediments diameter and transport coefficient of the Meyer-Peter and Müller formula) became of considerable influence when adding correlations.

This study can be generalized to other applications, such as the use of different sediment transport formulas, the study of suspended sediment transport or the influence of different physical phenomenons, for example waves (TOMAWAC module in the TELEMAC-MASCARET system).

Fig. 2 :

 2 Fig. 2: Steps for an uncertainty study

Fig. 3 :Fig. 4 :

 34 Fig. 3: Ashida et Michiue channel experiment representation

  and Müller transport coefficient Uniform[5,15] θc Shields critical parameter Uniform[0.03,0.06] λ Porosity Uniform[0.25,0.4] α S Skin friction coefficient Uniform[1.0,6.6] Φ S Angle of repose of the sediments (Slope effect on transport direction) Uniform[30,80] β2 Deviation parameter (Slope effect on transport amplitude) Uniform[0.1,5.0] αc Secondary currents coefficient Uniform[0.75,1.0] TABLE I: Uncertain parameters and their PDFs C. Sensitivity analysis 1) Monte-Carlo Sampling and statistical estimations: The Monte Carlo method requires random generation of input variables from their probability distributions. The resulted sampling of a given size N is a N × V matrix, V being the number of uncertain parameters. Each row of the matrix x i = (x 1 , ..., x V ) i represents a possible configuration for the coupled hydro-morphodynamic simulation. Corresponding realizations of the output are generated by successive deterministic simulations with each configuration of the inputs. Statistical estimators of the response Y = (Y 1 , ..., Y N ) = (M (x i )) i∈{1,...,N } can therefore be computed from the output as follows :

Fig. 5 :

 5 Fig.5: Bootstrap algorithm[START_REF] Kadiri-Ottmani | Traitement statistique d'un échantillon[END_REF]

Fig. 6 :

 6 Fig. 6: Channel case -Z f variances implied by 8 uncertain parameters First order and total Sobol indices are therefore calculated using the SALTELLI method, and their confidence intervals estimated via Bootstrap. The results are given for the point 86 in Figure 7. Similar results are observed for the other interest points.

Fig. 7 :

 7 Fig. 7: Channel case -Sobol indices for Point 86

Fig. 8 :

 8 Fig. 8: Channel case -Comparison of Sobol indices calculation using SALTELLI and PCE methods for Point 1

Fig. 9 :Fig. 10 :

 910 Fig. 9: Bifurcation -Variances of Z fSobol indices are estimated in interest points like 317 (Figure9) using the PCE method, as shown in Figure10.

Fig. 11 :

 11 Fig. 11: Channel case -Statistics of Z f for N = 30000

Fig. 12 :

 12 Fig. 12: Evolution of the bottom in the channel's center line

  .

Fig. 13 :Fig. 14 :

 1314 Fig. 13: Variation coefficient for point 86

Fig. 15 :

 15 Fig. 15: Bifurcation -Statistics for Z f -N = 12800

Fig. 16 :

 16 Fig. 16: Z f in point 317 for d=0.149mm -First (a) and last iteration (b)

Fig. 17 :

 17 Fig. 17: Bifurcation -reduced scale of variations for Z f

Fig. 18 :

 18 Fig. 18: Bifurcation -Variation coefficient point 317The variation coefficients converge at about N = 8000 for both cases, which signifies that the results of this study can be

Fig. 20 :

 20 Fig. 20: Sensitivity indices for point 86

Fig. 21 :

 21 Fig. 21: ANCOVA Indices Vs correlated ANCOVA indices

Fig. 22 :

 22 Fig. 22: Uncertainty propagation without (a) and with (b) correlations for the channel's case


Domenico Napoletani
email: napo-leta@chapman.edu

Marco Panza
email: marco.panza@univ-paris1-fr

Daniele C Struppa
email: struppa@chapman.edu

Daniele C Struppa Schmid

Forcing optimality and Brandt's principle

We argue that many optimization methods can be viewed as representatives of "forcing", a methodological approach that attempts to bridge the gap between data and mathematics on the basis of an a priori trust in the power of a mathematical technique, even when detailed, credible models of a phenomenon are lacking or do not justify the use of this technique. In particular, we show that forcing is implied in particle swarms optimization methods, and in modeling image processing problems through optimization. From these considerations, we extrapolate a principle for general data analysis methods, what we call 'Brandt's principle', namely the assumption that an algorithm that approaches a steady state in its output has found a solution to a problem, or needs to be replaced. We finally propose that biological systems, and other phenomena that respect general rules of morphogenesis, are a natural setting for the application of this principle.

Introduction

In a series of previous papers [START_REF] Napoletani | Agnostic science. Towards a philosophy of data analysis[END_REF][START_REF] Napoletani | Artificial diamonds are still diamonds[END_REF][START_REF] Napoletani | Processes rather than descriptions?[END_REF][START_REF] Napoletani | Is big data enough? A reflection on the changing role of mathematics in applications[END_REF] we described what we called the 'microarray paradigm' and we showed that there are specific methodological motifs that structure the approach of data analysis to scientific problems. By 'microarray paradigm' we referred to the belief that sufficiently large data collected from a phe-nomenon allows answering any question about the phenomenon itself. Answers are then found through a process of automatic fitting of the data to models that do not carry any structural understanding beyond the actual solution of the problem. This is a practice we suggested to label 'agnostic science'.

We argued as well in [START_REF] Napoletani | Agnostic science. Towards a philosophy of data analysis[END_REF] that, in data analysis, mathematics is "forced" onto the data. By this we mean that there are techniques expected to be useful, even when the assumptions under which these techniques should be applied do not appear to hold for the phenomenon under study. This process, which we called 'forcing', can be viewed as a direct, coarse, and willful attempt to bridge the gap between data and mathematics. This agnostic approach displays a role of mathematics in science that is essentially different from the traditional one: rather than offering a structured interface between our problems and the raw data, mathematics now provides a rich repository of techniques for forcing. The extent to which forcing is possible, and therefore the relation between specific classes of phenomena and specific mathematical techniques to force onto them, suggests a shift in the focus of our understanding: from the phenomenon itself and its inner structure, to the structure of the algorithmic processes that are privileged in data analysis. The consequence is that the link between reasonable heuristic assumptions and successful problem solving through data analysis can be broken without impact on the effectiveness of the problem solving itself.

Here, we will show the implications of this shift by dissecting one of its most pervasive aspects: the search for optimization and the belief that optimization as such is always useful and necessary. For our purposes, to solve an optimization problem means to find the minimum of a fitness function F(x) on a given domain S (usually a subset of R n for some n > 0). We explore how forcing appears in the ways the fitness function F(x) is built starting from a specific empirical problem, and in the methods used to approach the optimization problem itself.

We also propose that the way optimization techniques are used in practice often hints at a more basic methodological principle in data analysis, what we call 'Brandt's Principle', which can be expected to have wide applicability in problems from life and social sciences. The articulation of this suggestion goes through three sections. In Sect. 2 we explore how optimization can be seen as forcing. In Sect. 3, we draw the extreme consequence of regarding optimization as forcing; we suggest that, quite often, the local behavior of optimization algorithms is more important than its ultimate convergence to any optimal solution to a problem, and we introduce Brandt's principle, an operative principle of this approach. In Sect. 4, we speculate on the reasons of the effectiveness of this principle, particularly when applied to solve problems about what we suggest to call 'historical phenomena', i.e. phenomena significantly constrained by their past development. To this end, we first analyze a recently proposed principle of morphogenesis in developmental biology [START_REF] Minelli | A principle of Developmental Inertia[END_REF], and then explore its strong analogy with Brandt's principle.

Forcing optimality

When using optimization to solve an empirical problem, there are two levels at which we may have forcing.

At one level, it is possible that the optimization method is forced upon the empirical problem: though the fitness function may or may not encode significant information on the problem and the underlying phenomenon, the method used to solve the optimization problem does not actually find the global extrema of the function, but nevertheless the original empirical problem is solved by the (non-optimal) output of the algorithm. To show an example of this situation, we analyze in Sect. 2.1 a popular method of optimization in data analysis, particle swarms optimization (PSO). Specifically, we show that its effectiveness is not bound to the actual ability of this algorithm to find a globally optimal solution to a problem, rather, PSO can be seen as a collection of local explorations of the solution space.

At another level, the fitness function is forced upon the empirical problem, even when poorly related to it. The more complex the empirical problem, the more the associated fitness function will be ad hoc, and many fitness functions, not all leading to the same optimum solutions, are possible. And yet, despite the ad hoc nature of the fitness function, the empirical problem is successfully solved. To explore this second level of forcing, we show in Sect. 2.2 how important problems, such as image processing, lead to entire families of fitness functions with different choices of parameters.

In Sect. 2.2 we also show that the complexity of these problems and associated fitness functions is such that the whole variety of optimization methods used to solve the corresponding optimization problems often reduce to a local, point by point search for solutions, with no hope in general for the identification of the global optimum, similarly to what we argued in Sect. 2.1 for the restricted case of PSO methods. We can imply therefore that forcing the fitness function usually leads to forcing the optimization method itself.

We conclude that, in data analysis, what is generally relevant in the use of optimization is not a special significance of global optimum, but rather its capacity to give mathematical form to the problem itself, and its effectiveness in generating local search algorithms acting on the space of potential solutions to the problem. It is in light of this conclusion that we claim that in data analysis, the following is typically true:

Using optimization methods to solve empirical problems is a form of forcing.

In the remaining part of the present section, we will support this claim through examples. This will allow us to distinguish the two relevant forms of forcing optimality most frequently used in data analysis. We are then left to understand in Sect. 3 and Sect. 4 the reasons of the effectiveness of forcing optimality.

Forcing optimization methods

In one of the early examples of computational methods inspired by biological systems, Kennedy and Eberhart motivated an innovative optimization method in [START_REF] Kennedy | Particle swarm optimization[END_REF] by assuming that birds move in a flock to optimize their search for food. Modifying some preexisting models of birds flocking behavior [START_REF] Reynolds | Flocks, herds, and schools: A distributed behavioral model[END_REF], they built virtual particles that try to stay cohesive, while individually searching for optimality of a fitness function. The resulting searching method, called particle swarm optimization (PSO), turned out to be extremely popular in solving optimization problems, partly because of the attractiveness of a technique, such as PSO, that does not require much effort to be adapted to specific problems.

Here we want to show how the structure of this technique relies heavily on heuristic rules, to the point that PSO loses its original goal of finding a global optimum of a fitness function. PSO methods starts by setting (possibly randomly) position and velocity of multiple virtual particles, each of them associated with the evaluation of the fitness function at its position. The position and velocity of each particle is then slightly modified, taking into consideration both the location of its best evaluation up to that moment, and the location of the best evaluation of a set of neighboring particles. This updating process is justified by the hope that, as the particles move in the domain of the fitness function, and communicate their respective evaluations with each other, at least one of them will eventually settle its position to the sought for location of the global minimum. To which extent this assumption is true will be discussed shortly, after giving some details on the actual implementation of PSO methods.

Consider a swarm composed by N distinct particles, each of which is defined by a position x i = x i (t) ∈ R n and a velocity v i (t) ∈ R n , (i = 1, ..., N), depending on the value of a temporal parameter t. Let f (x), with x ∈ R n , be a fitness function that we wish to minimize on a domain S ⊂ R n , and let x i (0) ∈ R n and v i (0) ∈ R n be the initial position and the initial velocity vector of the i-th particle, respectively.

We define b i (t) as the best position that the i-th particle achieved up to time t, which means that we have f (b i (t)) ≤ f (x i (t )) for all t ≤ t and b i (t) = x i (T ) for some T ≤ t. Assuming that each particle has a well defined neighborhood N i (either a small set of particles within a given distance, or a set of particles that are a priori considered to be close to each other), we also define l i (t) as the best position within the neighborhood, which means that f (l i (t)) ≤ f (b j (t)) for all j ∈ N i and l i (t) = x j (T ) for some T ≤ t and j ∈ N i . We then update the velocity and position of each particle according to the following rules [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Dorigo | Particle swarm optimization[END_REF]:

v i (t + 1) = wv i (t) + ϕ 1 U 1 (t)(b i (t) -x i (t)) + ϕ 2 U 2 (t)(l i (t) -x i (t)) x i (t + 1) = x i (t) + v i (t + 1)
where w > 0 is a parameter called 'inertia weight' which gives a measure of the tendency of a particle to keep its direction of motion, ϕ 1 , ϕ 2 are also positive and are called acceleration coefficients, and U 1 ,U 2 are random diagonal matrices with values in the interval [0, 1) regenerated at each iteration. The choice of w, ϕ 1 , and ϕ 2 is crucial to ensure that the algorithm output does not diverge. The whole term ϕ 1 U 1 (t)(b i (t)x i (t)) is referred to as 'cognitive component', as it quantifies a sort of memory of the particle, that keeps track of the location of its own best past or current evaluation of the fitness function. The effect of this term is to contribute to the velocity vector v i (t + 1) a component that always points to b i (t), and that, the further x i (t) is from b i (t) the larger it is (not withstanding the variability of U 1 (t)). The term ϕ 2 U 2 (t)(l i (t)x i (t)) is referred to as 'social component', since it involves a local comparison within a set of neighboring particles, and it gives the particle a tendency to stay close to the location of the best past or current evaluation of the entire set of neighbors. Like the previous term, it contributes to the velocity vector a component that is always pointing in the direction of the best location ever found by the neighbors.

The purpose of local optimization is encoded in the structure of the PSO algorithm through these two components that decide the direction of the movement of the particle. The halting of the algorithm updating the positions and velocities of the particles is usually done either by setting a very large number of iterations before running the algorithm, or by terminating the iterative updating when the velocities of all or most particles get close to zero (presumably because they are close to the global minimum). In both these cases the estimate of the global minimum will be the minimum among the values of the fitness function at the location of the particles. The use of several particles allows to explore widely the fitness function on the domain S , while maintaining a certain coherence among their evolution. However, there is no evidence in general that running the PSO algorithm will indeed lead to some global optimum in S .

While there are many variants of PSO methods, they broadly respect the structure sketched above, and they turn out to be effective on a large variety of optimization problems, including many in image and signal processing, and in control theory [START_REF] Poli | Analysis of the publications on the applications of particle swarm optimisation[END_REF]. Here we broadly summarize an image processing application that approaches the issue of image denoising using PSO methods [START_REF] Bhutada | PSO-based learning of sub-band adaptive thresholding function for image denoising[END_REF]. This is a particularly interesting example because it will introduce a way to build fitness functions that we will further explore in Sect. 2.2.

Let I be a discrete image, given as a real valued function defined on a grid of N discrete points. It is often more useful to represent I as a weighted sum of simpler images, that encode some specific characteristics we may expect in the image, for example, such simpler images could be oscillating patterns with fixed frequency, or they could be sharp spikes that are nonzero only over a few points on the grid. Let us call such collection of images a dictionary G = {g 1 , ...g M }, and assume that I = ∑ M m=1 c m g m for some coefficients c m , so that, in a given dictionary G , our image can be represented by the list of coefficients c = {c 1 , ...c M }. We note that in general we need more dictionary elements than the size of the image, i.e. M > N (see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], chapter 5). In practice the image could be contaminated by noise, and in an idealized scenario we can assume that such noise is additive, that is, we have access only to the noisy image Ĩ = I +W , where W is a set of values drawn from some probability distribution. Under these conditions, the dictionary coefficients that we have access to are those of the noisy image Ĩ, i.e. c = { c1 , ... cM }.

How can we reconstruct the original image from these noisy coefficients? It turns out that, if the dictionary is well chosen, we can set to zero or shrink the smallest coefficients (threshold them) without losing the most important characteristics of the original image, for example its edges and texture. Wavelet representations are among the dictionaries that allow this thresholding operation [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Without going into the details of their structure, we can say that they entail a way to decompose the image that encodes information at large scales (say, the rough contours of the objects in the image) and information at much finer scales (information on local textures and noise). Then it is possible to build a threshold function T ( cm ) that, according to the magnitude of each cm , decides whether that coefficient needs to be replaced by zero, or, possibly, by some other value smaller than cm . For example, early in the development of this field, it was proposed the threshold function defined by [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF]. This turned out to be too crude in some applications where more information is available on the type of images of interest, and more complex threshold functions, depending non linearly on several extra parameters, were suggested. In particular one recent type, developed in [START_REF] Nasri | Image denoising in the wavelet domain using a new adaptive thresholding function[END_REF] and used in [START_REF] Bhutada | PSO-based learning of sub-band adaptive thresholding function for image denoising[END_REF], depends from three parameters p, q, r. We denote it here by T ( cm , p, q, r). This threshold function can be seen, for each choice of parameters p, q, r, as a piecewise smooth function very close to zero in a symmetric neighborhood close to the origin, and almost linear outside that neighborhood.

T ( cm ) = 0 if | cm | < C, T ( cm ) = cm if | cm | > C, for a carefully chosen value of C [
The parameters in T ( cm , p, q, r) are chosen to keep the new reconstructed image as close as possible to the original one, and one way to do that is to assume that the error (fitness) function E(p, q, r)

= I -∑ M m=1 T ( cm , p, q, r)g m 2 2
, is minimized, where we denote by * 2 2 the square of the L 2 norm of a discrete image (i.e. the sum of the squares of the values of the image at each of the N points on which it is defined). Alternatively, under some restrictions on the type of dictionary used, one could minimize the simpler error function E(p, q, r) = ∑ M m=1 |T ( cm , p, q, r)c m | 2 , as it is done in [START_REF] Bhutada | PSO-based learning of sub-band adaptive thresholding function for image denoising[END_REF]. Assuming we know one representative image with and without noise, it is possible to minimize the function E(p, q, r) above, with the hope that, once optimal parameters p, q, r are found, the threshold T ( cm , p, q, r) with that choice of parameters will effectively denoise other noisy images. Now, the complex dependance of the threshold function from its parameters causes a proliferation of local minima for the error function E. Having several local minima makes it very difficult to have a good, single initial estimate for the approximate location of the minimum, that is, a selection of values p 0 , q 0 , r 0 that are close enough to the globally optimal values. Such good estimate would be necessary if the global minimum on the domain S is searched with more traditional gradient descent methods ( [START_REF] Arfken | Mathematical Methods for Physicists[END_REF], Sect. 7.3), which essentially modify the location of an evaluation of a fitness function in the direction of the steepest negative change of the function. What is more troubling, these methods depend on an estimation of the derivative of the error fitness function, which, in turn, depends from the coefficients c m of the original image (with potentially very wide variations) making a numerical estimate of the derivative difficult.

PSO methods can overcome these difficulties thanks to the following two properties: they do not need a good initial estimate for the global minimum, since an entire population of initial evaluations is used instead; they do not need the computation of the derivatives of the error fitness function, and therefore their performance is robust with respect to variations in the image. In fact, in [START_REF] Bhutada | PSO-based learning of sub-band adaptive thresholding function for image denoising[END_REF] it is shown that PSO methods can be used successfully and quickly to minimize the error function E, either in the basic setting of uniform threshold across all the dictionary, or even the more computationally intensive sub-band setting, where the parameters are chosen independently for subgroups of the dictionary elements. The threshold parameters found by PSO are then shown to give effective image reconstructions, when confronted with other denoising techniques.

The preponderance of applications of PSO methods in image processing is not casual, as optimization problems derived in this setting often display the complex dependance from their parameters argued in the image denoising example (see Sect. 2.2). But despite the wide applicability of these methods, they are known to perform poorly when applied to combinatorial optimization problems [START_REF] Poli | Analysis of the publications on the applications of particle swarm optimisation[END_REF] where the true optimum is one of a finite number of discrete possibilities that, possibly, grows exponentially with respect to the number of the variables involved in the problem. However, these optimization problems are exactly those that are hardest to solve (many combinatorial optimization problems are NP-complete [START_REF] Wegener | Complexity Theory: Exploring the Limits of Efficient Algorithms[END_REF]), and for PSO methods to be a significant advance in optimization we would expect them to perform well when applied to these problems. As noted in the conclusion of [START_REF] Poli | Analysis of the publications on the applications of particle swarm optimisation[END_REF]: most people will not care as to whether their new tool is guaranteed to give the absolute best performance on problem. What they want is something simple and reliable. Finally, probably the PSO has, at the moment, in the mind of many people the sort of magical black box flavor that attracted so many researchers to other area of artificial/computational intelligence before. This black box belief was briefly discussed in [START_REF] Napoletani | Artificial diamonds are still diamonds[END_REF] when we suggested that the microarray paradigm has a counterpart, in the use of computational processes, in the belief that a sufficiently complex machine should be able to solve any scientific problem.

The work in [START_REF] Pedersen | Simplifying Particle Swarm Optimization[END_REF] suggests that PSO methods can be conceptually and effectively simplified by recognizing the heuristic nature of the social and cognitive components in their implementation. Under this interpretation, they are methods of space exploration whose implementation is guided by a fitness function (to be optimized). However, they do not use in an essential way the analytical properties of the fitness function (for example its derivatives) and they are not expected to eventually converge to a global minimum. It is argued instead that only the actual performance of PSO methods on benchmark problems, where there is a consensus of what constitutes a solution, can discriminate those variants that are faster and more robust.

This means that the effectiveness of PSO methods does not depend on their having been conceived as optimization techniques. In fact, what counts here is not the optimization, but the capacity of the algorithm to appropriately modify the relevant values (position and velocity of each particle). Even if a fitness function f needs to be appealed to, for the algorithm to work, once this is done, optimization is no more a central concern of the resulting recursive algorithmic procedure. A PSO algorithm is accepted as suitable simply as long as, by following it, some solution to the original empirical problem is found, in a way that does not exceed some predetermined time and resources constraints, rather than in the best possible way.

It is true that a choice of a fitness function to be used by PSO methods might not be effective for a specific problem, and this would require it to be replaced by another function, until an appropriate fitness function is chosen. Still, this does not imply that we have a trial and error procedure, or some sort of hypothetico-deductive approach. PSO methods are shown to be effective, in a very specific algorithmic form, for entire classes of very different empirical problems even though they may not truly optimize sufficiently complex fitness functions. What happens is that a fixed mathematical technique, or better, a well specified algorithm, is forced onto the problem; if trials and errors occur, it is at the level of the fitness function, and always within the framework of PSO methods, which is mostly fixed under the passage from each trial to another. This is not to say that the fitness function itself is never forced onto the problem. Rather the contrary is often true, both in the application of PSO, and in other optimization methods used in data analysis: the specific choice of the fitness function (and its corresponding extrema) does not hinge on a structured heuristic insight into the empirical problem, which prefigures another sort of forcing. This is what we shall see in the next section.

Forcing fitness functions

The way image denoising is phrased as an optimization problem in Section 2.1 is exemplary of the way in which an optimization framework is used to phrase other complex problems about natural phenomena, irrespective of the specific method used to solve the resulting optimization problem. The emphasis in Sect. 2.1 was on how to solve an optimization problem, already derived from an empirical problem. Here instead we explore, in the context of more sophisticated image processing problems, the difficulties in the preliminary process of translating the empirical problem into an optimization one by choosing a fitness function.

A first consideration in image processing is that a natural image has important features that are relevant to its recognition and classification, for example overall contours, or different uniform textures in different regions, or different objects superimposed to each other. Capturing this structure involves first of all having a compact representation for the image, for dictionaries that efficiently encode edge information ( [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], chapter 9), this requirement can be approximately satisfied by requiring the minimization of

E 1 = ∑ M m=1 |T ( cm )|
, where by T ( cm ) we indicate here symbolically the coefficients of the reconstructed image (without being concerned with the analytical shape of the function T ), and where we enforce that Ĩ = ∑ M m=1 T ( cm )g m . Minimizing the function E 1 has been shown to lead to a sparse, compact, representation for the reconstructed image where only a few coefficients are nonzero [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF]. Suppose now that we wish to have a sparse image reconstruction and a denoising of the image at the same time. A way to proceed would be to look at the minimization of the overall fitness function

E = λ E 1 + E 2 ,
where

E 2 = Ĩ -∑ M m=1 T ( cm )g m 2 2
(a modification of the fitness function seen in Sect. 2.1 with respect to the noisy image Ĩ) and λ > 0 determines the relative strength of the two fitness terms in E [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF]. What is important for our discourse, is that there is now an entire family of fitness functions, parameterized by λ . While under certain restrictive conditions on the dictionary it is possible to select heuristically a specific good value of λ [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], this is not the case in general, and we are left with the problem of choosing which, if any, of the fitness functions is truly the most suitable for the specific noisy image to represent sparsely.

More sophisticated image processing problems are associated to ever more complex fitness functions. Consider for example the problem of identifying incomplete contours of objects in an image. One of the most efficient method to solve it goes under the name of 'active snakes' [START_REF] Kass | Snakes: Active contour models[END_REF]. These algorithms are based on the idea of superimposing contours (snakes) on the image, and of subjecting these snake contours to some internal forces (minimizing the complexity of the shape of the contours) and some external forces (that maximize the fit of the contours on features of the image). Starting from some initial snake contour guess, these forces slowly change the shape of the snake on the image until it stabilizes on some significant feature within the image (usually the actual contour of an object).

The effect on the snake of internal and external forces is determined by a fitness function depending on several parameters that encode the relative strength of the forces, and to each choice of parameters (to adjust heuristically on the specific problem) corresponds a distinct fitness function. Effectively, we have an entire family of fitness functions, all suitable in principle to solve the problem of finding contours of objects in a given image. However, the dynamics of the snake on the image can be very complex, and different choices of fitness functions lead the snakes to stabilize on very different contours.

The tentative ways in which image processing problems translate into difficult optimization problems makes it clear that fitness functions are often ad hoc, and several parameters and fitness terms constraining the optimal solution need to be determined to make them plausible candidates to solve them. The significance of any specific fitness function is weakened, even though we need to have one such function to apply optimization methods. This is in line with the general fact that any application of forcing reduces the significance of our structural understanding of a phenomenon and a related problem. Which of these potentially distinct, but all equally plausible, fitness functions should we optimize? As Brandt notes in [START_REF] Brandt | Multiscale Scientific Computation: Review[END_REF]: this combination of penalty terms creates a monstrous minimization problem [...] It is extremely difficult to solve -and unnecessarily so.

The difficulty that Brandt highlights is due to the fact that for a general fitness function F there are no closed, analytical methods to solve the related optimization problem. The only way to look for a solution is often a local search, what he calls 'point by point minimization', an iterative method that, given a point x i , gives as output a point x i+1 in a neighborhood of x i such that F(x i+1 ) ≤ F(x i ). The iterative process is started by selecting an initial guess x 0 for the potential solution and it is terminated when the distance |x i+1x i | between successive outputs of the process is smaller than some preset small δ > 0, in which case we say that the point by point minimization has converged.

We have seen this local approach to optimization already in PSO methods, and indeed for all fitness functions used in the image processing problems above the solution is usually found with one of several recursive local optimization methods, such as gradient descent and Euler-Lagrange methods [START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Bresson | Fast Global Minimization of the Active Contour/Snake Model[END_REF], or interior point methods [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF][START_REF] Boyd | Convex Optimization[END_REF]. We suggest that for a sufficiently complex optimization problem and related fitness functions, arising from any empirical problem, a computational method that attempts to solve the problem is generally bound to be some type of point by point minimization. In particular, for a general fitness function, a point by point optimization process may not converge at all, may not converge quickly enough, or may converge to a point that does not minimize the fitness function on its domain; this limitation is intrinsic to these methods, rather than being incidental to some of the specific examples we have seen. Even more important, there is no guarantee that choosing slightly different fitness functions, with different parameters, will lead to similar solutions. In the next section we will see how "unnecessarily" difficult optimization problems and fitness functions can be reinterpreted and simplified by looking closely at the structure of point by point minimization processes.

Brandt's principle

In this section we suggest that both PSO and the other methods developed by forcing optimization onto empirical problems are governed by one methodological principle, which we call 'Brandt's principle' (derived and generalized from considerations of Achi Brandt in [START_REF] Brandt | Multiscale Scientific Computation: Review[END_REF]) which clarifies the role and appropriate usage of point by point minimization. We argue moreover that this principle shows a specific, operative way in which the link between structural understanding of a phenomenon and successful problem solving through data analysis can be broken, without apparent impact on the effectiveness of the problem solving itself. In particular, we show that the actual implementation strategies of Brandt's principle, while depending on partial, fragmented information on the phenomenon, are an expression of the microarray paradigm.

The starting point of our analysis is the strategy suggested in [START_REF] Brandt | Multiscale Scientific Computation: Review[END_REF][START_REF] Brandt | Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics[END_REF] to solve image processing problems. Those reviews show how direct multiscale algorithms, that satisfy a whole set of contrasting criteria at different stages, by focusing on different scales at each of these stages, can perform as well as global regularization techniques in finding solutions of some ill posed problems (i.e. problems where the solution is not uniquely defined, and/or such that the solution is not robust under variations of the input of the problem). The relevance of multiscale methods is not surprising for image analysis, given what we have said in Sect. 2.1 concerning the way information about the relevant image are encoded differently at different scale of detail (contours versus texture for example). The main focus and interest in [START_REF] Brandt | Multiscale Scientific Computation: Review[END_REF][START_REF] Brandt | Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics[END_REF] is to understand how to approach optimization problems pertaining to image processing and to partial differential equations with multiscale methods, and how to build, case by case, such methods. Here we would like to expand the scope of this insight to wider classes of problems and methods.

According to Brandt, in many cases, the search for optimization can be replaced by "a more universal and often far easier" approach, which just consists in admitting: a solution which is just the end product of a suitable numerical process not necessarily designed to satisfy, even approximately, any one governing criterion. ( [START_REF] Brandt | Multiscale Scientific Computation: Review[END_REF], page 60)

The problem with this highly agnostic approach (in our sense) is how to identify suitable numerical processes. A criterion Brandt suggests for this purpose is this: the amount of computational work [in a numerical process] should be proportional to the amount of real physical changes in the computed system. Stalling numerical processes must be wrong. [START_REF] Brandt | Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics[END_REF] It would follow that what truly distinguishes a solution of a problem, or decides that the time has come to switch to a different solving algorithm, is a nearsteady state of the algorithm output over successive iterations of the algorithm itself; stalling algorithms are wrong either because we have found a solution and we should terminate the algorithm (as it is the case for point by point optimization algorithms that have converged), or because we need to switch to a new algorithm. We state succinctly the radical idea that underlies the two quotations above in the following principle, which we name after Achi Brandt. And we claim that its domain of applicability is far more encompassing than the field of multiscale methods, and that it can be taken as the organizing principle of most data-driven computational methods:

Brandt's principle: An algorithm that approaches a steady state in its output has found a solution to a problem, or needs to be replaced. This principle proposes a fundamental shift away from the search of best solutions to problems. This shift underscores the data-driven nature of the methods in data analysis, and complements it. Since Brandt's principle does not have a built-in criterion to establish when a solution has been found (instead of having reached a stage in which one has to switch from an algorithm to another), the output of the relevant numerical process is to be checked, by external validation means, for its usefulness: a fuzzy concept in general, which nevertheless does have quite precise renderings in many specific applications. For many empirical problems, a precise task has to be achieved, and any method that can achieve such a task would be suitable. Think, for example, of the problem of having an autonomous, driverless car travel from one location to another in rough terrain; achieving this task, in any way, is the really interesting problem, because of the dangers of the car getting overtopped, or otherwise incapacitated. Finding the very best route in terms of shortest distance, or fuel efficiency is secondary. In general, the distinction we make here is between optimality, which usually designate some singular points in the space of solutions, and usefulness, which rather designate entire regions in the space of solutions. While Brandt's principle deemphasizes the search for optimal solutions to problems, it is satisfied by local, point by point optimization methods. We have already seen in Sect. 2 that wide classes of computational methods, including PSO and image processing methods, can ultimately be expressed in this way. Moreover, most data analysis methods attempt to fit a model on data by minimizing some suitable error function [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Because the fitness functions that correspond to such methods are usually both ad hoc and complex, local minimization techniques are preferred in this context. This is true in particular for artificial neural networks, currently some of the most powerful data analysis techniques ( [START_REF] Izenman | Modern Multivariate Statistical Techniques: Regression, Classification, and Manifold Learning[END_REF], chapter 10). As we pointed out in [START_REF] Napoletani | Processes rather than descriptions?[END_REF], the success of artificial neural networks can be ascribed in great part to the effectiveness of a particular local recursive process to estimate their parameters, the backpropagation algorithm. In turn, this algorithm is essentially a form of gradient descent algorithm ( [START_REF] Hastie | The Elements of Statistical Learning[END_REF], Sect. 11.4).

In light of the pervasiveness of local optimization techniques in data analysis, we can therefore claim that most of its methods satisfy Brandt's principle when implemented in an effective algorithmic form. Moreover, some of the most important types of classification methods, such as committee machines ([20] chapter 14) can be shown to satisfy Brandt's principle without an explicit local optimization implementation. This is particularly significant when we consider that boosting algorithms, singled out in [START_REF] Napoletani | Agnostic science. Towards a philosophy of data analysis[END_REF] as one of the most significant embodiments of the microarray paradigm, are a type of committee machines ( [START_REF] Izenman | Modern Multivariate Statistical Techniques: Regression, Classification, and Manifold Learning[END_REF], Sect. 14.3). In fact, committee machines are built starting from more basic, structurally similar classification algorithms (classifiers). A standard committee machine switches among these basic classifiers, each attempting to find a solution for a slightly different version of some initial classification problem, and it eventually combines them into a single, potentially more accurate classifier (we refer to [START_REF] Napoletani | Agnostic science. Towards a philosophy of data analysis[END_REF] for a general methodological analysis of classification problems and classifiers). Crucially, the decision of switching among the basic classifiers is made on the basis of the convergence to a stable tentative solution to the classification problem at each iteration of the process. A similar, sequential approach is also used in hybrid classification methods [START_REF] Delimata | Rough Sets and Intelligent Systems -Professor Zdzisław Pawlak in Memoriam[END_REF]. The difference with respect to committee machines is that the basic classifiers utilized by hybrid methods do not necessarily share a common structure.

We did not emphasize so far the aspects of an algorithm that are uniquely related to a specific problem. Brandt talks about "suitable" algorithms, but what does it mean in the context of a specific problem? To begin to answer this question, we recall that, in Sect. 2.2, the fitness functions for image processing problems are built piece by piece on the basis of partial assumptions on the images. While individually these assumptions do encode some understanding of the structure of images, the resulting fitness functions, with all their parameters to be chosen, do not offer a transparent way to interpret the solution to a given problem, and their complexity does not generally allow for a convergence to the optimum. Fitness functions are forced on the problem and they can, at best, be justified as rough quantitative a priori constraints on the form of a solution.

More generally, in image processing the a priori constraints are our assumptions on what constitutes a natural image or a contour, and these assumptions affect the choice of image dictionaries. It is possible however to make this choice, and establish our assumptions on images, automatically, for example by allowing sparsity to be the governing criterium of image representation, and using large databases of natural images to search for the best dictionary elements that, on average, optimize sparsity on a given image of the database [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF][START_REF] Field | Wavelets, vision and the statistics of natural scenes[END_REF]. This process is computationally intensive, but it can be achieved in many cases, and the resulting choices of dictionaries can then be used for further image processing problems. Of course, the assumption of sparsity is not justified from the data per se, but is required by our limited resources.

We note also that the use of specific prior assumptions on a phenomenon to force optimization is similar to the basic approach of pattern theory [START_REF] Grenander | Pattern Theory: From Representation to Inference[END_REF] and Bayesian analysis [START_REF] Berger | Statistical Decision Theory and Bayesian Analysis[END_REF] (other very powerful theories to model and classify complex phenomena), where the solution to empirical problems is contingent on the identification of appropriate priors on the probability distributions of the stochastic processes used to model the underlying phenomena.

From this analysis, we can conclude that the building blocks of methods that respect Brandt's principle are based on the unavoidable inclusion of partial pieces of prior information (what, following a customary usage, we can call 'priors', tout court) on a phenomenon. This is an important point that we need to clarify, as apparently the presence of priors proper to a particular phenomenon may be seen to be in contradiction with the microarray paradigm, and the whole agnostic approach to data analysis. One could argue indeed that, if for each problem we should identify the exact prior information that is necessary to solve it, we would have no reason to say that the microarray paradigm is at work.

Still, as we have seen in Sect. 2.1, even the use of specific priors (fitness functions in that case) does not guarantee that we can understand the reason a specific output of an optimization method solves the relevant empirical problem. If this is not the case, and it happens very often, the problem solving is still agnostic in our sense, and the microarray paradigm applies. In other terms, the fact that some priors are proper to a certain problem does not mean that their identification depends on a structural understanding of the subjacent phenomenon.

In data analysis we assume a huge quantity of data about a phenomenon, and from these data we can extract many significant features of the phenomenon, in an entirely data driven way if necessary, to the effect that these features are generally not structurally related to each other. More than that, it is often quite hard to guess in advance which subgroup of the features so detected will be useful to solve a given problem. This large number of priors can be partially structured, for example in the form of plausible guesses on the mathematical nature of the data (as in Bayesian analysis, where we may have information on the shape of probability distributions about the phenomenon), or it can be simply in the form of large collection of raw data from the phenomenon. In both cases, it is the sheer size of priors that set us squarely within the microarray paradigm: we need as many priors as possible to solve a problem in a data driven way, and to have an excess of partially structured priors is as opaque and agnostic as to have none.

We finally note that effective hybrid classification methods differ significantly in the order and the type of the basic classifiers they use. Similarly, image processing methods based on optimization allow for significantly different, equally suitable fitness functions. This suggests that, when building a method based on Brandt's principle, the specific way the priors are put is not essential: what counts is the ability to distinguish when an algorithm produces dynamically changing outputs, and when instead it gives nearly steady outputs. And the reason for the usefulness of the alternating process advocated by Brandt's principle is that it does not explore the solution space randomly, but according to priors that, no matter how they are combined, encode preferential ways to look at the problem itself. Therefore the expectations on algorithms governed by Brandt's principle are the following: they have to be easy to build, starting from a large number of weak assumptions on a phenomenon and a problem about it; their nearly steady outputs have to be exceptional and clearly identifiable; and switching among steady states outputs has to provide a fast exploration of the solution space of the problem at hand.

Data analysis of historical phenomena

We have seen in Sect. 3 how the process of forcing optimality and its effectiveness can be reinterpreted in light of Brandt's principle, and we have argued that this principle governs most agnostic, data-driven methods. In this Section we explore the appropriateness of using such methods, consistent with Brandt's principle, to solve problems about historical phenomena (introduced in [START_REF] Napoletani | Is big data enough? A reflection on the changing role of mathematics in applications[END_REF]). We perform this analysis by looking at one of the organizing principles that have been suggested in recent years to make sense of biological processes, the so called principle of developmental inertia [START_REF] Minelli | A principle of Developmental Inertia[END_REF]. We will show the generality of developmental inertia as an organizing principle of historical phenomena, subject to morphogenesis (a term used here in the general and original meaning of change of characteristics). Finally we will draw parallels between Brandt's principle and the principle of developmental inertia to find evidence of the relevance of the former (and of data analysis in general) for the solution of problems about historical phenomena.

In [START_REF] Napoletani | Is big data enough? A reflection on the changing role of mathematics in applications[END_REF] we suggested that biology and social sciences could be preferred domains of exploration of data analysis, as sciences concerned with "historical phenomena", i.e. phenomena significantly constrained by their past, historical development. We briefly argued moreover that it is possible to describe historical phenomena as "those phenomena whose development can only be constrained locally (in time and/or space) by (potentially multiple) optimization processes acting on subsets of variables, and in such a way that the functions to be optimized change over long periods of time" [START_REF] Napoletani | Is big data enough? A reflection on the changing role of mathematics in applications[END_REF]. This characterization was motivated by an analysis of fitness landscapes in evolutionary biology. There are important relations between these phenomena and Brandt's principle, in particular in the way the latter allows to reinterpret point by point minimization algorithmic processes. To understand these relations, we first explore the general structure of historical starting with ideas from developmental biology.

In [START_REF] Minelli | A principle of Developmental Inertia[END_REF] Alessandro Minelli, in order to conceptualize the way organisms structure themselves, suggests the principle of developmental inertia; generally stated, it asserts that we can see:

[biological] developmental processes...as deviations from local self-perpetuation of celllevel dynamics ( [START_REF] Minelli | A principle of Developmental Inertia[END_REF], page 119).

Similarly to the concept of inertial conditions in physics, the main purpose of the notion of developmental inertia is to identify an appropriate singular state in the relevant processes, the null (or inertial) state of these processes, this is the state that the system would stay in, if no perturbations (external to this very state) had occurred. The intrinsic nature of such null state is less relevant than the nature of the perturbation acting on it, allowing the system to evolve along different lines. It is argued in [START_REF] Minelli | A principle of Developmental Inertia[END_REF] that:

[developmental] inertial conditions [...] do not represent origin, in the ordinary meaning of the word, but only a convenient "zero" term of comparison for the study of something that happens in time-a segment of history, without prejudice of what happened before it.

The identification of suitable inertial conditions in biological systems leads to fruitful reinterpretations of asymmetry, segmentation and regeneration in complex full-grown organisms, exactly by pointing out the appropriate, null state of biological development. A particularly simple example can be seen in the context of embryo development, where the appropriate inertial state of embryos is the reproduction, symmetrical in space and indefinite in time, of the same basic tissue structure. However, in a real system, developmental inertia is constantly perturbed ( [START_REF] Minelli | A principle of Developmental Inertia[END_REF], page 123), which allows complex, inhomogeneous and asymmetrical organisms to form.

Developmental inertia may be a key principle in conceptualizing biological developmental processes, but it would need to have much wider applications to be an organizing principle of general historical phenomena. Since social phenomena are arguably, together with biological ones, the most important type of historical phenomena, the principle of developmental inertia should apply to these phenomena as well. Showing case by case the usefulness and validity of this principle in this new context would require a large survey in itself (which cannot be performed here). However, it would appear that, if we look at a social system as a collection of interacting individual agents with their own decision rules for behavior [START_REF] Epstein | Generative Social Science: Studies in Agent-Based Computational Modeling[END_REF], collective agreement among such individuals can be assumed to be an inertial state. For example, social norms have been shown to have a tendency to become entrenched and to be self-enforced [START_REF] Epstein | Learning to Be Thoughtless: Social Norms and Individual Computation[END_REF], to the effect that their entrenchment and their spreading across individuals could be taken, locally in time and space, as an inertial state in a society.

A better, more general route to justify the principle of developmental inertia for social phenomena is to rely on the broad, strong homology of social and biological systems. This homology is made explicit in accounts such as the "living systems theory" detailed in [START_REF] Miller | Living Systems[END_REF], that emphasize the notion of process in all living systems, whether biological or social [START_REF] Bailey | Living Systems Theory and Social Entropy Theory[END_REF]. Assuming its truth, it is likely that general principles from biology transfer to social phenomena, and we can conclude that the principle of developmental inertia applies to a wide range of historical phenomena, both biological and social in nature.

We also note that, at any given time, there is accumulation of structures in the states of historical phenomena; for example, accumulation of tissues differentiation in the embryo, and accumulation of legacy social norms in societies. This accumulation of structures gives clues about the inertial states that occurred along the history of a phenomenon: we can identify a structure as distinct from the whole (for example a specific tissue in an embryo) only because some inertial state (the repeated proliferation of that tissue) is replaced by some other inertial state (the differentiation into a new tissue). This suggests focusing on those characteristics of an historical phenomenon that can be ascribed to some unperturbed developmental inertial state that occurred in its past. We shall refer to them as the developmental structures of such a phenomenon.

In light of the broad applicability of the principle of developmental inertia to historical phenomena, it is striking how similar this principle is to Brandt's principle. Both principles identify appropriate null states: self-perpetuating cell dynamics for developmental inertia; near-steady state outputs of algorithms for Brandt's principle. Both principles see the breakdown and deviation from a null state as essential: to morphogenesis processes in complex organisms in the case of developmental inertia; to computational processes that solve problems in the case of Brandt's principle. Moreover, to have a near-steady state output of an algorithm implies that its (dynamical) computational processes are nearly repeating themselves, so that we can speak of self-perpetuation of dynamical processes also in the context of Brandt's principle. We believe this strong similarity is fundamental to the understanding of the effectiveness of data analysis, and we highlight it in the following proposition:

The principle of developmental inertia, as applied to historical processes, and Brandt's principle, as applied to computational processes, are homologous.

For historical phenomena, there is no compact description of the totality of developmental structures, so that it is necessary to identify and collect as many individual structures as we can to solve problems, which is in line with the microarray paradigm. Indeed, each developmental structure of an inertial state essentially differs from, and cannot be reduced to, the developmental structure due to another, contemporary or subsequent, inertial state. If it is necessary to focus on distinct developmental structures to solve a problem, these structures have to be identified, and harnessed, independently of each other. The author of [START_REF] Minelli | A principle of Developmental Inertia[END_REF] talks of distinct centers of developmental dynamics that are "everything everywhere", so that the proliferation of developmental structure does not reduce itself to a simple, unifying description. However, the collection of all such developmental structures does encode, to some extent, the state of the historical phenomenon.

Computational methods that respect Brandt' principle, or 'Brandt's methods', are uniquely suitable to take advantage of the incoherent proliferation of developmental structures, since each of these structures can be used as a prior for a corresponding algorithm, a building block for Brandt's methods. Alternating among these algorithms, when their output reaches a steady state, may be the preferred way to successfully explore the solution space of a problem. For general problems about historical phenomena, a suitable Brandt's method may be as efficient a way to search for solutions as any other method. The lack of global finality or optimality of developmental structures frustrates any attempt at faster exploration of the solution space.

As we have seen in Sect. 3, Brandt's principle offers a fruitful, theoretical scaffolding for forcing optimization and, more generally, for data analysis. The homology of Brandt's principle with the principle of developmental inertia suggests something more: data analysis, rather than being simply an heuristic, preliminary set of tools, could actually be the privileged way to approach historical phenomena and their problems with quantitative, theoretical tools.