Pierre Montalbano
email: pierre.montalbano@inrae.fr

Simon De Givry
email: simon.de-givry@inrae.fr

George Katsirelos

Multiple-choice knapsack constraint in graphical models

Keywords: graphical model, cost function network, knapsack problem

Graphical models, such as cost function networks (CFNs), can compactly express large decomposable functions, which leads to efficient inference algorithms. Most methods for computing lower bounds in Branch-and-Bound minimization compute feasible dual solutions of a specific linear relaxation. These methods are more effective than solving the linear relaxation exactly, with better worst-case time complexity and better performance in practice. However, these algorithms are specialized to the structure of the linear relaxation of a CFN and cannot, for example, deal with constraints that cannot be expressed in extension, such as linear constraints of large arity. In this work, we show how to extend soft local consistencies, a set of approximate inference techniques for CFNs, so that they handle linear constraints, as well as combinations of linear constraints with at-mostone constraints. We embedded the resulting algorithm in toulbar2, an exact Branch-and-Bound solver for CFNs which has demonstrated superior results in several graphical model competitions and is state-of-theart for solving large computational protein design (CPD) problems. We significantly improved performance of the solver in CPD with diversity guarantees. It also compared favorably with integer linear programming solvers on knapsack problems with conflict graphs.

Introduction

A Graphical Model (GM) may express an arbitrary complex function on several variables as a combination of smaller local functions on subsets of the variables. GMs have been used to reason about logic and probabilities. A deterministic GM can represent a Constraint Satisfaction Problem (CSP) where each local function is a constraint evaluating to true (satisfied) or false (unsatisfied) and the combination operator is Boolean conjunction. It can also represent a Cost Function Network (CFN), where each function evaluates to a cost and the combination operator is addition [START_REF] Cooper | Valued Constraint Satisfaction Problems[END_REF]. A probabilistic GM represents a probability distribution on random variables. Local functions may correspond to conditional probability distributions as in Bayesian Networks (BN) or potentials as in Markov Random Fields (MRF) [START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF]. They are combined by multiplication. In the following, we focus on CFNs. It can be shown that finding the Maximum A Posteriori assignment on MRFs (MAP/MRF) or the Most Probable Explanation on BNs (MPE/BN) can be cast as finding a solution of minimum cost in an appropriate CFN [START_REF] Cooper | Graphical models: Queries, complexity, algorithms (tutorial)[END_REF]. By allowing infinite costs to represent infeasibility, CFNs can be seen as a strict generalization of CSPs.

Exact methods to solve GMs/CFNs mostly rely on Branch-and-Bound (B&B) algorithms [START_REF] Marinescu | And/or branch-and-bound for graphical models[END_REF][START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF]. These methods have proved useful in many GM applications, such as resource allocation [START_REF] Cabon | Radio Link Frequency Assignment[END_REF], image analysis [START_REF] Haller | Exact map-inference by confining combinatorial search with LP relaxation[END_REF], or computational biology [START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Allouche | Operations Research and Simulation in Healthcare, chap. Cost Function Networks to Solve Large Computational Protein Design Problems[END_REF]. For those, it has been shown to outperform other approaches, including Integer Linear Programming (ILP), MaxSAT and Constraint Programming (CP) [START_REF] Hurley | Multi-language evaluation of exact solvers in graphical model discrete optimization[END_REF].

CFNs have no native way to express linear constraints. This is in large part due to the algorithms used to compute lower bounds in B&B, which require that all constraints be expressed in extension. In many cases, having the ability to add such constraints would significantly improve the usefulness of CFNs. For example, when searching for diverse solutions, the Hamming distance constraint is naturally expressed as a linear constraint. There are ways to work around this [START_REF] Ruffini | Guaranteed diversity & quality for the weighted CSP[END_REF][START_REF] Ruffini | Guaranteed diversity and optimality in cost function network based computational protein design methods[END_REF] but, as we show later, they come with a non-trivial performance penalty. The lack of linear constraints is even more severe when the constraints have large coefficients, as in the knapsack problem with a conflict graph (KPCG). In this case, there is no workaround for the lack of linear constraints and CFN technology cannot be applied.

Contributions. Here, we show how to extend soft local consistency algorithms, a set of approximate inference techniques for CFNs, to deal with Pseudo-Boolean linear constraints (PB constraints for short), i.e., linear constraints over 0/1 variables. In the presence of unary cost functions (a cost function coupling a cost to each value), these correspond to knapsack constraints. We additionally consider the combination of PB constraints with Exactly-One (EO) or At-Most-One (AMO) constraints, which correspond to multiple-choice knapsack constraints [START_REF] Pisinger | Knapsack problems[END_REF], allowing finite-domain variables.

This new ability enables more modeling options for GM/CFN users, which we demonstrate by applying it to generating diverse solutions for Computational Protein Design (CPD). Here, the objective function has quadratic and linear terms that can be decomposed in a sum of binary cost functions on pairs of variables (the quadratic terms) and unary cost functions on single variables (linear terms). Searching for diverse solutions introduces linear constraints in the model (see Section 5.1). Our approach also compared favorably with a stateof-the-art ILP solver on knapsack problems with conflict graphs.

Related work

Problems defined by PB constraints are a generalization of the SAT problem. Solvers for PB SAT typically use SAT-inspired constraint learning techniques, either by direct translation to Conjunctive Normal Form (CNF) [START_REF] Sakai | Construction of an ROBDD for a PB-constraint in band form and related techniques for PB-solvers[END_REF] or by generalizing the clause learning mechanism to PB constraints [START_REF] Elffers | Divide and conquer: Towards faster pseudo-boolean solving[END_REF]. These solvers typically do not compute lower bounds during search and have to rely on conflict reasoning only to prove bounds. A notable exception is RoundingSAT [START_REF] Devriendt | Learn to relax: Integrating 0-1 integer linear programming with pseudo-boolean conflict-driven search[END_REF], which uses a Linear Programming (LP) solver to compute bounds during search and learn constraints from bound violations, but limits the number of iterations given to the LP solver in order to keep the runtime overhead of the solver reasonable. This is in contrast to our approach, which uses a suboptimal LP solver, but places no resource bounds on it. Also, PB solvers are usually restricted to a linear objective, whereas our approach can combine PB constraints with non-linear quadratic (or more) cost functions. PB solvers can also exploit the presence of AMO or EO constraints to strengthen propagation of PB constraints [START_REF] Ansótegui | Automatic detection of at-most-one and exactly-one relations for improved SAT encodings of pseudo-boolean constraints[END_REF][START_REF] Bofill | An mdd-based SAT encoding for pseudoboolean constraints with at-most-one relations[END_REF].

ILP solvers are well suited to solve CFNs, given the local polytope. Their LP solving is not limited to a specific form of LP, like soft local consistency algorithms such as Existential Directional Arc Consistency (EDAC) [START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF] and Virtual AC (VAC) [START_REF] Cooper | Soft arc consistency revisited[END_REF] are, therefore they have no issue reasoning with other linear constraints, as well as combinations with AMO/EO constraints. However, previous evaluations [START_REF] Hurley | Multi-language evaluation of exact solvers in graphical model discrete optimization[END_REF] showed that the size of the linear program that specifies that local polytope is often too large even for such highly optimized implementations and therefore they perform worse than a dedicated CFN solver in such problems.

On the CFN side, there has been work on clique constraints [START_REF] De Givry | Clique cuts in weighted constraint satisfaction[END_REF], a special case of PB constraints. Dlask and Werner [START_REF] Dlask | Bounding linear programs by constraint propagation: Application to Max-SAT[END_REF][START_REF] Dlask | On relation between constraint propagation and blockcoordinate descent in linear programs[END_REF] have shown how to handle arbitrary LPs using BCA algorithms, based on a generalization of VAC. However, despite recent advances [START_REF] Trösser | Relaxation-aware heuristics for exact optimization in graphical models[END_REF], BCA algorithms remain too costly for use at every node of a B&B. Many (soft) global constraints can be described by a set of (soft) linear constraints, but require an LP solver [START_REF] Lee | Consistency techniques for polytime linear global cost functions in weighted constraint satisfaction[END_REF]. In addition, maintaining (weak) EDAC and the coupling with the other local cost functions can be costly in practice. This was also the case for other soft global constraints exploiting flow-based or dynamic programming algorithms [START_REF] Lee | Consistency techniques for flow-based projection-safe global cost functions in weighted constraint satisfaction[END_REF][START_REF] Lee | Polynomially decomposable global cost functions in weighted constraint satisfaction[END_REF]. In our approach, we propose a simple and effective soft local consistency called Full ∅-Inverse Consistency for PB constraints. Finally, we can decompose linear constraints using cost functions of arity 3 and intermediate variables [START_REF] Allouche | Filtering decomposable global cost functions[END_REF], similar to CNF encodings used by PB solvers. However, the size of the domains of the intermediate variables increases linearly with the value of the coefficients of the PB constraints.

3 Preliminaries Definition 1. A Cost Function Network (CFN) P is a tuple (X, D, C, ) where X is a set of variables, with finite domain D(x) for x ∈ X. C is a set of constraints. Each constraint c ∈ C is defined over a subset of variables called its scope (scope(c) ⊆ X). is a maximum cost indicating a forbidden assignment.

The size of the scope of a constraint is its arity. Unary (resp. binary) cost functions have arity 1 (resp. 2). A partial assignment τ is an assignment of all the variables x i in its scope (scope(τ )) to a value of its domain D(x i ). The set of all the partial assignments on a scope S is denoted τ (S). A constraint over a scope S is denoted c S . The cost of a partial assignment τ for a constraint c S is denoted c S (τ ) with S ⊆ scope(τ ). Without loss of generality, we assume all costs are positive integers, bounded by , a special constant signifying infeasibility. Hence if c S (τ ) = then the assignment τ is not a feasible solution. A constraint c S is hard if for all τ ∈ τ (S), c S (τ ) ∈ {0, }, otherwise it is soft. A CFN P that contains only hard constraints is a constraint network (CN). In the following, we use the term cost function interchangeably with the term constraint. An assignment τ with scope(τ ) = X is a complete assignment. The cost of a complete assignment τ is given by c P (τ ) = c S ∈C c S (τ ). The Weighted Constraint Satisfaction Problem (WCSP) asks, given a CFN P , to find a complete assignment minimizing c P (τ ). This task is NP-hard [START_REF] Cooper | Valued Constraint Satisfaction Problems[END_REF]. When the underlying CFN is a CN, the problem is the CSP, which we call crisp CSP here. In the following, we use WCSP to refer both to the optimization task and the underlying CFN.

In this paper, we assume there exists exactly one unary constraint for each variable and we say that the unary cost of x i = v for some v ∈ D(x i ) is c i (v). We also assume the existence of a constraint c ∅ with empty scope, which represents a constant in the objective function and, since there exist no negative costs, it is a lower bound on the cost of all possible assignments.

Exact methods to solve GMs/CFNs mostly rely on Branch-and-Bound (B&B) algorithms [START_REF] Marinescu | And/or branch-and-bound for graphical models[END_REF][START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF]. At every node of the B&B tree, the solver computes a bound and closes the node if that bound is higher than the cost of the incumbent solution or if it represents infeasibility. Typical bounding algorithms compute either static memory-intensive bounds [START_REF] Dechter | Mini-buckets: A general scheme for bounded inference[END_REF] or memory-light ones [START_REF] Cooper | Soft arc consistency revisited[END_REF] better suited to dynamic variable orderings. The latter, on which we focus here, are called Soft Arc Consistencies (SAC) because they reason on each non-unary cost function one by one, in a generalization of propagation in CSP.

Soft arc consistencies use c ∅ as the lower bound and compute a reparameterization of the instance with a higher c ∅ . A reparameterization P of a WCSP P is a WCSP with an identical structure, i.e., one where there exist constraints over the same scopes, the costs assigned by each individual cost function may differ, but c P (τ ) = c P (τ ) for all complete assignments τ .

Procedure MoveCost(c S1 , c S2 , τ 1 , α): Move α units of cost between the tuple τ 1 of scope S 1 and tuples τ 2 that extend τ 1 in scope S 2 Data:

Scopes S1 ⊂ S2 Data: τ1 ∈ τ (S1) Data: cost α to move 1 c S 1 (τ1) ← c S 1 (τ1) + α 2 foreach τ2 ∈ τ (S2) | τ2[S1] = τ1 do 3 c S 2 (τ2) ← c S 2 (τ2) -α
All reparameterizations that we study here are computed as a sequence of local Equivalence Preserving Transformations (EPTs). Let S 1 ⊂ S 2 be two scopes with corresponding cost functions c S1 and c S2 . Procedure MoveCost describes how a cost α moves between the corresponding cost functions. To see its correctness, observe if τ 1 is used in a complete assignment, then exactly one extension of τ 1 to S 2 will be used. Therefore, the sum of c S1 and c S2 remains unaffected whether the cost α is attributed to τ 1 in c S1 or to all of its extensions τ 2 in c S2 . As an example, it is clear that adding a cost α on c x (a) and subtracting a cost α on c {x,y} ({x = a, y = b}) for all b ∈ D(y) preserves problem equivalence. Indeed, paying α when we assign x = a (cost function c x (a) = α) or when we assign x = a and y = b (∀b ∈ D(y)) (cost function c {x,y} ({x = a, y = b}) = α, ∀b ∈ D(y)) is equivalent. As a matter of terminology, when α > 0, cost moves from the larger arity cost function c S2 to the smaller arity c S1 and the move is called a projection, denoted project(c S1 , c S2 , τ 1 , α). When α < 0, cost moves to the larger arity cost function c S2 and the move is called an extension, denoted

extend(c S1 , τ 1 , c S2 , -α), equivalent to MoveCost(c S1 , c S2 , τ 1 , α). When S 1 = ∅ and |S 2 | = 1, with S 2 = {x i }, the move is called a unary projection, denoted unaryP roject(c i , α), equivalent to MoveCost(c ∅ , c i , ∅, α).
We never perform extensions from c ∅ , so it monotonically increases during the run of an algorithm and as we descend a branch of the search tree.

Finding which cost moves lead to an optimal reparameterization, which means one that derives the optimal increase in the lower bound, is not obvious. It has been shown that any reparameterization can be derived by a set of local cost moves [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF] and that the optimal reparameterization (with α rational)and, equivalently, the optimal set of cost moves -can be found from the optimal dual solution of the following linear relaxation of the WCSP [START_REF] Cooper | Soft arc consistency revisited[END_REF], whose feasible region is called the local polytope:

min c S ∈C,τ ∈τ (S) c S (τ ) × y τ s.t. y τ1 = τ2∈τ (S2),τ2[S1]=τ1 y τ2 ∀c S1 , c S2 ∈ C, S 1 ⊂ S 2 , τ 1 ∈ τ (S 1 ), |S 1 | ≥ 1 τ ∈τ (S) y τ = 1 ∀c S ∈ C, |S| ≥ 1
However, solving this LP to optimality is often prohibitively expensive because the worst-case complexity of an exact LP algorithm is O(N 2.5 ) [START_REF] Vaidya | Speeding-up linear programming using fast matrix multiplication[END_REF], with N ∈ O(ed + nd) for binary WCSPs, where e is the number of distinct binary cost functions, n is the number of WCSP variables and d is the maximum domain size. The poor asymptotic complexity matches empirical observation [START_REF] Hurley | Multi-language evaluation of exact solvers in graphical model discrete optimization[END_REF]. Moreover, the particular structure of this LP does not allow for a more efficient solving algorithm, as it has been shown that solving LPs of this form is as hard as solving any LPs [START_REF] Prusa | Universality of the local marginal polytope[END_REF]. Instead, work has focused on producing good but poten-tially suboptimal feasible dual solutions. Various algorithms have been proposed for this, going all the way back to Schlesinger [START_REF] Schlesinger | Sintaksicheskiy analiz dvumernykh zritelnikh signalov v usloviyakh pomekh (Syntactic analysis of two-dimensional visual signals in noisy conditions)[END_REF], who first expressed the problem as linear optimization and gave a specific algorithm for optimizing the dual. Since Schlesinger, a long line of algorithms has been pursued both in areas like image analysis [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF][START_REF] Werner | A Linear Programming Approach to Max-sum Problem: A Review[END_REF][START_REF] Sontag | Tightening LP relaxations for MAP using message-passing[END_REF][START_REF] Komodakis | MRF energy minimization and beyond via dual decomposition[END_REF][START_REF] Sontag | Efficiently searching for frustrated cycles in MAP inference[END_REF][START_REF] Tourani | Taxonomy of dual block-coordinate ascent methods for discrete energy minimization[END_REF], where Block-Coordinate Ascent (BCA) algorithms were developed, and constraint programming [START_REF] Schiex | Arc consistency for soft constraints[END_REF][START_REF] Larrosa | On arc and node consistency in weighted CSP[END_REF][START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF][START_REF] Zytnicki | Bounds Arc Consistency for Weighted CSPs[END_REF][START_REF] Cooper | Soft arc consistency revisited[END_REF], where they are called soft local consistencies. Notably, the strongest algorithms from both lines of research, such as TRWS [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF] and VAC [START_REF] Cooper | Soft arc consistency revisited[END_REF] converge on fixpoints with the same properties.

We do not describe all the existing local consistency algorithms but we need the following consistency properties: Definition 2. A WCSP P is Node Consistent (NC) [START_REF] Larrosa | On arc and node consistency in weighted CSP[END_REF] if for every variable x i ∈ X there exists a value v ∈ D(x i ) such that c i (v) = 0 and for every value

v ∈ D(x i ), c ∅ + c i (v ) < .
In the following, we assume that a WCSP is NC before our propagator runs. Definition 3. A WCSP P is ∅-Inverse Consistent (∅IC) [START_REF] Zytnicki | Bounds Arc Consistency for Weighted CSPs[END_REF] if for every cost function c S ∈ C there exists a tuple τ ∈ τ (S) such that c S (τ ) = 0. Definition 4. A WCSP P is Existential Arc Consistent (EAC) [START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF] if it is NC and for every x i ∈ X there exists a value v ∈ D(x i ) such that c i (v) = 0 and for every cost function c S ∈ C,

x i ∈ S, |S| > 1, there exists a tuple τ ∈ τ (S) verifying τ [x i ] = {v} (i.e., x i = v in τ ) and c S (τ ) + xj ∈S c j (τ [x j ]) = 0. Value v is called an EAC support.
This last definition applies only to binary cost function networks. 3 A weaker notion of EAC has been defined on global cost functions in order to avoid cost oscillation [START_REF] Lee | Consistency techniques for flow-based projection-safe global cost functions in weighted constraint satisfaction[END_REF]. Given a variable x i , it relies on a partition of the unary cost functions c j (τ [x j ]), x j ∈ X such that each part is associated to some non-unary cost function c S related to x i (x i ∈ S).

We follow another weakening approach related to ∅IC. We strengthen the previous definition to take into account unary costs as in EAC.

Definition 5. A WCSP is Full ∅-Inverse Consistent (F∅IC) if for every cost function c S ∈ C there exists τ ∈ τ (S) such that c S (τ ) + xj ∈S c j (τ [x j ]) = 0.
Compared to existing notions of consistency, F∅IC is weaker than T-DAC [START_REF] Allouche | Tractability-preserving transformations of global cost functions[END_REF]. It is also weaker than EAC on binary cost function networks, but it is incomparable with weak EAC [START_REF] Lee | Consistency techniques for flow-based projection-safe global cost functions in weighted constraint satisfaction[END_REF] on non-binary networks. Let ∀u ∈ D(x), α u = min v∈D(y) (c {x,y} ({x = u, y = v}) + c y (v)) and ∀v ∈ D(y), β v = max u∈D(x) (α u -c {x,y} ({x = u, y = v})). We apply project(c x , c {x,y} , {x = u}, α u ) for each value u ∈ D(x) and extend(c y , {y = v}, c {x,y} , β v ) for each value v ∈ D(y). These cost moves will result in adding a cost β v -α u to every tuple in c {x,y} . We have α a = 0, α b = α c = 1 and β a = 1, β b = β c = 0. All the costs remain positive (proof in [START_REF] Larrosa | In the quest of the best form of local consistency for weighted CSP[END_REF]). The reparameterized cost functions are c x (a) = c x (b) = c x (c) = 1, c {x,y} ({x = a, y = a}) = c {x,y} ({x = c, y = c}) = 1, the rest being equal to 0. We can now increase c ∅ by 1 using unaryP roject(c x , 1). The resulting WCSP is EAC.

For each of the consistencies we defined above, there exist corresponding algorithms that compute parametrization that satisfy them in polynomial-time 4 . Given the connection to linear programming, these reparameterizations map to feasible dual solutions of the local polytope. However, these algorithms rely on all constraints being expressed in extension, meaning that for all constraints the cost of every partial assignment must be explicitly written. This is not the case for many constraints that are typically used in modeling in CP, namely global constraints, i.e., those whose definition does not imply a fixed arity. In order to enforce these soft local consistencies in instances that contain global constraints, we need to define bespoke algorithms. In contrast with crisp CSPs, these algorithms must do more than prune values that appear in no feasible solution. They must compute a reparameterization such that the constraint satisfies the appropriate consistency, F∅IC here.

Here, we deal with pseudo-Boolean (PB) linear constraints and their generalizations. These are constraints of the form xi∈S w i x i C, where S is a scope, all x i ∈ S are Boolean variables, w i and C are constants and ∈ {<, ≤, =, ≥, >}. A PB constraint is normalized if w i , C ≥ 0 and is ≥. Any PB constraint can be written as a combination of normalized PB constraints. It is possible to detect in linear time whether this constraint is satisfiable in a crisp CSP, by testing if xi∈S max(0, w i ) ≥ C. It is also possible to detect values that appear in no solutions by computing all partial sums of |S| -1 variables, in linear time.

A PB constraint is an at-most-one (AMO) constraint if it has the form

xi∈S x i ≤ 1, normalized as xi∈S -x i ≥ -1.
It is an exactly-one (EO) constraint if it has the form xi∈S x i = 1.

Pseudo-Boolean constraints in CFNs

The specific constraint we consider here is a pseudo-Boolean constraint xi∈S w i x i ≥ C along with a partition of its variables into sets A 1 , . . . , A k such that there exists an EO constraint among the variables of each partition A i .

Reformulations. This formulation allows us to express PB constraints over multi-valued variables. Let S be a scope over a set of WCSP variables with arbitrary domains, and w iv weights for each value. The constraint xi∈S,v∈D(xi) w iv x iv ≥ C, where x iv is the 0/1 variable which takes the value 1 if x i = v, matches the pattern described above, with partitions

A i = {x iv | v ∈ D(X i )}.
Finally, this formulation admits the case where there exists an AMO constraint over some partitions: we add another 0/1 variable in each such partition and give it weight 0, so that this partition now has an EO constraint.

Constraint representation. We will focus here on F ∅IC as the soft consistency we aim to enforce. But first, we need an appropriate encoding that can represent the state of the constraint after a series of cost moves to and from unary cost functions, without storing a cost for each of the exponentially (in the arity of the constraint) many tuples. Observe first that the cost of any given tuple starts out at 0 for allowed tuples and for tuples that violate the constraint. After some cost moves, the cost of each tuple is the sum of costs that have been moved to or from the values it contains. Therefore, it can be expressed as a linear function. Let δ iv be the total cost that has been moved between the constraint and the corresponding unary cost and δ ∅ the cost we have moved from this constraint to c ∅ . Therefore, initially δ ∅ = 0 and δ iv = 0 for all i, v. We use the following integer program as the representation of the constraint.

min xi∈S,v∈D(xi) δ iv x iv -δ ∅ (1)
s.t. xi∈S,v∈D(xi) w iv x iv ≥ C (2) v∈D(xi) x iv = 1, ∀x i ∈ S (3)
x iv ∈ {0, 1},

∀x i ∈ S, v ∈ D(x i ) (4)
We call this ILP ∅ . The main property of ILP ∅ is that the cost of any feasible complete assignment is equal to the cost of the corresponding tuple in c S after any sequence of cost moves. Hence, opt(ILP ∅ ) > 0, if and only if c S is not ∅IC, and we can move some cost to c ∅ : project(c ∅ , c S , ∅, opt(ILP ∅ )).

However, for the purposes of detecting violations of F ∅IC, it is not enough to look at the cost of tuples of the constraint, as we must also take unary costs into account. Therefore, while ILP ∅ remains the representation of the constraint, the propagator considers the problem with the modified objective min xi∈S,v∈D(xi)

(δ iv + c i (v))x iv -δ ∅ (5)
Let this problem be ILP F ∅ . c S is F ∅IC if and only if opt(ILP F ∅ ) = 0. In the following, we write p iv = δ iv + c i (v) for compactness, when it does not matter how much of the coefficient came from δ iv and how much came from c i (v). In contrast with ILP ∅ , if opt(ILP F ∅ ) > opt(ILP ∅ ), we cannot move opt(ILP F ∅ ) units of cost to c ∅ . Instead, we first have to move some cost from unary cost functions into the constraint before we can project it to c ∅ . In this case, the composition of p iv from δ iv and c i (v) is significant.

Unfortunately, ILP ∅ and ILP F ∅ have the knapsack problem as a special case, hence it is NP-hard to determine whether a PB constraint is ∅IC or F ∅IC. Therefore, we detect only a subset of cases where the constraint is not F ∅IC by relaxing the integrality constraint (4) into 0 ≤ x iv ≤ 1 and solving the resulting linear programs, called LP ∅ and LP F ∅ , respectively. This forgoes the guarantee that opt(LP F ∅ ) = 0 if and only if the constraint is F ∅IC, and satisfies only the 'only if' part. More simply, if opt(LP F ∅ ) > 0 then the constraint is not F ∅IC, and similarly for LP ∅ and ∅IC.

LP F ∅ has a special structure. It is a Multiple-Choice Knapsack Problem (MCKP) [START_REF] Pisinger | Knapsack problems[END_REF], or a knapsack problem with special ordered sets [START_REF] Johnson | A note of the knapsack problem with special ordered sets[END_REF]. These can be solved more efficiently than arbitrary LPs, a fact that we use in our propagator.

Solving the Knapsack LP

We obtain an optimal solution x * of the primal LP F ∅ by applying Pisinger's greedy algorithm [START_REF] Pisinger | Knapsack problems[END_REF]. This gives a x * in time O(N log N ) 5 , with N = |x * |, such that either x * has no fractional value or it has exactly two fractional values. In the latter case, the WCSP variable x k ∈ S, verifying ∃s, s ∈ D(x k ) such that 0 < x * ks , x * ks < 1, is called a split class and x ks , x ks are the split variables. We denote by o = xi∈S,v∈D(xi) p iv x * iv -δ ∅ , the optimal solution cost of LP F ∅ . Consider now the dual of LP F ∅ :

max C × y cc + xi∈S y i (6)
s.t.

y cc × w iv + y i ≤ p iv ∀x i ∈ S, v ∈ D(x i ) y cc ≥ 0
Where y cc is the dual variable corresponding to the capacity constraint and y i corresponds to the EO constraint of x i . From the optimal primal solution, it is easy to compute the optimal dual solution. Let x k be the split class, x ks , x ks the split variables and for i = k, define the variable x is as the variable used in the optimal solution, i.e., x * is = 1.

y cc = p ks -p ks w ks -w ks y k = p ks -y cc × w ks = p ks -y cc × w ks y i = p is -y cc × w is ∀x i ∈ S \ {x k }
From the dual solution y, we compute the reduced cost rc y (x iv ) of every variable x iv , i.e., the slack of the dual constraint that corresponds to x. When context makes it clear, we omit y and write rc(x iv ).

The reduced cost of a variable x can be interpreted as the amount by which we must decrease the coefficient of x in the objective function in order to have x > 0 in the optimal solution. We explain later that this implies that we can project some cost to unary cost functions.

In the specific case of LP F ∅ , we have:

rc(x ks ) = rc(x ks ) = 0 rc(x is ) = 0 ∀x i ∈ S \ {x k } rc(x iv ) = p iv -y cc × w iv -y i ∀x i ∈ S, v = s Observation 1 Consider the linear program LP F ∅ which is identical to LP F ∅ but has p iv = p iv -rc(x iv ). Then opt(LP F ∅ ) = opt(LP F ∅ ).
Proof. The optimal solution x * of LP F ∅ has the same cost o in LP F ∅ and LP F ∅ , as the coefficients of the variables that are greater than 0 are unchanged. The optimal dual solution x * remains feasible in LP F ∅ , as the slack in the dual of LP F ∅ matches exactly the reduction in the right-hand side. Moreover, as the dual objective did not change, it has the same cost and matches the primal cost, so opt(LP

F ∅ ) = o = opt(LP F ∅ ).
Example

x iv = 1 ∀x i ∈ {x 1 , x 2 } 0 ≤ x iv ≤ 1 ∀x i ∈ {x 1 , x 2 }, v ∈ D(x i )
Pisinger's algorithm gives the optimal primal solution x * = {0, 1, 0, 7 12 , 5 12 } with cost o = 55 + 7 12 × 47 + 5 12 × 95 = 122. We deduce the following dual optimal solution : y cc = 2, y 1 = 55 -2 × 14 = 27, y 2 = 47 -2 × 16 = 15. The following reduced costs are obtained : rc(x 12 ) = rc(x 21 ) = rc(x 22 ) = 0 and rc(x 11 ) = 5, rc(x 13 ) = 10, we deduce that replacing the previous objective function by the following one does not change the cost of the optimal solution: min 35x 11 + 55x 12 + 75x 13 + 47x 21 + 95x [START_REF] De Givry | Clique cuts in weighted constraint satisfaction[END_REF] We observe that the solution x * = {0, 1, 0, 7 12 , 5 12 } is still optimal.

Propagation

Given a PB constraint and the associated unary costs, it is possible to increase the lower bound by at least opt(LP F ∅ ). Our goal is to extend as little cost as possible from the unary cost functions in order to make opt(LP ∅ ) = o = opt(LP F ∅ ) and then project o to c ∅ .

Procedure TransformPB(c S , y cc , y i , o)

Data: c S : PB constraint Data: ycc, yi, o: optimal dual solution of LP F ∅ 1 for all the variables xiv do

2 ci(v) ← ci(v) -ycc × wiv -yi + δiv 3 δiv ← ycc × wiv + yi 4 c∅ ← c∅ + o 5 δ∅ ← δ∅ + o
If we move |c i (v) -rc(x iv )| between the constraint and each unary cost function and value, then opt(LP ∅ ) = o and we can project o to c ∅ . Indeed we have |c i (v) -rc(x iv )| = |(y cc × w iv + y i ) -δ iv |, we thus obtain the EPTs performed by Procedure TransformPB.

Theorem 1. Algorithm TransformPB preserves equivalence.

Proof. Recall that iv = c i (v) + δ iv and that rc(x iv ) ≥ 0. If c i (v) -rc(x iv ) ≥ 0 then the cost move is an extension of less than c i (v), it is valid. If c i (v)-rc(x iv ) < 0 then the cost move is a projection, while the cost of any solution x with x iv = 1 is at least o -c i (v) + rc(x iv ). This operation is also valid.

Finally, to check that our sequence of EPTs justifies the increase of c ∅ by o, we compute the optimum of LP ∅ . From Observation 1, opt(LP ∅ ) = o, which means we can project o to c ∅ and increase δ ∅ to bring opt(LP ∅ ) = opt(LP F ∅ ) = 0.

We can improve on this by observing that the integer optimum must be integral. Therefore, we can increase c ∅ by o . In this case, it is also necessary to round up all cost moves. By rounding up, we can no longer rely on Observation 1, but it still holds that opt(LP ∅ ) = 0. We also approach ∅IC by verifying that for any value x ab we have δ ab + min xi∈S\xa,v∈D(xi) (δ iv + c i (v))x iv -δ ∅ = 0. If this is not the case, we can project a positive cost to c a (b).

Procedure Propagate is the entry point to the propagator. It enforces domain consistency on the PB constraint, then solves LP F ∅ . If there is more than one optimal solution we prefer the one minimizing the reduced cost of the EAC support of each variable. Finally, it uses Procedure TransformPB, to perform cost moves. Example 3. Returning to Example 2, where c S is the PB constraint with EO partitions over two WCSP variables x 1 and x 2 , we had the following reduced costs: rc(x 12 ) = rc(x 21 ) = rc(x 22 ) = 0, rc(x 11 ) = 5, rc(x 13 ) = 10, the optimal cost was 122. We deduce the following cost moves:

-extend(c 1 , {x 1 = 1}, c S , 35) -extend(c 1 , {x 1 = 2}, c S , 55) -extend(c 1 , {x 1 = 3}, c S , 75) -extend(c 2 , {x 2 = 1}, c S ,
x iv = 1 ∀x i ∈ {x 1 , x 2 } 0 ≤ x ij ≤ 1 ∀x i ∈ {x 1 , x 2 }, v ∈ D(x i )
The optimal solution is x * = {0, 0, 1, 1, 0} and its cost is o = 76 + 47 -122 = 1. We deduce the dual optimal solution y cc = 1, y 1 = 52, y 2 = 31 with reduced costs rc(x 11 ) = rc(x 13 ) = rc(x 21 ) = 0 and rc(x 12 ) = 19, rc(x 22 ) = 24. We carry out the following cost moves: extend(c 1 , {x 1 = 1}, c S , 21), extend(c 1 , {x 1 = 2}, c S , 11), extend(c 1 , {x 1 = 3}, c S , 1), project(c 2 , c S , {x 2 = 2}, 24), project(c ∅ , c S , ∅, 1), with δ 11 = 56, δ 12 = 66, δ 13 = 76, δ 21 = 47, δ 22 = 71, δ ∅ = 123.

Experimental results

We implemented our approach in toulbar2, an exact WCSP solver in C++, 6 winner of past UAI-2008, 2014 competitions. toulbar2 default variable ordering heuristic is the weighted degree heuristic [START_REF] Boussemart | Boosting systematic search by weighting constraints[END_REF], in order to gain information from PB constraints, we adapted an explanation-based weighted degree for linear inequality presented by Hebrard and Siala [START_REF] Hebrard | Explanation-based weighted degree[END_REF]. For all the tests we imposed a time limit of 30 minutes (except for CPD with 1 hour) on a single core of an Intel Xeon E5-2680 v3 at 2.50 GHz and 256 GB of RAM. We compared our PB propagator with other modeling approaches in protein design. We compared toulbar2 to state-of-the-art ILP solver cplex 20.1 on knapsack problems with conflict graphs. We also compared toulbar2 on pseudo-Boolean Competition 2016 (previously out of reach by toulbar2) but the results were not competitive with recent PB solvers (not reported here for the lack of space).

Sequence of diverse solutions for CPD

A protein is a chain of simple molecules called amino acids. This sequence determines how the protein will fold into a specific 3D shape. The Computational Protein Design (CPD) [START_REF] Allouche | Computational protein design as an optimization problem[END_REF] problem consists of identifying the sequence of amino acids that should fold into a given 3D shape. This problem can be modeled as a CFN 7 with unary and binary cost functions representing the energy of the protein but the criteria only approximate the reality, thus producing a sequence of diverse solutions increases the chance of finding the correct real sequence of amino acids. Each time a solution is found, a Hamming distance constraint is added to the model to enforce the next solution to be different from the previous ones. This Hamming distance can be directly encoded as a PB linear constraint, in the form of Eq. ( 2), with EO partitions associated to domains (Eq. ( 3)). For each variable, a negative weight of -1 is associated to the value found in the last solution (other values having a zero-weight) and the weighted sum in Eq. ( 2) must be greater than or equal to -(|X|-ζ), where ζ corresponds to the required minimum Hamming distance.

This has been implemented in toulbar2 and compared to previous automatabased encoding approaches (ternary, hidden, and dual encodings from [START_REF] Ruffini | Guaranteed diversity and optimality in cost function network based computational protein design methods[END_REF]) on 30 instances [START_REF] Traoré | A new framework for computational protein design through cost function network optimization[END_REF]. 8 Selected instances have from 23 to 97 residues/variables with 6 https://github.com/toulbar2/toulbar2 version 1.2. 7 Other paradigms such as ILP or Max-SAT have been tested but the experimental results using their corresponding state-of-the-art solvers were inferior to the CFN approach using toulbar2 -divm=(0 for dual, 1 for hidden, 2 for ternary, and 3 for the PB encoding)). Figure 1 reports the solving time of each encoding. The dual and ternary encodings failed to give 10 diverse solutions for one instance, while the PB and hidden encodings didn't. Moreover the PB encoding is faster for 29 instances and it solves 23 of them in less than 30 seconds while dual, hidden, and ternary encodings solve respectively 12, 13, and 4 instances in less than 30 seconds. Note that we are computing a greedy sequence of solutions, the different encodings do not return the exact same sequence (except for 7/30 instances). We also compared for each instance the number of backtracks and time (not reported here) of the previous toulbar2 default encoding (dual) with the PB encoding. In all the instances the PB encoding needs fewer backtracks than the dual encoding and except for one instance, the PB encoding is also faster. Automata-based encodings have the flaw of introducing extra variables that can disturb the variable ordering heuristic (by default, min domain size over weighted degree [START_REF] Boussemart | Boosting systematic search by weighting constraints[END_REF]) 9and local consistency algorithm (by default, EDAC during search, except partial F∅IC for PB constraints). While the PB encoding directly encodes the Hamming distance, it is heavier to propagate as we can see by comparing the number of backtracks per second (170 for PB encoding and 1060 for dual encoding). 1. Number of solved instances (left) and number of times a solver found the best solution within the time limit (right) for six different classes of KPCG.

Knapsack problem with a conflict graph

We compare here toulbar2 and cplex on Knapsack with Conflict Graph (KPCG) [START_REF] Bettinelli | A branch-and-bound algorithm for the knapsack problem with conflict graph[END_REF][START_REF] Coniglio | A new combinatorial branch-and-bound algorithm for the knapsack problem with conflicts[END_REF], a knapsack problem combined with binary constraints representing conflicts between pairs of variables. We use 6 different classes C1,C3,C10,R1, R3,R10. In three of them the weight and the profit of each variable are correlated (class C) otherwise the profit is random between [START_REF] Allouche | Operations Research and Simulation in Healthcare, chap. Cost Function Networks to Solve Large Computational Protein Design Problems[END_REF]100] (class R). The numbers 1, 3, 10 correspond to a multiplying coefficient of the capacity, which has the effect of making the instances harder as the multiplier increases. In each class half of the instances have capacity 150, weights are uniformly distributed in [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF]100], and the number of Boolean variables varies between 120, 250, 500, and 1000. For the other half, the capacity is 1000, weights are uniformly distributed in [250, 500], and the number of Boolean variables varies between 60, 120, 349, and 501. Additionally, the density of the conflict graph varies from 0.1 to 0.9. In total, each class has 720 instances. We used a direct encoding for toulbar2. For cplex, we tried with both tuple and direct encodings (tuple encoding corresponds to the local polytope with integer variables) [START_REF] Hurley | Multi-language evaluation of exact solvers in graphical model discrete optimization[END_REF]. Table 1 reports the number of instances solved by each solver. toulbar2 was more efficient than cplex with the tuple encoding and competitive with cplex using the direct encoding for four out of six classes. Moreover, toulbar2 finds the best solutions for the largest number of instances in every class.

Conclusion and future work

It is now possible to model pseudo-Boolean linear constraints in deterministic and probabilistic graphical models. This provides greater modeling flexibility and allows a WCSP solver like toulbar2 to solve more problems, such as computational protein design problems with diversity guarantee or knapsack problems with conflict graphs. One of the weaknesses of our approach is that the algorithm fundamentally produces a suboptimal solution to the linear program, because it propagates the pseudo-Boolean linear constraints one by one and does not take into account other constraints (except at-most-one constraints). There are several ways to improve this, including adapting work previously done in this context on Lagrangian relaxation [START_REF] Komodakis | MRF energy minimization and beyond via dual decomposition[END_REF] or an approach closer to VAC [START_REF] Dlask | Bounding linear programs by constraint propagation: Application to Max-SAT[END_REF]. It also opens up possibilities for other uses of linear constraints in the WCSP framework, such as the generation of cuts.

Example 1 .

 1 Consider two variables x, y with D(x) = D(y) = {a, b, c} and three cost functions c x , c y , c {x,y} such that the only non-zero costs are c x (a) = c y (a) = 1, c {x,y} ({x = b, y = b}) = c {x,y} ({x = b, y = c}) = c {x,y} ({x = c, y = b}) = 1, and c {x,y} ({x = c, y = c}) = 2.

Theorem 2 .

 2 Procedure Propagate runs in O(nd log nd) time where n is the number of WCSP variables involved and d the maximum domain size. Proof. Pisinger's algorithm dominates the complexity, as it runs in O(N log N ), where N is the number of LP variables. In our case, N = nd, so it takes O(nd log nd) time. Domain consistency on the linear inequality can be performed in linear time. Finally, Procedure TransformPB iterates once over all variables and values and performs constant time operations on each. Hence, the total complexity is O(nd log nd). Procedure Propagate(c S ) Data: c S : PB constraint with EO partitions 1 DomainConsistency(c S ) 2 (ycc, yi, o) = DualSolve(LP F ∅ ) 3 TransformPB(c S , ycc, yi, o)

  47) extend(c 2 , {x 2 = 2}, c S , 95) project(c ∅ , c S , ∅, 122) It implies the resulting costs: δ 11 = 35, δ 12 = 55, δ 13 = 75, δ 21 = 47, δ 22 = 95, δ ∅ = 122. The unary costs after these operations are c 1 (2) = c 2 (1) = c 2 (2) = 0, c 1 (1) = 5, c 1 (3) = 10. If we construct the table of possible assignments of LP ∅ obtained after the extensions, we can see that c S ({x 1 = 1, x 2 = 2}) = 8, c S ({x 1 = 2, x 2 = 2}) = 28, c S ({x 1 = 3, x 2 = 1}) = 0, c S ({x 1 = 3, x 2 = 2}) = 48, and all the other assignments don't satisfy the constraint. We observe that the optimal solution is 0, hence our extensions justify the increase of c ∅ . Now assume that other EPTs outside the PB constraint have modified the unary costs: c 1 (1) → c 1 (1) + 16 = 21, c 1 (2) → c 1 (2) + 30 = 30, c 1 (3) → c 1 (3) -9 = 1. We want to compute a new lower bound for the PB constraint by solving LP F ∅ : min 56x 11 + 85x

Fig. 1 .

 1 Fig. 1. Cactus plot of CPU solving time (log scale) for different encodings of Hamming distance constraints on CPD.

  12 + 76x 13 + 47x 21 + 95x 22 -122 s.t. 4x 11 + 14x 12 + 24x 13 + 16x 21 + 40x 22 ≥ 40

	v∈D(xi)

An extension to ternary cost functions has been proposed[START_REF] Sánchez | Mendelian error detection in complex pedigrees using weighted constraint satisfaction techniques[END_REF] but it requires managing all scope intersections and not only unary cost functions.

E.g., EDAC[START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF], an extension of EAC property, is maintained in O(ed 2 max(nd, )) for a WCSP with n variables, maximum domain size d, and e binary cost functions.

The Dyer-Zemel algorithm[START_REF] Dyer | An o(n) algorithm for the multiple-choice knapsack linear program[END_REF][START_REF] Zemel | An o(n) algorithm for the linear multiple choice knapsack problem and related problems[END_REF] can compute a solution in O(N) time, but we have not yet implemented it.

Additionally, the PB constraint provides finer-grain weights using explanations[START_REF] Hebrard | Explanation-based weighted degree[END_REF] when linear coefficients are not all equal as it is the case in the KPCG benchmark.

This research was funded by the French "Agence Nationale de la Recherche" through grants ANR-18-EURE-0021 and ANR-19-P3IA-0004.


Elisabet Rothenberg

Stephan C Bischoff

Ronan Thibault
email: ronan.thibault@chu-rennes.fr

Prof Ronan Thibault

Keywords: Texture modification, International Dysphagia Diet Standardisation Initiative (IDDSI), Terminology

mentioning several different systems for describing the consistency of food and drink [START_REF] Thibault | ESPEN guideline on hospital nutrition[END_REF] was to highlight that there is currently no uniform terminology used in this area. Different systems for texture modification exsist, terminology regarding different textures differ, and in the routine practices, prescriptions vary between countries and health care systems [START_REF] Rothenberg | Reply to -Letter to the Editor by Pedrolli C entitled 'IDDSI: Worth or not?[END_REF]. Thus, our intention was merely to mention IDDSI and refer to IDDSI website without further description. Although there are several different systems for describing textures of food and drinks, we acknowledge that IDDSI has achieved widespread use internationally in recent years and, for that reason, we would like to mention this system of texture classification [START_REF] Rothenberg | Reply to -Letter to the Editor by Pedrolli C entitled 'IDDSI: Worth or not?[END_REF].

To make it clearer for the readers, we propose to add to the paragraph 3.10 [1] the underlined words or sentences, so that it should be read as follows:

"Texture modification of food and/or drink is an important intervention used so that people with dysphagia can swallow effectively and safely. However, the different names for and number of levels of modification and the characteristics used within and across countries might increase the risk to patient safety. One internationally recognized standardized system for evaluating and describing different levels of texture modified food and thickened fluids is the International Dysphagia Diet Standardisation Initiative (IDDSI), which provides a common terminology for food textures and drink thickness (https://iddsi.org). In countries where IDDSI is not implemented, a variety of non-standardised descriptors for different textures can be used. We propose to describe them as follows:

J o u r n a l P r e -p r o o f 4 -Liquidized/thin puree; homogenous consistency that does not hold its shape after serving.

-Thick puree/soft and smooth; thickened, homogenous consistency that holds its shape after serving and does not separate into a liquid and solid component during swallowing, i.e., cohesive.

-Timbal: homogenous smooth consistency that is omelette-like in texture and made from smooth purees mixed with egg and then baked. Timbal holds its shape after serving, is not sticky and does not separate into a liquid and solid component after serving or during swallowing, i.e., cohesive. Can be eaten with a spoon or fork.

-Finely minced; soft diet of cohesive, consistent textures requiring some chewing (particle size most often described as 0.5 * 0.5 cm).

-Modified normal; normal foods of varied textures that require chewing, avoiding particulate foods that pose a choking hazard (particle size most often described as 1.5 * 1.5 cm)."

In conclusion, we hope that you and many others will profit from this guideline and thank again for the comments.

Acknowledgement

None.

J o u r n a l P r e -p r o o f

We thank very much Riquelme et al. for their interest in our recently published 'ESPEN guideline on hospital nutrition' [1]. We have noticed their concern regarding the different textures mentioned in the paragraph 3.10 entitled 'Texture modified food and thickened fluids' [1].

We would like making clear that the different textures described are not specifically referring to the International Dysphagia Diet Standardisation Initiative (IDDSI). The reason for

Funding statement

This ESPEN hospital nutrition guideline was solely financed by ESPEN, the European Society for Clinical Nutrition and Metabolism.

Author contributions

ER, SCB & RT wrote the letter; all three authors approved the final submitted version of the manuscript.

Conflict of interest

The authors declare no conflict of interest with regard to this letter.