Douglas S Gonçalves
email: douglas.goncalves@ufsc.br

Antonio Mucherino
email: mucherinoantonio.mucherino@irisa.fr

Carlile Lavor
email: clavor@ime.unicamp.br

Leo Liberti
email: liberti@lix.polytechnique.fr

B Leo Liberti

Recent advances on the interval distance geometry problem

Keywords: Distance geometry, Discretization, Molecular conformation

come

Introduction

Given a positive integer K and a simple weighted undirected graph G = (V, E, d),where d maps edges {u,v}∈E to positive interval weights [d({u,v}), d({u,v})], the Distance Geometry Problem (DGP) [START_REF] Liberti | Euclidean distance geometry and applications[END_REF] is the problem of finding a realization of the graph G in a K -dimensional Euclidean space. In other words, the DGP requires the identification of a map x : V → R K , satisfying the distance constraints:

d({u,v}) ≤ x(u) -x(v) ≤ d({u,v}), ∀{u,v}∈E, (1)
where • denotes the Euclidean norm.

A solution for ( 1) is called a realization or an embedding. In order to simplify the notation, we will use x u := x(u) and d uv := d(u,v) := d({u,v}) hereafter.

In structural biology, the problem of identifying molecular conformations from a given list of distance restraints between atom pairs is a DGP in dimension K = 3. This problem is also known in the scientific literature as the Molecular Distance Geometry Problem (MDGP).

In this particular application, the distances may be exact (i.e. d uv = duv ) or represented by a positive real-valued interval (i.e. duv > d uv > 0).

Exact distances are related to the chemical bonds whereas interval ones can be provided by experimental techniques. Such techniques include Nuclear Magnetic Resonance (NMR) [START_REF] Almeida | An overview on protein structure determination by NMR: historical and future perspectives of the use of distance geometry methods[END_REF], Förster resonance energy transfer (FRET) [START_REF] Bizien | A brief Survey of State-of-the-Art BioSAXS[END_REF] and mass spectrometry (MS) cross-linking [START_REF] Chen | Architecture of the RNA polymerase II-TFIIF complex revealed by cross-linking and mass spectrometry[END_REF].

The DGP is NP-hard [START_REF] Saxe | Embeddability of weighted graphs in k-space is strongly NP-hard[END_REF] and there exist several approaches to this problem (see [START_REF] Liberti | Euclidean distance geometry and applications[END_REF][START_REF]Distance Geometry: Theory, Methods and Applications[END_REF] and Sect. 1.1), where the DGP is reformulated as a global optimization problem on a continuous search domain, whose objective function is generally a penalty function of the distance constraints. More recently, a discrete approach to the DGP was proposed [START_REF] Liberti | Molecular Distance Geometry Methods: from Continuous to Discrete[END_REF], where the continuous domain of the optimization problem is transformed into a discrete domain.

Literature review

Distance Geometry (DG) has played a prominent part in Global Optimization (GO) insofar as it has important applications to science (e.g. protein conformation) and engineering (localization of sensor networks, structural rigidity, control of unmanned underwater vehicles and robotic arms), and it is naturally cast as a system of nonconvex constraints (Eq. ( 1))

in terms of continuous decision variables. In general, DGPs are reformulated as a minimization of constraint violations. Such reformulations have the property that the optimal objective function value is zero for feasible instances, and strictly positive for infeasible ones. Various approaches have been proposed in this journal for the general case [START_REF] Dong | A linear-time algorithm for solving the molecular distance geometry problem with exact inter-atomic distances[END_REF][START_REF] Dong | A geometric build-up algorithm for solving the molecular distance geometry problem with sparse distance data[END_REF][START_REF] An | Solving large scale molecular distance geometry problems by a smoothing technique via the Gaussian transform and d.c. programming[END_REF][START_REF] Huang | Some properties for the Euclidean distance matrix and positive semidefinite matrix completion problems[END_REF][START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF][START_REF] Liberti | Double variable neighbourhood search with smoothing for the molecular distance geometry problem[END_REF][START_REF] Moré | Distance geometry optimization for protein structures[END_REF][START_REF] Wu | An updated geometric build-up algorithm for solving the molecular distance geometry problem with sparse distance data[END_REF][START_REF] Zou | A stochastic/perturbation global optimization algorithm for distance geometry problems[END_REF], and many others on the application to protein conformation [START_REF] Costa | Discretization orders for protein side chains[END_REF][START_REF] Ferguson | Conformational searches for the global minimum of protein models[END_REF][START_REF] Grand | The application of the genetic algorithm to the minimization of potential energy functions[END_REF][START_REF] Locatelli | Minimal interatomic distance in morse clusters[END_REF][START_REF] Maranas | Global minimum potential energy conformations of small molecules[END_REF]54]. In this paper we focus on the case where the input graph is rigid, which implies that the search process has an inherently combinatorial side.

Over the years, the solution to MDGPs (DGPs arising in structural biology) have been typically attempted by employing tools such as ARIA [START_REF] Linge | ARIA: automated NOE assignment and NMR structure calculation[END_REF], CYANA [START_REF] Güntert | Automated NMR structure calculation with CYANA[END_REF] and UNIO [START_REF] Guerry | CASD-NMR 2: robust and accurate unsupervised analysis of raw NOESY spectra and protein structure determination with UNIO[END_REF],

which are all based on the Simulated Annealing (SA) meta-heuristic [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF].

While molecular conformations are generally obtained by the above methods and successively stored in databases such as the Protein Data Bank (PDB) [START_REF] Berman | The Protein Data Bank[END_REF], a second class of methods based on Nonlinear Programming (NLP) solution techniques has emerged in the last decades. A well-known example is the DGSOL algorithm [START_REF] Moré | Distance geometry optimization for protein structures[END_REF], which employs a homotopy method based on locally solving progressively finer Gaussian smoothings of the original problem. A third class of methods is based on Euclidean Distance Matrix Completion [START_REF] Alfakih | Solving Euclidean distance matrix completion problems via semidefinite programming[END_REF][START_REF] Dokmanic | Euclidean distance matrices: essential theory, algorithms, and applications[END_REF]. This is the case for the EMBED algorithm [START_REF] Crippen | Distance Geometry and Molecular Conformation[END_REF], which aims to fill in the missing distance bounds by constraint propagation of triangle and tetrangle inequalities. Thereafter, a candidate distance matrix (named dissimilarity matrix) is sampled from the completed interval distance matrix, and atom coordinates are obtained by matrix decomposition [START_REF] Dattorro | Convex Optimization and Euclidean Distance Geometry[END_REF][START_REF] Schoenberg | Remarks to Maurice Fréchet's article "Sur la définition axiomatique d'une classe d'espaces distanciés vectoriellement applicable sur l'espace de Hilbert[END_REF]. Since the dissimilarity matrix is not guaranteed to be a Euclidean Distance matrix, some of the original constraints might be violated. The last phase therefore consists in minimizing the constraint violation by local minimization, using the obtained embedding as an initial point.

A fourth class is centered around the so-called Build-Up algorithm [START_REF] Dong | A linear-time algorithm for solving the molecular distance geometry problem with exact inter-atomic distances[END_REF][START_REF] Wu | An updated geometric build-up algorithm for solving the molecular distance geometry problem with sparse distance data[END_REF]. These methods are based on the ancient idea of triangulation, used by humankind ever since navigation existed. In the context of distance geometry, where a point position is determined by the distances to it rather than the angles they subtend, this is known as "trilateration". Build-Up algorithms in dimension K = 3 attempt to place an unknown point x i (for some i ≤ n)by identifying at least four other points with known positions, and having known distances to

x i . When dealing with proteins and experimental data, the assumption of having four known exact distances to any given point may be excessively strong [START_REF] Malliavin | Distance geometry in structural biology: new perspectives[END_REF]. We point out, however, that some variants of the Build-Up algorithm overcome this limitation. For example, in order to address the uncertainty of the given distance values, the extension presented in [START_REF] Sit | Solving a generalized distance geometry problem for protein structure determination[END_REF] takes into account atomic coordinates and an unknown radius representing the uncertainty. Another variant [START_REF] Wu | Rigid versus unique determination of protein structures with geometric buildup[END_REF] partly addresses the requirement of unknown vertices having at least four adjacent vertices with known positions. This variant can find multiple valid realizations, but appears to lack the ability to finding all possible incongruent realizations.

A fifth and very important class of methods is based on solving a Semidefinite Programming (SDP) relaxation of the DGP [START_REF] Biswas | Semidefinite programming based algorithms for sensor network localization[END_REF][START_REF] Krislock | Explicit sensor network localization using semidefinite representations and facial reductions[END_REF][START_REF] So | Theory of semidefinite programming for sensor network localization[END_REF]. In particular, [START_REF] Krislock | Explicit sensor network localization using semidefinite representations and facial reductions[END_REF] exploits the cliques in the graph to reduce the size of the SDP formulation (also see [START_REF] Alipanahi | Determining protein structures from NOESY distance constraints by semidefinite programming[END_REF]). This method was shown to be able to solve NMR instances containing real data and to reconstruct conformation models that are very close to the ones available on the PDB.

The authors of this paper are among the researchers who proposed and worked on a sixth class of methods based on a combinatorial algorithm called Branch & Prune (BP) [START_REF] Liberti | A branch-and-prune algorithm for the molecular distance geometry problem[END_REF].

Protein graphs share some common properties: for example, they can be decomposed into a backbone subgraph and many side chains subgraphs [START_REF] Schlick | Molecular Modelling and Simulation: An Interdisciplinary Guide[END_REF] (these can be realized separately and then put together [START_REF] Santana | Side chain placement using estimation of distribution algorithms[END_REF]). The backbone subgraph is larger than the subgraphs related to side chains, and hence most difficult to realize. However, it also defines an order on the atoms with certain topological properties, which we formally discuss below (informally, we can say that every atom in this order has at least three predecessors which are also adjacent in the graph structure). Under this assumption, the search domain of the underlying optimization problem can be reduced to a discrete set with a tree structure [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF][START_REF] Mucherino | The discretizable distance geometry problem[END_REF], which can be searched by the BP algorithm [START_REF] Liberti | A branch-and-prune algorithm for the molecular distance geometry problem[END_REF]. If the distances are exact, BP can find all realizations of a given protein backbone graph. Although an exhaustive search in the conformation tree is worst case exponential [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF], numerical experiments have shown that BP behaves polynomially in protein-like instances [START_REF] Liberti | The discretizable molecular distance geometry problem seems easier on proteins[END_REF]. In fact, it can be proved that, in such cases, the problem is Fixed Parameter Tractable (FPT) [START_REF] Liberti | On the number of realizations of certain Henneberg graphs arising in protein conformation[END_REF]. In computational experiments, the parameter value could always be fixed at a single constant, which explains the polytime behaviour. For protein backbone instances with exact distances, BP is one of the fastest available methods, one of the most reliable, and the only one which can certifiably find all incongruent realizations.

An adaptation of the BP to the interval distance setting was proposed in [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF], where intervals are replaced with a finite set of discrete points. We refer to this BP adaptation as the interval BP (iBP). This algorithm was tested on real protein instances in [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF]. Although Several other approaches for solving DGPs can be found in the scientific literature. The interested reader can refer to [START_REF] Lavor | Recent advances on the discretizable molecular distance geometry problem[END_REF][START_REF] Liberti | Euclidean distance geometry and applications[END_REF][START_REF]Distance Geometry: Theory, Methods and Applications[END_REF].

Aim of this paper

O u rm a i nm o t i v a t i o ni nt h i sw o r ki st oi m p r o v et h eiBP algorithm proposed in [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF]f o r solving MDGPs with interval data. For this purpose, we identify the main limitations of this discrete approach in presence of interval distances and propose a new variant of iBP to find approximate solutions for interval MDGPs.

The identification of the barriers against the successful application of iBP in real settings represents an important step towards a combinatorial methodology with the following properties:

-it is specifically suitable for solving the protein conformation problem from distance restraint data;

-it can work with uncertain data, specified as interval distances provided by experimental techniques;

-it can potentially find all incongruent realizations of a given instance.

The paper is organized as follows. In Sect. 2, we define the subclass of DGP instances describing protein backbone graphs: we discuss assumptions, discretization orders, the iBP algorithm variant, pruning devices, and the parameterization of the coordinates. Section 3 presents a method, based on some interval distance constraints, which is able to reduce the set of candidate positions for certain vertices before the iBP branching phase. Section 4 addresses the main limitations in handling larger molecules in presence of interval data and presents a heuristic for finding approximate realizations. The computational results in Sect. 5

illustrate the improvements due to the proposed approaches.

A combinatorial approach

Let G = (V, E, d) be a simple weighted undirected graph representing an instance of the MDGP. In the following, vertices of V will represent atoms of the given molecule and {u,v}∈ E if the distance between the atoms u and v is available. The map d relates each edge {u,v}∈

E to a positive interval weight [d(u,v), d(u,v)].
The MDGP asks to find a realization

x : V → R 3 (see Introduction), i.e. a molecular conformation in three-dimensional space such that:

d(u,v) ≤ x u -x v ≤ d(u,v), ∀{u,v}∈E. (2)
Recall that d(u,v) and d(u,v) denote, respectively, the lower and upper bounds for the

distance d(u,v)(with d(u,v) = d(u,v)if d(u,v)is exact).
We also suppose that the given set of distances is realizable in R 3 .

In order to discretize the search domain, MDGP instances need to satisfy some particular assumptions. The main requirement is that the atoms need to be sorted in a way that there are at least three reference atoms for each of them (aside, obviously, from the first three). We if only one of the three distances is allowed to take values into a certain interval, then the intersection gives two arcs of a circle, generally disjoint, where sample points can be chosen [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF]. In both last situations, the discretization can be performed. More details are given in the next section.

iBP algorithm and discretization orders

Let G = (V, E, d) be an instance of the MDGP and let us suppose that there exist an order for the atoms v ∈ V , so that we can assign a numerical label i ∈{ 1, 2,...,|V |} to each of them. At each recursive call of the iBP algorithm, candidate positions for the current atom i are computed using the positions of the previously placed reference atoms and their distances to the atom i.

When the distances between i and its references are exact, the intersection of three spheres needs to be computed. If the reference atoms {a, b, c} are not collinear, then such an intersection results in at most two points. When this situation is verified for all atoms i > 3, then the search domain has the structure of a binary tree [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF].

However, if one of the three reference distances, say d ci , is an interval, then the two spheres centered at x a and x b need to be intersected with a spherical shell centered at x c . As a result, the intersection gives two candidate arcs (see Fig. 1). These arcs are over the dashed circle C defined by the intersection of the two spheres. When the intersection consists of two arcs, a finite number D of sample positions should be selected from each of them [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF]. This way, we still have a discrete set of possible positions for the atom i.

Therefore, the discretization strongly depends on an order for the vertices (atoms) of G satisfying specific properties. Definition 1 formalizes the assumptions mentioned above.

Definition 1

The interval Discretizable DGP in dimension 3 (iDDGP 3 )

Given a simple weighted undirected graph G = (V, E, d),whereE ′ ⊂ E is the subset of edges for which their weights are exact distances, we say that G represents an instance of the iDDGP 3 if there exists a total order on the vertices of V verifying the following conditions: Even if the tree grows exponentially (in the worst case scenario), the pruning devices allow iBP to focus the search on the feasible parts of the tree. The easiest and most efficient pruning device is the Direct Distance Feasibility (DDF) criterion [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF], which consists in verifying the ǫ-feasibility of constraints involving distances between the current atom i and previously placed atoms:

(a) G C = (V C , E C ) ≡ G[{1, 2,
1. a < i, b < i, c < i;
d(h, i) -ǫ ≤ x h -x i ≤ d(h, i) + ǫ, ∀{h, i}∈E, with h < i and h / ∈{a, b, c}. (3)
The distances involved in the above constraints are called pruning distances.

Protein backbone model: discretization orders and pruning devices

A necessary preprocessing step for solving MDGPs by this discrete approach consists in finding suitable atomic orders allowing each atom v to have at least three reference atoms.

We name such orders discretization orders [START_REF] Cassioli | Discretization vertex orders in distance geometry[END_REF]. In previous works, discretization orders have been either handcrafted [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF], or obtained by looking for paths on pseudo de Bruijn graphs consisting of cliques of G [START_REF] Mucherino | A pseudo De Bruijn graph representation for discretization orders for distance geometry[END_REF], or even automatically detected by a greedy algorithm [START_REF] Lavor | Discretization orders for distance geometry problems[END_REF][START_REF] Mucherino | On the identification of discretization orders for distance geometry with intervals[END_REF].

In fact, if we consider distances defined by bond lengths and bond angles as exact, along with the peptide plane geometry, it is possible to find orders for the protein backbone (and also for side chains [START_REF] Costa | Discretization orders for protein side chains[END_REF]) satisfying the assumptions required for the discretization. This preprocessing step can be performed efficiently, in polynomial time [START_REF] Mucherino | On the identification of discretization orders for distance geometry with intervals[END_REF], so that the necessary assumptions can be fulfilled by graphs related to proteins. However, we point out that when some additional assumptions are imposed to the searched orders, such as the consecutivity of the reference vertices, this problem becomes NP-hard [START_REF] Cassioli | Discretization vertex orders in distance geometry[END_REF]. In this work, we consider a model for the protein backbone as depicted in Fig. 2.

With the backbone atoms N, C α , and C, we also considered the attached H and H α ,andwe have included only the C β and H β atoms to represent the side chain of each amino acid. There are two exceptions for this amino acid model: the glycine, where C β and H β are replaced by one H, and the proline, which has a missing hydrogen. The backbone model in Fig. 2 only considers 3 amino-acids, but it can be repeated for all amino-acids in a longer protein sequence, because of the regular pattern defining the protein backbone.

The order depicted in Fig. 2 is the one used in our numerical experiments. On the basis of the model in Fig. 2 for the protein backbone geometry, it is possible to conceive other pruning devices [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF][START_REF] Mucherino | Influence of pruning devices on the solution of molecular distance geometry problems[END_REF] to be integrated with DDF (see, Eq. 3),basedonthe following considerations:

-Helices in proteins can be either right or left-handed. The former situation is statistically more common, because of side chains steric constraints. In this work, we do not consider side chains explicitly, but we suppose that it is possible to understand, from an analysis of the protein sequence, whether right-handed or left-handed helices are expected to be present. We call this pruning device as the chirality-based device: in some situations, it can allow for placing uniquely some atoms during the execution of the search. -The oxygen atoms in Fig. 2 are included in the model because they participate in hydrogen bonds. Each oxygen can be placed uniquely by using the exact distances with the other atoms of the peptide plane.

Computing coordinates for candidate positions

The method employed to compute the candidate positions at each recursive call of Algorithm 1 has a fundamental importance. While looking for candidate atomic positions for the atom i,it is supposed that the reference atoms {a, b, c} are already positioned. These reference atoms define a local coordinate system centered at a [START_REF] Gonçalves | Discretization orders and efficient computation of cartesian coordinates for distance geometry[END_REF]62]. This coordinate system is illustrated in Fig. 3.

Let v 1 be the vector from b to a and v 2 be the vector from b to c.Thex-axis for the system in a can be defined by v 1 , and the unit vector in this direction is x = v 1 / v 1 . Moreover, the vectorial product v 1 × v 2 gives another vector that defines the z-axis, whose corresponding unit vector is ẑ. Finally, the vectorial product x ×ẑ provides the vector that defines the y-axis (let the unit vector be ŷ). These three unit vectors are the columns of a matrix U a = x ŷ ẑ , whose role is to convert directly vertex positions from the coordinate system defined in a to the canonical system.

Once the matrix U a has been computed, the canonical Cartesian coordinates for a candidate position for the vertex i can be obtained by:

x i (ω i ) = x a + U a ⎡ ⎣ -d ai cos θi d ai sin θi cos ω i d ai sin θi sin ω i ⎤ ⎦ , (4)
where θi and ω i are the angles related to the spherical coordinates of vertex i.

We will use the symbol θi in order to refer to the angle formed by the two segments (i, a)

and (a, b), and we will use the symbol ω i to refer to the angle formed by the two planes defined by the triplets (a, b, c) and (b, a, i) (see Fig. 1). The cosine of the angles θi and Recall from Sect. 2.1 that if the three reference distances are all exact, then the three spherical shells are in fact three spheres, whose intersection gives 2 points, with probability 1 [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF]. The two points x + i and x - i correspond to two possible opposite values, ω + i and ω - i , for the angle ω i . When one of the three distances is instead represented by an interval (see Definition 1), the third sphere becomes a spherical shell, and the intersection provides two curves (see Fig. 1). These two curves correspond to two intervals,

[ω + i , ω+ i ] and [ω - i , ω- i ],
for the angle ω i . In order to discretize these intervals, a certain number of points, say D, can be chosen from the two curves.

As shown in [START_REF] Gonçalves | Discretization orders and efficient computation of cartesian coordinates for distance geometry[END_REF], the generalized procedure for the computation of atomic coordinates in Algorithm 1, based on equation (4), is very stable when working on MDGP instances related to real proteins. Moreover, equation ( 4) is also at the basis of an important technique that can be used to reduce the feasible arcs obtained by sphere intersection. This technique for arc reduction was firstly proposed in [START_REF] Gonçalves | An adaptive branching scheme for the branch & prune algorithm applied to distance geometry[END_REF]. Another approach for arc reduction, based on Clifford algebra, is presented in [START_REF] Lavor | Clifford algebra and the discretizable molecular distance geometry problem[END_REF].

Pruning distances and arc reduction

When candidate atomic positions, at each recursive call of the iBP algorithm (see Algo- increase the width of the tree; on the other hand, too small values can generate trees where no solutions can be found (all branches are pruned, because all positions, at a certain layer, are not compatible with pruning distances).

In [START_REF] Gonçalves | An adaptive branching scheme for the branch & prune algorithm applied to distance geometry[END_REF], an adaptive scheme was proposed for tailoring the branching phase of the iBP algorithm so that all computed candidate positions are feasible at the current layer. The idea is to identify, before the branching phase of the algorithm, the subset of positions on the two candidate arcs that is feasible with respect to all pruning distances to be verified on the current layer.

Let us suppose that, at the current layer i, the distance d ci is represented by the interval

[d ci , dci ]
. By using Equation ( 4), two intervals for the angle ω i can be identified:

[ω + i , ω+ i ]⊂ [0, π] and [ω - i , ω- i ]⊂[π, 2 π],
such that the distance constraints

x a -x i (ω i ) =d ai , x b -x i (ω i ) =d bi , d ci ≤ x c -x i (ω i ) ≤ dci , (5)
are satisfied.

However, there may be pruning distances, at layer i, that could be exploited for tightening these two arcs. Let us suppose there is an h ∈{j < i | j / ∈{a, b, c}}, such that the distance d hi is known and lies in the interval [d hi , dhi ]. The solution set of the inequalities

d hi ≤ x h -x i (ω i ) ≤ dhi (6)
consists of intervals for ω i that are compatible with the distance d hi .

A discussion about how to solve the inequalities (6), by using Eq. ( 4), is presented in details in [START_REF] Gonçalves | An adaptive branching scheme for the branch & prune algorithm applied to distance geometry[END_REF].

The feasible positions for the atom i can be therefore obtained by intersecting the two previously computed arcs (in bold in Fig. 1), and several spherical shells, each of them defined by considering one pruning distance between i and h < i. For each available pruning distance, other inequalities (6) can be defined and new arcs on the circle C may be identified. The final subset of C which is compatible with all available distances can be found by intersecting the arcs obtained for each pruning distance with the two initial disjoint arcs, given by Eq. ( 5).

After considering all pruning distances, i.e., after performing all intersections, the final result provides a list of arcs on C that are feasible with all the distances that can be verified at the current layer. All positions that can be taken from these arcs are feasible at the current layer: all of them generate a new branch and may serve as a reference for computing new candidate positions on deeper layers of the tree. In order to integrate the iBP algorithm with this adaptive scheme, there are two main changes to be performed on Algorithm 1. On line 8, the adaptive scheme needs to be invoked for taking into consideration the information about the pruning distances. Moreover, the use of the DDF pruning device has become unnecessary, and it should not be considered at line 20 of Algorithm 1.

It is important to remark that this adaptive scheme is not supposed to speed up the execution of the search, but rather to help in defining search trees that can actually contain solutions.

Without the use of this adaptive scheme, all sample positions selected from the two arcs obtained with the discretization may be discovered to be infeasible as soon as the DDF pruning device is invoked. The other extreme situation is instead the one where the adaptive scheme can allow us to select a subset of sample positions that all bring to the definition of a solution. Naturally, the second situation is desirable, even if, in terms of complexity, it tends to increase the total computational cost.

Limitations of the current approach: finding approximate realizations

For DGP instances where all available distances are exact, the presented discrete approach is extremely efficient, allowing for example to realize graphs having thousands of vertices in few seconds with a standard computer [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF].

However, for iDDGP 3 instances, there are some difficulties encountered by the iBP algorithm [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF], even for finding one solution. Such limitations, related to the presence of interval data in both discretization and pruning distances, are discussed in this section and a heuristic to overcome such barriers is proposed.

Sampled distances and embeddability

Recent computational experiments have shown that taking equidistant sample points on the feasible arcs (or equidistant samples from the interval distance, see Algorithm 1 in Sect. 2.1), even after the intersection with the available pruning distances (see Sect. 3), is not enough to allow the iBP algorithm to solve some MDGPs within a predefined precision. The sampled distances are taken independently in each layer of the tree and, in particular for small D values, it is not likely that they are compatible with each other and with other pruning distances available at deeper layers.

The underlying issue is related to the embeddability of a given set of distances. Suppose that we are positioning the atom i and that the interval distance [d ci , dci ] is used in the discretization. Even if we assume that there exists a distance value d * ci which is compatible with the other distances in E leading to a solution, we cannot ensure that, with a finite number D of samples taken from [d ci , dci ], the compatible distance d * ci is actually sampled.

In order to illustrate this fact, consider the following example where five points in R 3 are embedded (Fig. 4).

Suppose that the straight lines represent exact distances, and let the black lines be the According to the Cayley-Menger conditions [START_REF] Liberti | Euclidean distance geometry and applications[END_REF][START_REF] Sippl | Cayley-Menger coordinates[END_REF], for this set of distances to be real- ) is represented by the blue curve in Fig. 5.

It is easy to see that, unless the grid is sufficient refined (number of samples D is sufficient large), a valid pair of distances (d 2 14 , d 2 25 ) can be sampled with probability zero.

Long-range distance restraints

Long-range distance restraints are related to atoms that are at least four amino-acids apart in the protein sequence. Even if far in the protein sequence, some atom pairs may be in condition to be detected by an experimental technique. For example, if we consider NMR, it is typical to detect distances between atoms that are very far in the sequence, but quite close in space (≤5Å).

In case of all available distances are exact, the pruning distances efficiently guide the search in the binary tree corresponding to the discretized search space [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF][START_REF] Liberti | The discretizable molecular distance geometry problem seems easier on proteins[END_REF]. However, when interval distances are present, the search tree is no longer binary, because D samples are taken from each feasible arc. Moreover, the DDF pruning criterion To illustrate this fact, we depict in Fig. 6 the available entries of the upper triangular part of the distance matrices for two instances belonging to our set: 1FJK and 2E2F.

Notice that the almost-band structure close to the main diagonal is a consequence of the assumptions concerning the discretization. In fact, the distances between pairs (i -1, i)

and (i -2, i) are generally derived from the bond lengths and angles, while the distances (i -3, i) can be generally obtained from the analysis of the torsion angle among the quadruplet of atoms (i -3, i -2, i -1, i). Moreover, the distance (i -3, i) may be also estimated by applying an experimental technique such as NMR. Other distances that are far from the main diagonal of the matrix can be obtained applying an experimental technique.

Although 1FJK has more atoms than 2E2F, the former instance can be easily solved by iBP in a few seconds, whereas the latter cannot be solved in less than one minute (using D ≤ 20). The difficulty in solving 2E2F is related to the presence of long-range pruning distances: there are several entries in its distance matrix that are far from the diagonal.

Approximate realizations

The presence of interval distances implies uncertainty on the atomic positions obtained by sampling points in the intersection between spheres and spherical shells: even a small error introduced at tree layer i can have a relevant propagation until the layer j >> i and, when the pruning distance is finally tested, it is likely that the propagated error leads to infeasible positions for atom j.

Thus, an error introduced during the intersection discretization in a certain tree layer, might make every sampled candidate position infeasible with pruning distances in a further layer. This phenomenon is more evident when considering long-range distance restraints. One possibility to avoid pruning out all branches of the search tree, in order to obtain approximate solutions to the problem, is to relax the distance constraints related to long-range distances.

We define the set where M is a positive integer used to identify long-range distance restraints. Our relaxation consists in avoiding the application of the DDF feasibility test (Eq. 3), as well as the intersection scheme (Sect. 3), to pruning distances in L.

L ={{i, j}∈E ||i -j|≥M}, ( 7
Naturally, when such pruning distances are neglected, some information is lost and this can have an impact on the found solutions. In fact, long-range distance restraints are the main responsible for the global fold. Thus, in order to mitigate this effect, we introduce another pruning criterion based on the partial Mean Distance Error (MDE) at the current layer k:

PMDE k (X ) = 1 |J k | {i, j}∈J k ⎡ ⎣ max d i, j -x i -x j , 0 d i, j + max x i -x j -di, j , 0 di, j ⎤ ⎦ , (8)
where

J k ={{i, j}∈E | i ≤ k ∧ j ≤ k}.
Let n =|V | and note that J n = E. It is common to measure the quality of a realization by the Mean Distance Error measure:

MDE(X ) = PMDE n (X ).
Thus, by monitoring the PMDE k (X ) for k < n, we can control the quality of partial realizations. This suggests the PMDE pruning device:i fa tl a y e rk, PMDE k (X )>ε, then the candidate partial realization should be pruned. We set ε > ε,whereε is the tolerance used in DDF (Eq. 3).

When this new pruning device is introduced, a solution found by Algorithm 1 is actually an approximate solution in the sense that it satisfies all distances in E \ L (with tolerance ε), while some distances in L can be violated. However, the total MDE value for such a solution remains relatively small, because of the new pruning device based on [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF]. By applying this scheme, together with the chirality and peptide plane constraints (see Sect. 2.2), we expect that the fold of the obtained conformation mimics the fold of the true protein. This is the case for the set of instances used in the computational experiments.

Computational experiments

In this section we present some computational results on a set of artificially generated MDGP instances. Our aim is to assess the improvements on iBP (Algorithm 1) due to the integration of a set of recently proposed techniques: the pruning devices based on chirality and peptide plane geometry, described in Sect. 2.2; the arc reduction technique presented in Sect. 3;and the partial MDE pruning device introduced in Sect. 4.3.

The instances that we consider in our experiments were generated as it follows. The protein conformations were extracted from the PDB: by using the coordinates of a known conformation, all pairwise distances between atom pairs of the backbone were computed.

Then, only a small subset of all distance pairs is kept for defining an instance. The distances related to bond lengths and those that can be obtained from bond angles are considered as exact, as well as distances between atoms belonging to the same peptide plane (see

Assumptions concerning distances between hydrogens

During the generation of our instances, we rely on the premise that experimental techniques, such as NMR spectroscopy, are able to give information about all distances between hydrogen atoms that are close in space [START_REF] Lavor | On the computation of protein backbones by using artificial backbones of hydrogens[END_REF]. Moreover, these distances can be supposed to be more precise than other ones. Statistics on such distances [START_REF] Billeter | Sequential resonance assignments in protein 1 H nuclear magnetic resonance spectra. Computation of sterically allowed proton-proton distances and statistical analysis of proton-proton distances in single crystal protein conformations[END_REF][START_REF] Wüthrich | Pseudo-structures for the 20 common amino acids for use in studies of protein conformations by measurements of intramolecular proton-proton distance constraints with Nuclear Magnetic Ressonance[END_REF], with the geometry of consecutive peptide planes, validate this assumption.

We will consider therefore that all distances between hydrogens belonging to the same or to two consecutive amino-acids are available, and we suppose that they lie in an interval having width 0.1Å. Besides these distances, responsible for the local geometry, we also consider distances between hydrogens that belong to amino acids that are far in the protein sequence, but close in space. These distances are responsible for the global fold. We include these distances in our generated instances whenever they are smaller than 5Å and consider that imprecisions lead to an interval of width 1Å.

Hydrogen bonds H -O, responsible for stabilizing α-helices and β-strands, are also considered. If the distance between H and O of distinct amino-acids is greater than 1.3Å and less than 3.5Å, such a distance is included in our instances as an interval of width 1Å. All intervals have a predefined width and their extremes are randomly generated in a way that the interval contains the true distance.

Numerical results

Let us refer to the algorithm presented in [START_REF] Lavor | The interval branch-and-prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF]asiBP, while we will name "New iBP" the algorithm integrated with the new method for the computation of candidate positions (see Sect. 2.3), with the technique for arc reduction (see Sect. 3), with the chirality and peptide plane pruning devices (see Sect. 2.2), and with the pruning device introduced in Sect. 4.3.

In both iBP variants, the tolerance used in the experiments for the DDF criterion (Eq. 3)

is ε = 0.001. In new iBP, for the PMDE-based pruning device, we used ε = 0.01 and set M = 40 in definition of L (see Eq. 7). We gradually increased the number of samples D taken from the feasible arcs until the first solution is found in less than 60 s (timeout).

The numerical experiments were run in a Intel MacBook Pro, 2Ghz, 2GB RAM, and the Algorithm 1 was implemented in C programming language, compiled using GNU GCC with flag -O3.

Table 1 shows a comparison between iBP and New iBP. The number of amino acids (aa), atoms (|V |) and available distances (|E|) are given for each instance. For the two versions of iBP, the performance is evaluated by the minimum number of samples D (taken from interval arcs) necessary to find one solution, the number of recursive calls and the CPU time in seconds. The quality of the realization is assessed through the MDE. The character "*" means that the instance could not be solved in less than one minute for any value of D ≤ 20.

We can notice that the arc reduction technique presented in Sect. 3 is very effective in reducing the number of sample positions that we need to extract from the arcs in order to obtain at least one solution. This is an important improvement because it is not known a priori how many samples are sufficient to allow iBP to find a conformation. We can observe iBP uses the relaxed pruning criterion PMDE, it cannot ensure a better MDE for all the instances. For some of them we observe a better MDE which is a consequence of the other considered pruning devices. Although we relaxed some distance constraints by using PMDE, the chirality constraints helped in improving the local geometry, resulting in a better MDE.

For those instances with a worse MDE, like 2KXA, 2E2F or 2KIQ, the PMDE relaxation was, in some sense, the way to "pass-through" the long-range distance constraints and find a realization in a affordable time. Additionally, we remark that an MDE value around 10 -3 is able to guarantee a sufficient detection of the global fold of the protein. In fact, the realizations found by the New iBP are not so far from the true ones. The quality of such realizations is discussed in the next subsection.

Quality of a realization and practical usage

While looking at Table 1, a natural question emerges: how good are the realizations X with MDE(X ) ≈ 10 -3 when compared to the "true" protein ?

Since we have relaxed the distance constraints related to long-range distances, in principle, we cannot ensure that the underlying molecule is recovered. However, we will illustrate that the realization found by the New iBP gives a very good approximation of the true conformation.

First, let us take a look at the instance 2KXA. Figure 7 shows the realizations found (first found solutions) by the original iBP and the New iBP. Although the MDE of the first is smaller than the second, 10 -6 against 10 -3 , we can see that the conformations are roughly the same, except by partial reflections. The New iBP produced a right-handed helix because it contains, in its list of pruning devices, the chirality-based device. Now, let us consider the instance 2E2F. According to Table 1, the MDE for the realization found by the New iBP is approximately 10 -2 . By superimposing the realizations found with the true one from the PDB (first model), see Fig. 8, we obtain a RMSD value equal to 0.7 Å (according to TM-align [START_REF] Zhang | TM-align: a protein structure alignment algorithm based on TM-score[END_REF]).

Therefore, although the realization found by New iBP does not fit perfectly with the true conformation, it is close enough to identify its global fold, and it also can be used as a smart

Conclusion and future work

We collected in this paper the most recent and promising advances in solving the MDGP with our combinatorial approach. The main contributions presented in this paper can be summarized as follows:

1. identification of the main barriers against the successful application of the discrete approach to real MDGPs with interval data;

2. another pruning devices based on chirality and peptide planes, whose easy implementation is allowed by the model and discretization order presented in Sect. 2.2;

  u n c o r r e c t e d p r o o f J Glob Optim

123

  Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small u n c o r r e c t e d p r o o f J Glob Optimthis BP variant shows promise, its practical applicability is currently limited by the choice of discretization points.

Fig. 1

 1 Fig. 1 The two feasible arcs obtained by intersecting two spheres and one spherical shell

Algorithm 1

 1 3}] is a clique and E C ⊂ E ′ ; (b) ∀i ∈{4,...,|V |}, there exists {a, b, c} such that 123 Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small u n c o r r e c t e d p r o o f J Glob Optim The iBP algorithm.

2 .Fig. 2 A

 22 Fig. 2 A model for the protein backbone and a possible discretization order

  The first three atoms, N-C α -H α , of the first amino-acid can be used as initial clique (see assumption (a) in Definition 1) for the discretization order because the involved distances are defined by bond lengths and angles, that can be considered as exact[START_REF] Crippen | Distance Geometry and Molecular Conformation[END_REF]. Analogously, taking into account the peptide plane distances and the distances between hydrogens provided by 123 Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small u n c o r r e c t e d p r o o f J Glob Optim NMR, it is not hard to verify that assumptions (b.1)-(b.2) of Definition 1 aresatisfiedbythe order given in Fig. 2.

  For the carbon C, in fact, we can get only one (instead of two) possible positions by using N-C α -H α as reference atoms. An analogous reasoning can be applied to C β . The chirality defines the orientation of the tetrahedron formed by C, C β ,N ,C α ,H α ,w h e r eC α is the chiral center, and it can be used to avoid unnecessary branching; -The tetrahedron around C α forms a clique as well as the peptide planes[START_REF] Alipanahi | Determining protein structures from NOESY distance constraints by semidefinite programming[END_REF]. Such local structures define rigid regions of the protein backbone. Using the peptide plane clique, it is possible to find a unique position for C α . It is also possible to place N uniquely, because its relative orientation with respect to H, C and C α of the same peptide plane can be computed by taking into account the van der Waals minimum distance;

Fig. 3

 3 Fig. 3 The reference vertices a, b and c induce a system of coordinates

2 ab + d 2 ai -d 2 bi 2 d

 2222 ω i can be computed by exploiting the positions of the reference vertices a, b and c,aswellas the available distances d ai , d bi and d ci . Thus, cos ω i = cos θc,b,i -cos θa,b,i cos θa,b,c sin θa,b,i sin θa,b,c , where we consider the positive values for the sines, and cos θi = cos θb,a,i = d ab d ai .

rithm 1 )

 1 , are computed by intersecting two spheres with one spherical shell, a continuous set of positions is obtained, which generally corresponds to two disjoint arcs, related to two intervals for the corresponding torsion angle values. During a typical run of Algorithm 1, every time the reference distance d ci is represented by an interval, D equidistant samples are taken from each arc [34]. As a consequence, 2 D atomic positions are generated in total, and 2 D new branches are added to the tree, at the current layer, for every branch at the upper level. After their computation, the feasibility of each atomic position is verified. On the one hand, too large D values can drastically 123 Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small u n c o r r e c t e d p r o o f J Glob Optim

123Fig. 4

 4 Fig. 4 Realization of five points in R 3

  exact distances used in the discretization. The dashed blue lines are the interval distances (used to compute the possible positions of atoms 4 and 5) and the red straight line represents one pruning distance (that can be used to validate the possible positions for the atom 5). The associated distances are the following: d 12 = d 23 = d 24 = d 34 = d 35 = d 45 = 1, d 13 = √ 2, d 14 = √ x ∈[0.5, 2], d 15 = √ 3, d 25 = √ y ∈[0.5, 2].

Fig. 5 14 and y = d 2 25

 51425 Fig. 5 Solution set for the five point Cayley-Menger determinant with d 2 14 and d 2 25 as missing distances

  [START_REF] Almeida | An overview on protein structure determination by NMR: historical and future perspectives of the use of distance geometry methods[END_REF] becomes much less effective when the bounds[d, d] are loose, resulting in a large number of active nodes in the tree, which increases exponentially the cost of exploring a whole subtree. Furthermore, since other interval distances are also employed in the discretization, the sampled positions in the feasible arcs for previous atoms are only approximations for their true positions, and such a sequence of approximate positions may lead to an infeasibility at a further layer. For this reason, the longest-range pruning distances may fail to be verified (eveniftheyarerepresentedbyaninterval).

123Fig. 6

 6 Fig. 6 Available distances for the instances 1FJK (left) and 2E2F (right)

) 123 Journal: 10898 -

 12310898 JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small Author Proof u n c o r r e c t e d p r o o f J Glob Optim

  Fig 2). Torsion angles on the protein backbones give rise to the definition of interval distances, related to the minimal and maximal extension of such torsion angles. Distances between pairs of hydrogens are also included, as specified in the next subsection.123 Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small Author Proof u n c o r r e c t e d p r o o f J Glob Optim

  that the number of calls and CPU time were reduced in 9 out of 11 instances. It is also worth to mention that the pruning criteria based on peptide plane geometry and chirality helped the new version of iBP in reducing the number of calls in some instances and improving the global fold as well. Concerning the MDE, the original iBP seems to be more stable, although it fails to solve four of the instances (within the specified timeout). On the other hand, since the New 123 Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small u n c o r r e c t e d p r o o f J Glob Optim

123Fig. 9

 9 Fig. 9 Four incongruent realizations for 1AQR. The conformations in the bottom are reflections of the top ones. All four conformations differ by partial reflections

3 .

 3 a new pruning device that "relaxes" long-range distance constraints (Sect. 4.3)w h i c h allows us to obtain approximate realizations.Computational experiments on artificially generated instances showed the effectiveness of all above mentioned points, when integrated in the iBP algorithm. We are in fact able to find approximate realizations for protein backbones up to 64 amino acids in an affordable time, 123 Journal: 10898-JOGO Article No.: 0493 TYPESET DISK LE CP Disp.:2016/12/24 Pages: 21 Layout: Small

Table 1

 1 Numerical results on artificially generated instances from the PDB

	PDB ID Instance		iBP from [34]N e w iBP
		aa |V || E|	D	Calls	Time	MDE	D Calls	Time MDE
	2JMY	15	120	660	13 37,658	0.13	3e-06 5	2983	0.01	1e-16
	2KXA	24	177	973	10 215,669	0.92	5e-06 3	5064	0.01	6e-03
	1DSK	28	222 1210	14 31,309	0.13	4e-06 4	53,890	0.14	1e-06
	2PPZ	36	287 1522	9	2,372,242 11.34 2e-06 3	442,965	1.87	4e-08
	1AQR	40	310 1596	*	*	*	*	4	114,671	0.20	6e-03
	2ERL	40	324 1792	14 1,495,282 6.14	4e-06 3	10,410	0.03	1e-03
	2E2F	41	315 1716	*	*	*	*	3	19916	0.06	9e-03
	1FJK	52	417 2306	12 115,426	0.73	4e-06 4	925,090	3.07	2e-06
	2JWU	56	448 2416	*	*	*	*	4	226870	0.81	1e-02
	2KIQ	57	455 2452	20 1,217,945 12.79 6e-06 4	317,136	1.12	7e-04
	2LOW	64	497 2650	*	*	*	*	3 3,738,152	8.79	2e-07

Acknowledgements DG and CL are thankful to the Brazilian research agencies FAPESP and CNPq for partial financial support. LL was partially supported by the "Bip:Bip" project within the ANR "Investissement d'Avenir" program. AM thanks University of Rennes 1 for financial support. AM and DG also acknowledge Brittany Region (France) for partial financial support.

Once the first solution X 1 is found by the "New iBP",asetoffeasiblee xactdistances for the distances that were originally represented by intervals can be selected. This set of distances defines a DGP instance with exact distances which contains X 1 in its finite solution set. Moreover, by solving such an instance with the basic BP algorithm (for exact distances), we can compute all other feasible conformations that can be obtained from X 1 by partial reflections [START_REF] Lavor | The discretizable molecular distance geometry problem[END_REF][START_REF] Liberti | Euclidean distance geometry and applications[END_REF]. This procedure gets rid of the flexings 1 in the molecule, but only in this case the solution set is finite.

Applying this scheme to a modified 1AQR instance, where hydrogen distances between consecutive amino-acids were removed and the threshold was lowered to 4.5 Å, four incongruent conformations are obtained, as depicted in Fig. 9.

We claim that, even though interval distances pose some difficulties to the extension of this combinatorial approach, it is still possible to explore all the (discrete) conformational space obtained with discretization. Henceforth, we propose our New iBP as an exploratory tool to enumerate protein conformations that satisfy most of the given distance restraints, that can be further improved by local minimization procedures. 1 Continuous motions of part of the structure preserving all distance restraints. u n c o r r e c t e d p r o o f J Glob Optim and with reasonable precision, that can be further improved by using our solutions as starting points for a local minimization solver.

In the presented experiments, the two compared versions of the iBP algorithm were both used for identifying only one solution to the problem. However, as remarked before and illustrated in Sect. 5.3,theiBP algorithm can potentially enumerate the entire solution set of a discretizable MDGP. Research is currently focused on efficiently enumerating all conformations belonging to the search tree. Due to the presence of interval distances, many solutions may belong to the same cluster/ensemble of conformations. Hence, the next step is to define a method to classify the solutions in equivalence classes and integrate iBP with an scheme able to pick only one representative conformation from each incongruent ensemble.


Olivier Wang
email: olivier.wang@polytechnique.edu

Leo Liberti
email: leo.liberti@polytechnique.edu

Controlling some statistical properties of Business Rules programs

Keywords: Distribution learning, mixed-integer programming, statistical goals, business rules

Business Rules programs encode decision-making processes using "if-then" constructs in a way that is easy for non-programmers to manipulate. A common example is the process of automatic validation of a loan request for a bank. The decision process is defined by bank managers relying on the bank strategy and their own experience. Bank-side, such processes are often required to meet goals of a statistical nature, such as having at most some given percentage of rejected loans, or having the distribution of requests that are accepted, rejected, and flagged for examination by a bank manager be as uniform as possible. We propose a mathematical programming-based formulation for the cases where the goals involve constraining or comparing values from the quantized output distribution. We then examine a simulation for the specific goals of (1) a max percentage for a given output interval and (2) an almost uniform distribution of the quantized output. The proposed methodology rests on solving mathematical programs encoding a statistically supervised machine learning process where known labels are an encoding of the required distribution.

Introduction

Business Rules (BR) are a "programming for non programmers" paradigm that is often used by large corporations to store industrial process knowledge formally. BR replaces the two most abstract concepts of programming, namely loops and function calls, by means of an implicit outer loop and meta-variables used within a set of easy-to-manage "if-then" type instructions. BR interpreters are implemented by all BR management systems, e.g. [START_REF][END_REF]. BR programs are often used by corporations to encode their policies and empirical knowledge: given some technical input, they produce a decision, often in the form of a YES/NO output. Corporations often require their internal processes to perform according to a prescribed statistical behavior, which could be imposed because of strategy or by law. This required behavior is typically independent of the BR input data. The problem is then to parametrize the BR program so it will behave as prescribed on average, while still providing meaningful YES/NO answers on given inputs.

In [START_REF] Wang | Controlling the average behaviour of business rules programs[END_REF] we studied a simplified version of the problem where the statistical behavior was limited to a given mean. In this paper we provide a solution methodology for a more general (and difficult) case, where the statistical behavior is described by a given discrete distribution. We achieve this goal by encoding a Machine Learning (ML) procedure by means of a Mathematical Program (MP) of the Mixed-Integer Linear Programming (MILP) type. The ML procedure relies on non-input specific labels that encode the given knowledge about the distribution. Controlling the statistical behavior of a complex process such as a BR program is a very hard task, and to the best of our knowledge this work is the first of its kind in this respect. Methodologically speaking, we think our MILP formulation is also innovative in that it encodes an ML training process having labels which, instead of applying to individual inputs, apply to the entire input distribution at once. Such an ML process bypasses the usual difficulties of trying to label the training set data, thereby being more practical for industrial applications.

The motivation for this study is a real industrial need expressed by IBM (which co-funds this work) with respect to their BR package ODM. Our previous paper [START_REF] Wang | Controlling the average behaviour of business rules programs[END_REF] laid some of the groundwork, limited to the most basic statistical indicator (the mean of a distribution). Though that was a necessary step to the current work, the methodology described herein is the first to actually address the need expressed by industry: we feel this is one of the main feature that sets this work apart from our previous work. We still rely on MILP-based methodology, but now the input is a whole discrete distribution, the cardinality of which largely determines the size of the new MILP formulations presented below. Our tests show that this has an acceptable impact on empirical solution complexity.

As an experimental illustration, we consider the two cases of the statistical behavior being (1) a maximum percentage of a certain output value and (2) the output values being distributed in a fashion close to the uniform distribution, for integer outputs. We provide an optimization based approach to solving the learning problem for each of those cases, then examine some test results.

Preliminaries

We formally represent a BR program as an ordered list of sentences of the form:

if cond(p, x) then x ← act(p, x) end if
where p is a control parameter vector (with c components) which encodes a possible "tuning" of the program (e.g. thresholds which can be adjusted by the user), x ∈ X ⊆ R d is a variable vector representing intermediate and final stages of computation, cond is a boolean function, and act a function with values in X.

We call rule such a sentence, condition an expression cond(p, x) and action an instruction x ← act(p, x), which indicates a modification of the value of x. We write the final value of the variable x as x f = P (p, q), where P represents the BR program and q is an input parameter vector representing a problem instance and equal to the initial value of x. Although in general BR programs may have any type of output, we consider only integer outputs, since BR programs are mostly used to take discrete decisions. We remark that p, x are symbolic vectors (rather than numeric vectors) since their components are decision variables.

BR programs are executed in an external loop construct which is transparent to the user. Without getting into the details of BR semantics, the loop executes a single action from a BR whose condition is True at each iteration. Which BR is executed depends on a conflict resolution strategy with varying complexity. De Sainte-Marie et al. [START_REF] Paschke | RIF production rule dialect (second edition)[END_REF] describe typical operational semantics, including conflict resolution strategy, for industrial BR management systems. In this paper, the list of rules is ordered and the loop executes the first BR of the list with a condition evaluating to True at each iteration. The loop only terminates once every condition of the BRs is False. We proved in [START_REF] Wang | The learnability of business rules[END_REF] that there is a universal BR program which can simulate any Turing Machine (TM), which makes the BR language Turing-complete.

We consider the problem where the q ∈ Q are the past, known instances of the BR program, and the outputs P (p, q) of those instances are divided into N evenly sized intervals [H 0 , H 1 ], . . . , [H N -1 , H N ], forming a quantized output distribution. Denoting ν 1 (p), . . . , ν N (p) the number of outputs in these categories, we can formalize the problem as:

min p,x p -p 0 1 C (ν 1 (p), . . . , ν N (p)) (1)
where • 1 is the L1 norm and C is a constraint or set of constraints. While this formulation uses the number of outputs rather than the probabilities themselves, the relation between the two is simply a ratio of 1/m, where m = card(Q) is the number of training data points.

In this paper, we suppose that P 1 and P 2 are BR programs with a rule set {R r | r ≤ ρ} containing rules of the form:

if L r ≤ x ≤ G r then x ← A r x + B r end if with L r , G r , B r ∈ R
and A r ∈ {0, 1} d×d . We note R = {1, . . . , ρ} and D = {1, . . . , d}. We discuss the concrete example of banks using a BR program in order to decide whether to grant a loan to a customer or not. The BR program depends on a variable vector x and initializes its parameter vector (a component of which is an income level threshold) to p 0 . A BR program P 1 is used to decide whether a first bank will investigate the loan request further or simply accept the automated decision taken by an expert system, and therefore has a binary output value. This bank's high-level strategy requires that no more than 50% of loans are treated automatically, but P 1 currently treats 60%. Another bank instead uses a BR program P 2 to accept, reject, or assign a bank manager to the loan request, and therefore has a ternary return value, represented by an integer in {0, 1, 2}. That bank's strategy requires that the proportion of each output is {1/3, 1/3, 1/3}, but it is currently {1/4, 1/4, 1/2}. Our aim is in each case to adjust p, e.g. modifying the income level, so that the BR program satisfies the bank's goal regarding automatic loan treatment. This adjustment of parameters could be required after a change of internal or external conditions, for example.

The first scenario can be formulated as:

min p,x p -p 0 1 E q∈Q P 1 (p, q) ≤ g (2)
where P 1 has an output in {0, 1}, g ∈ [0, 1] is the desired max percentage of 1 outputs, the q ∈ Q are the past known instances of the BR program, • 1 is the L1 norm, p, q must satisfy the semantics of the BR program P (p, q) when executed within the loop of a BR interpreter and E is the usual notation for the expected value.

Similarly, the second scenario where P 2 has an output in {1, . . . , N } and the desired output is as close to a uniform distribution as possible can be formalized as:

min p,x p -p 0 1 ∀s, t ∈ {1, . . . , N }, ν s -ν t ≤ 1 (3)
Note that the solution to this problem is not always a truly uniform distribution, simply because there is no guarantee that m is divisible by N . However, it will always be as close as possible to a uniform distribution, since the constraint imposes that all the outputs will be reached by either floor(m/N ) or ceil(m/N ) data points. Again, we use whole numbers (of outputs in a given interval) instead of frequencies to be able to employ integer decision variables. Such problems could be solved heuristically by treating P 1 or P 2 as a blackbox, or by replacing it by means of a simplified model, such as e.g. a low-degree polynomial. We approach this problem as in [START_REF] Wang | Controlling the average behaviour of business rules programs[END_REF]: we model the algorithmic dynamics of the BR by means of MIP constraints, in view to solving those equations with an off-the-shelf solver. That this should be possible at all in full generality stems from the fact that Mathematical Programming (MP) is itself Turing-complete [START_REF] Liberti | Mathematical programming: Turing completeness and applications to software analysis[END_REF].

We make a number of simplifying assumptions in order to obtain a practically useful methodology, based on solving a Mixed-Integer Linear Programming (MILP) reformulation of these equations using a solver such as CPLEX [START_REF]IBM: ILOG CPLEX 12.2 User's Manual[END_REF]:

1. We suppose Q is small enough that solving the MILP is (relatively) computationally cheap. 2. We assume finite BR programs with a known bound (n -1) on the number of iterations of the loop for any input q (industrial BR programs often have a low value of n relative to the number of rules). This in turn implies that the values taken by x during the execution of the BR program are bounded.

We assume that M 1 is an upper bound of all absolute values of all p, q, and x, as well as any other values appearing in the BR program. It serves as a "big M " for the MP described in the rest of the paper.

3. We assume that the conditions and actions of the BR program give rise to constraints for which an exact MILP reformulation is possible. In order to have a linear model, each BR must thus be "linear", i.e. have the form:

if L ≤ x ≤ G then x ← Ax + B end if with L, G, B ∈ R d and A ∈ {0, 1} d×d . In general, A h,k may have values in R if it
is not a parameter and x h has only integer values.

Related Works

We follow the formalism used in [START_REF] Wang | Controlling the average behaviour of business rules programs[END_REF] pertaining to Business Rules (BR) programs and their statistical behavior.

Business Rules (also known as Production Rules) are well studied as a knowledge representation system [START_REF] Clancey | The epistemology of a rule-based expert system: a framework for explanation[END_REF][START_REF] Davis | Production rules as a representation for a knowledge-based consultation program[END_REF][START_REF] Lucas | Principles of Expert Systems[END_REF], originating as a psychological model of human behavior [START_REF] Newell | Production systems: Models of control structures[END_REF][START_REF] Newell | Human Problem Solving[END_REF]. They have further been used to encode expert systems, such as MYCIN [START_REF]Rule Based Expert Systems: The Mycin Experiments of the Stanford Heuristic Programming Project[END_REF], [START_REF] Shortcliffe | Computer-Based Medical Consultations: MYCIN[END_REF], EMYCIN [START_REF]Rule Based Expert Systems: The Mycin Experiments of the Stanford Heuristic Programming Project[END_REF][START_REF] Scott | The EMYCIN Manual[END_REF], OPS5 [START_REF] Brownston | Programming Expert Systems in OPS5: An Introduction to Rule-based Programming[END_REF][START_REF] Forgy | OPS5 User's Manual[END_REF], or more recently ODM [START_REF][END_REF] or OpenRules [START_REF] Openrules | OpenRules User Manual[END_REF]. On business side of things, they have been defined broadly and narrowly in many different ways [START_REF] Knolmayer | Business rules[END_REF][START_REF] Kolber | Defining business rules -what are they really? Project Report 3, The Business Rules Group[END_REF][START_REF] Ross | Principles of the Business Rule Approach[END_REF]. We consider Business Rules as a computational tool, which to the best of our knowledge has not been explored in depth before.

Supervised Learning is also a well studied field of Machine Learning, with many different formulations [START_REF] Bakir | Predicting Structured Data (Neural Information Processing)[END_REF][START_REF] Liu | Learning to rank for information retrieval[END_REF][START_REF] Settles | Active learning literature survey[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. A popular family of algorithms for the classification problem uses Association Rules [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF][START_REF] Malioutov | Exact rule learning via boolean compressed sensing[END_REF]. Such Rule Learning is not to be confused with the problem treated in this article, which is more a regression problem than a classification problem. There exist many other algorithms for Machine Learning, from simple linear regression to neural networks [START_REF] Atiya | Learning Algorithms for Neural Networks[END_REF] and support vector machines [START_REF] Cortes | Support-vector networks[END_REF]. When the learner does not have as many known output values as it has items in the training set, the problem is known as Semi-Supervised Learning [START_REF] Chapelle | Semi-Supervised Learning[END_REF]. Similarly, there has been research into machine learning when the matching of the known outputs values to the inputs is not certain [START_REF] Brodley | Identifying mislabeled training data[END_REF]. A previous paper has started to explore the Learning problem when the known information does not match to a single input [START_REF] Wang | Controlling the average behaviour of business rules programs[END_REF].

Learning Goals with Histograms

In the rest of this paper, we concatenate indices so that (L r

) k = L rk , (G r ) k = G rk , (A r ) h,k = A rhk and (B r ) k = B rk . We assume that rules are feasible, i.e. ∀r, k ∈ R × D, L k ≤ G k .
In the rest of this section, we suppose that the dimension of p is c = 1, making p a scalar, and that p takes the place of A 111 . Similar sets of constraints exists for when the parameter p takes the place of a scalar in B r , L r or G r . Additional parameters correspond to additional constraints that mirror the ones used for the first parameter.

This formalization is taken from [START_REF] Wang | Controlling the average behaviour of business rules programs[END_REF], in which we have also proved that the set of constraints described in Fig. 1 models the execution of such a BR program. The iterations of the execution loop are indexed by i ∈ I = {1, . . . , n} where n -1 is the upper bound on the number of iterations, the final value of x corresponds to iteration n. We use an auxiliary binary variable y ir with the property: y ir = 1 iff the rule R r is executed at iteration i. The other auxiliary binary variables y U ir and y L ir are used to enforce this property. We note (C1), (C2), etc. the constraints related to the evolution of the execution and (IC1), (IC2), etc. the constraints related to the initial conditions of the BR program:

-(C1): represents the evolution of the value of the variable x -(C2): represents the property that at most one rule is executed per iteration -(C3): represents the fact that a rule whose condition is False cannot be executed -(C4)-(C6) represent the fact that only the first rule whose condition is True can be executed -(IC1) through (IC3) represent the initial value of a -(IC4) represents the initial value of x.

∀i ∈ I\{n} x i+1 = r∈R (arx i + Br)yir + (1 - r∈R yir)x i (C1) ∀i ∈ I r∈R yir ≤ 1 (C2) ∀(i, r) ∈ I×R Lr -M (1 -yir)e ≤ x i ≤ Gr + M (1 -yir)e (C3) ∀(i, r, k) ∈ I×R×D x i k ≥ G rk -M y U irk -M r <r y ir (C4) ∀(i, r, k) ∈ I×R×D x i k ≤ L rk + M y L irk + M r <r y ir (C5) ∀(i, r) ∈ I×R 2d -1 + yir ≥ k∈D (y U irk + y L irk ) (C6)
∀r ∈ {2, . . . , ρ} ar = Ar (IC1)

a111 = p (IC2) ∀(h, k) ∈ D 2 \ {1, 1} a 1hk = A 1hk (IC3) x 1 = q (IC4) ∀i ∈ I x i ∈ X ∀r ∈ R ar ∈ {0, 1} d×d ∀(i, r, k) ∈ I×R×D yir, y U irk , y L irk ∈ {0, 1}
Fig. 1. Set of constraints modeling the execution of a BR program (e ∈ R d is the all-one vector).

A MIP for learning quantized distributions

The Mixed-Integer Program from Fig. 2 models the problem from Eq. 1. We index the instances in Q with j ∈ J = {1, . . . , m}. We also limit ourselves to solutions which result in computations that terminate in less than n -1 rule executions. As modifying the parameter means modifying the BR program, the assumptions made regarding the finiteness of the program might not be verified otherwise.

We note

O = {1, . . . , N }, such that ∀t ∈ O, ν t = card{j ∈ J | x 1 n,j ∈ [H t-1 , H t ]}.
We enforce this definition of ν t by using an auxiliary binary variable s tj with the property:

s tj = 1 iff x 1 n,j ∈ [B t-1 , B t ].
The other auxiliary binary variables s U tj and s L tj are used to enforce this property. The constraints are mostly similar to the ones in Fig. 1. We simply add the goal of minimizing the variation of the parameter value and the constraints C (ν 1 (p), . . . , ν N (p)) from Eq. 1. The new constraints are:

-(C7) represents the need for the computation to have terminated after n -1 executions -(C8)-(C12) represents the definition of ν 1 , . . . , ν N -(IC4') represents (IC4) with an additional index j.

That solving the MIP in Fig. 2 also solves the original Eq. 1 is a direct consequence of the fact that the constraints in Fig. 1 simulate P (p,q). The proof is simple since (C8) through (C12) trivially represent the definition of ν 1 , . . . , ν N . A similar MIP can be obtained when p has values in different part of the BRs, from which a more complex MILP is obtained for when p is non-scalar. However, this formulation is still quite abstract, as it depends heavily on the form of C . In fact, it can almost always be simplified given a particular constraint over the quantized distribution, as we see in the rest of this section.

A MILP for the Max Percentage Problem

A constraint programming formulation of Eq. 2 is the Mixed-Integer Linear Program (MILP) described in Fig. 3. In the case of the Max Percentage problem, we can linearize the MIP in Fig. 2 as well as remove some superfluous variables, since only one of the ν t is relevant.

We now note e = (1, . . . , 1) ∈ R d the vector of all ones. We use the auxiliary variables w ∈ R I×J×R and z ∈ R I×J×R×D 2 such that w ijr = (A r x ij + B rx i,j )y ijr (i.e. w ijr = A r x i,j + B r -x i,j , the difference between the new and the old values of x j ) and z ijrhk = a rhk x i,j k . Any constraints numbered as before fulfills the same role. The additional constraints are:

-(C1' 1 ), (C1' 2 ), (C1' 3 ), (C1' 4 ) and (C1' 5 ) represent the linearization of (C1) from Fig. 1 -(C8') represents the goal from Eq. 2, that is a constraint over the average of the final values of x. It replaces C (ν 1 , . . . , ν N ) and all the constraints used to define ν t from the MIP in Fig. 2. The MILP from Fig. 3 finds a value of p that satisfies Eq. 2. This is again derived from the fact that Fig. 1 simulates a BR program, and from the trivial proof that (C1' 1 ), (C1' 2 ), (C1' 3 ), (C1' 4 ) and (C1' 5 ) represent the linearization of (C1).

∀(t, j) ∈ O×J Ht-1 -M (1 -stj) ≤ x n,j 1 ≤ Ht + M (1 -stj) (C8) ∀(t, j) ∈ O×J x n,j 1 ≥ Ht -M s U tj (C9) ∀(t, j) ∈ O×J x n,j 1 ≤ Ht-1 + M s L tj (C10) ∀(t, j) ∈ O×J stj ≥ s U tj + s L tj (C11) ∀t ∈ O νt = j∈J stj (C12) ∀j ∈ J x 1,j = q j (IC4') ∀(i, j) ∈ I×J x i,j ∈ X ∀k ∈ R ak ∈ {0, 1} d×d p ∈ {0, 1} ∀(i, j, r, k) ∈ I×J×R×D yijr, y U ijrk , y L ijrk ∈ {0, 1} ∀(t, j) ∈ O×J stj, s U tj , s L tj ∈ {0, 1} ∀t ∈ O νt ∈ N

A MILP for the Almost Uniform Distribution Problem

As before, we exhibit in Fig. 4 a MILP that solves Eq. 3. Any constraints numbered as before fulfills the same role. The additional constraints are:

-(C8") through (C10") represent the adaptation of (C8) through (C10) to the relevant case of integer outputs -(C13) represents the equivalent to C from Eq. 3.

This MILP is obviously equivalent to solving Eq. 3, since it is for the most part a straight linearization of the MIP in Fig. 2.

Implementation and Experiments

We use a Python script to randomly generate samples of 100 instances of P 1 and P 2 for different numbers of control parameters c, each instance having a minimize p,a,x,y,y U ,y L ,w,z p 0 -p subject to (C2), (C3), (C4), (C5), (C6), (C7), (IC1), (IC2), (IC3)

∀(i, j) ∈ I\{n}×J x i+1,j = r∈R wijr + x i,j (C1' 1 ) ∀(i, j) ∈ I×J×R -M yijre ≤ wijr ≤ M yijre (C1' 2 ) ∀(i, j, r, h) ∈ I×J×R×D k∈D zijrhk + Brh -x i,j h -M (1 -yijr) ≤ wijrh ≤ k∈D zijrhk + Brh (C1' 3 ) -x i,j + M (1 -yijr)e ∀(i, j, r) ∈ I×J×R -M ar ≤ zijr ≤ M ar (C1' 4 ) ∀(i, j, r, h, k) ∈ I×J×R×D 2 x i,j k -M (1 -arhk) ≤ zijhk ≤ x ij k (C1' 5 ) ∀j ∈ J r∈R ynjr = 0 (C7) j∈J x n,j 1 ≤ mg (C8') ∀(i, j, r, k1, k2) ∈ I×J×R×D 2 x i,j , zijrhk ∈ X ∀(i, j, r) ∈ I×J×R wijr ∈ R d ∀r ∈ R ar ∈ {0, 1} d×d p ∈ {0, 1} ∀(i, j, r, k) ∈ I×J×R×D yijr, y U ijrk , y L ijrk ∈ {0, 1}
Fig. 3. MILP formulation for solving Eq. 2.

corresponding set of inputs with d = 3, n = 10 and m = 100. The number of control parameters serves as an approximation of the complexity of the BR program to optimize: a more complex program will have more buttons to adjust, thus increasing the complexity, yet be more likely to have the goal be reachable at all, i.e. have the MILP be feasible. We define the space X as X ⊆ R × R × Z.

The BR programs are sets of ρ = 10 rules, where L r , G r , B r are vectors of scalars in an interval range and A r are d × d matrices of binary variables. In P 1 , we use range = [0, 1] and in P 2 , we use range = [0, 3]. All input values q are generated using a uniform distribution in range.

We use these BR programs to study the computational properties of the MILP. The value of M used is customized according to each constraint, and is ultimately bounded by 6 and 16 in P 1 and P 2 respectively (strictly greater than five times the range of possible values for x). We write the MILP as an AMPL model, and solve it using the CPLEX solver on a Dell PowerEdge 860 running CentOS Linux. minimize p,b,x,y,y U ,y L ,w,s,s L ,s g ,ν p 0 -p subject to (C1' 1 ), (C1' 2 ), (C1' 3 ), (C1' 4 ), (C1' 5 ), (C2), (C3), (C4), (C5), (C6), (C7), (C11), (C12) (IC1), (IC2), (IC3), (IC4')

∀(t, j) ∈ O×J t -M (1 -stj) ≤ x n,j 1 ≤ t + M (1 -stj) (C8") ∀(t, j) ∈ O×J x n,j 1 ≥ t -M s U tj (C9") ∀(t, j) ∈ O×J x n,j 1 ≤ t + M s L tj (C10") ∀(t, τ ) ∈ O 2 -1 ≤ νt -ντ ≤ 1 (C13) ∀(i, j) ∈ I×J x i,j ∈ X ∀r ∈ R ar ∈ {0, 1} d×d p ∈ {0, 1} ∀(i, j, r, k) ∈ I×J×R×D yijr, y U ijrk , y L ijrk ∈ {0, 1} ∀(t, j) ∈ O×J stj, s U tj , s L tj ∈ {0, 1} ∀t ∈ O νt ∈ N Fig. 4
. MILP formulation for solving Eq. 3.

The Max Percentage Problem

We observe the proportion of solvable instances of P 1 for c between 5 and 10 and c = 15 in Tab. 1. We use the MILP in Fig. 3 to solve Eq. 2 with the goal set to g = 0.5. An instance is considered solvable if CPLEX reports an integer optimal solution or a (non-)integer optimal solution. We separate the instances where the optimal value is 0 from the others, as those indicate that the randomly generated BR program already fulfill the goal condition. We expect around fifty of those for any value of c.

In Fig. 5, we observe both the success rate and the average solving time when considering only the non-trivial, non-timed out instances of P 1 . The success rate increases steadily, as expected. The solving time seems to indicate a non-linear increase for c greater than 6, even with its values being somewhat unreliable due to the small sample. Knowing that average industrial BRs are more complex than our toy examples, regularly having thousands of rules, this approach to the Maximum Percentage problem does not seem applicable to industrial cases.

The Almost Uniform Distribution Problem

We observe the proportion of solvable instances of P 2 for c between 5 and 10 and c = 15 in Tab. 2. We use the MILP in Fig. 4 we separate instances where the goal is already achieved before optimization, identifiable by being solved quickly with a value of p = p 0 , i.e. an optimal value of zero. In Fig. 6, we display the success rate and average solving time over the non-timed out, non-presolved instances for all three values of c. We observe a sharply non-linear progression, with the average problem taking about nine minutes with 15 control parameters. Knowing that average industrial BRs are much more complex than our toy examples, regularly having thousands of rules, we conclude that this method can only be used infrequently, if at all.

Conclusion, Discussion and Future Work

We have presented a learning problem of unusual type, that of supervised learning with statistical labels. We have further explored a particular subset of those problems, those where the labels apply to a quantized output distribution. This new approach is easily applied to practical applications in industry where control parameters must be learned to satisfy a given goal. We have given a mathematical programming algorithm that solves such a learning problem given a linear BR program. Depending on the specific learning problem, the mathematical program might be easy or difficult to solve. We examined two example learning problems with practical applications for which the learning is equivalent to solving a MILP.

We observe that, though one could detect a visual similarity in the plots presented in Fig. 5 and6, we believe that this similarity is only apparent. In fact, the error bars (which measure the standard deviation of the solution time over the instance subclass corresponding to a given size of parameters) point out that the "hard cases" (with high values of solution times) are also the cases where the error bars are longest. In other words, this "similarity" is simply a result of outliers in the corresponding peaks.

The experimental results indicate the general feasibility of this type of approach. It is clear that, due to the exponential nature of Branch-and-Bound (BB, the algorithm solving the MILPs), the performance will scale up poorly with the size of the BR program: but this can currently be said of most MILPs. This issue, which certainly requires more work, can possibly be tackled by pursuing some of the following ideas: more effective BB-based or formulation-based heuristics (also called mat-heuristics in the literature), cut generation based on problem structure, and decomposition. The latter, specifically, looks promising as the structure of the BR program is, up to the extent provided by automatic translation based on parsing trees, carried over to the resulting MILP.

Other avenues of research are in extending this statistical learning approach in other directions, e.g. learning other moments, or given quantiles in continuous distributions. Statistical goal learning problems are an apparently unexplored area of ML that has eminently practical applications.

  minimize p,a,x,y,y U ,y L ,s,s U ,s L ,ν p 0 -p subject to (C1), (C2), (C3), (C4), (C5), (C6), (IC1), (IC2), (IC3) C (ν1, . . . , νN )

Fig. 2 .

 2 Fig. 2. Mixed-Integer Program solving Eq. 1.

Fig. 5 .

 5 Fig. 5. Average solution time over P1 for varying values of c in seconds.

Fig. 6 .

 6 Fig. 6. Average solution time over non-trivial solvable P2 for varying values of c.

Table 1 .

 1 to solve Eq. 3 with N = 2. Again, Experimental values for the maximum percentage problem.

	Number of control parameters c	5	6	7	8	9	10 15
	Trivial solvable instances (objective = 0)	52 53 49 49 58 48 46
	Non-trivial solvable instances (objective = 0)	5	6	5	13	6	6	8
	Infeasible instances	43 43 40 36 31 35 14
	Timed out instances	0	0	7	2	5	11 32

Table 2 .

 2 Experimental Values for the Almost Uniform Distribution problem.

	Number of control parameters c	5	6	7	8	9	10 15
	Trivial solvable instances (objective = 0)	8	2	1	1	4	7	4
	Non-trivial solvable instances (objective = 0)	9	2	8	5	4	15 32
	Infeasible instances	83 96 91 93 92 77 63
	Timed out instances	0	0	0	1	0	1	1

Acknowledgments. The first author (OW) is supported by an IBM France/ANRT CIFRE Ph.D. thesis award.