Clement Viricel

David Simoncini

Sophie Barbe

Thomas Schiex

Clément Viricel

Guaranteed Weighted Counting for Affinity Computation: Beyond Determinism and Structure

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Graphical models [START_REF] Dechter | Reasoning with probabilistic and deterministic graphical models: Exact algorithms[END_REF] are sparse representations of highly dimensional multivariate distributions that rely on a factorization of the distribution in small factors. When variables are discrete, graphical models cover a variety of mathematical models that represent joint discrete distributions (or functions) that can be either Boolean functions (e.g., in propositional satisfiability SAT and constraint satisfaction CSP), cost functions (as in partial weighted MaxSAT and Cost Function Networks [START_REF] Cooper | Soft arc consistency revisited[END_REF]) or probability distributions (in stochastic models such as Markov Random Fields and Bayesian networks).

Typical queries on such graphical models are either optimization or counting queries (or a mixture of these). In optimization queries, we look for an assignment that maximizes the joint function, i.e., a model in SAT, a solution in CSP or a Maximum a posteriori assignment (MAP) in a Markov Random Field (MRF). All these problems have an associated NP-complete decision problem.

Counting problems are central in stochastic graphical models because they capture the computation of marginal probabilities on subsets of variables and the computation of the normalizing constant Z that is required to define a probability distribution from the non-normalized distribution of Markov Random Fields. This difficult problem requires a summation over an exponential number of elementary terms and is #P-complete [START_REF] Valiant | The complexity of computing the permanent[END_REF]. As shown by [START_REF] Toda | On the computational power of pp and ⊕p[END_REF], one call to a #-P oracle suffices to solve any problem in the Meyer-Stockmeyer polynomial hierarchy in deterministic polynomial time, an indication that it could be outside of the PH.

Computing Z is a central problem in statistics (e.g., for parameter estimation in MRFs), for Bayesian network processing (to account for evidence) and is also crucial in statistical physics where it is called the partition function. A typical domain where partition function computation can be extremely useful is computational protein design. Indeed, the affinity of a protein for a specific target molecule can be estimated by modeling both molecules as MRFs representing physics force fields and by computing the two partition functions: one for the bound protein and target and another for the same molecules in unbound state [START_REF] Silver | Efficient computation of small-molecule configurational binding entropy and free energy changes by ensemble enumeration[END_REF].

For these reasons, various approaches have been designed to tackle this problem. The Mean-Field algorithm [START_REF] Jaakkola | Advanced mean field methods: theory and practice[END_REF], Tree-reweighted Belief Propagation [START_REF] Wainwright | A new class of upper bounds on the log partition function[END_REF] as well as more recent proposals [START_REF] Liu | Bounding the partition function using holder's inequality[END_REF] have been proposed, but they do not offer any formal guarantee on the quality of the approximation they produce, except in very special cases. Monte-Carlo methods including Markov Chain Monte Carlo methods [START_REF] Gilks | Markov chain monte carlo[END_REF] offer asymptotic convergence /s, but convergence is impractically slow. Indeed, there are recent significant examples showing that the time needed for Monte Carlo methods to converge can be easily under-estimated [START_REF] Simoncini | Guaranteed discrete energy optimization on large protein design problems[END_REF]. Practical MCMC based tools also rely on heuristics that destroy these theoretical guarantees. More recent stochastic methods exploiting universal hashing functions offer "Probably Approximately Correct" (PAC) estimators [START_REF] Ermon | Taming the curse of dimensionality: Discrete integration by hashing and optimization[END_REF][START_REF] Chakraborty | Distribution-aware sampling and weighted model counting for sat[END_REF]. Here, a bound δ on the probability that the estimation does not lie within a (1 + ε) ratio of the true value is set and a corresponding estimation produced.

Finally, different methods, mostly based on SAT-solvers, have been defined that can perform exact weighted model counting (#SAT) with deterministic guarantees, a problem to which the problem of computing Z can be easily reduced. To avoid the exponential blowup in the number of terms to add, solvers providing deterministic guarantees rely on two independent ideas: exploiting determinism (zero weights) to prune regions of the space that do not contribute to the sum, and exploiting independence which may be detected at the graphical model structure level, as captured by its tree-width, but also at a finer level as context-sensitive independence [START_REF] Sang | Solving bayesian networks by weighted model counting[END_REF] 3 . Independence enables caching of intermediate counts that can be factored out and lead to exponential time savings at the cost of memory. The very same ideas are also exploited in knowledge compilers that may compile graphical models or SAT formulas to languages on which counting becomes easy [START_REF] Chavira | On probabilistic inference by weighted model counting[END_REF][START_REF] Oztok | A top-down compiler for sentential decision diagrams[END_REF].

In this paper, we explore the possibility of preserving the deterministic guarantees of exact solvers and explore a new source of pruning that may be present even when determinism or independence are too limited to allow for exact counting: detecting and pruning regions for which it is possible to prove, at limited cost, that they contain an amount of weight which is too small to significantly change the computed value of Z. Instead of providing a PAC guarantee, our algorithm provides an approximation of the normalizing constant that is guaranteed to lie within a ratio of 1 + ε of the true value (with probability 1), a guarantee that none of the PAC randomized algorithms above can provide in finite time.

Our initial motivation for computing Z lies in Computational Protein Design (CPD). The aim of CPD is to design new proteins that have desirable properties which are not available in the existing catalog of known proteins. One of these properties is the affinity between a protein and another molecule (such as another protein, a peptide, an amino-acid, a small organic molecule, etc. . . ). The binding affinity gives an indication of the likelihood that two molecules will prefer to bind together rather than remain dissociated and thus that a protein will be likely to bind to another molecule of interest. Proteins can be described as a set of bound atoms subjected to a number of atom scale forces captured by a pairwise force field defining a Markov Random Field [START_REF] Pearlman | Amber, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules[END_REF]. From this MRF, the binding affinity can be estimated by computing the ratio of the partition functions of the molecules in bound and unbound states [START_REF] Silver | Efficient computation of small-molecule configurational binding entropy and free energy changes by ensemble enumeration[END_REF][START_REF] Georgiev | The minimized dead-end elimination criterion and its application to protein redesign in a hybrid scoring and search algorithm for computing partition functions over molecular ensembles[END_REF].

In the rest of the paper, after introducing our notations and the binding affinity computation problem, we present the Z * ε algorithm, a variant of Branch and Bound targeted at counting instead of optimizing. Z * ε relies on the availability of a local upper bound on Z. We then consider different simple, fast, safe and incremental upper bounds on Z, integrate them in Z * ε and compare them to exact counting tools on two categories of benchmarks: general benchmarks extracted from the UAI and Probabilistic Inference (PIC'2011) challenges and partition function computation problems appearing as subproblems of binding affinity computation on real proteins. Surprisingly, despite a very limited caching strategy, the resulting algorithm is able to outperform exact solvers on a variety of problems and is especially efficient on CPD-derived problems. Because Z * ε relies on a new source of pruning, its underlying principle and associated bounds can be immediately used to improve existing SAT-based counters using Max-SAT bounds, which are closely related to local consistencies in Cost Function Networks [START_REF] Larrosa | Resolution in Max-SAT and its relation to local consistency in weighted CSPs[END_REF][START_REF] Bonet | Resolution for max-sat[END_REF][START_REF] Larrosa | A logical approach to efficient max-sat solving[END_REF].

Background

A Markov Random Field defines a joint probability distribution over a set of variables as a factorized product of local functions, usually denoted as potential functions.

Definition 1 A discrete Markov Random Field (MRF) is a pair (X, Φ) where X = {1, . . . , n} is a set of n random variables, and Φ is a set of potential functions. Each variable i ∈ X has a finite domain D i of values that can be assigned to it. A potential function φ S ∈ Φ, with scope S ⊆ X, is a function φ S : D S → R ∪ {∞} where D S denotes the Cartesian product of all D i for i ∈ S.

The energy or potential of an assignment t ∈ D X is denoted as

E(t) = φ S ∈Φ φ S (t[S])
where t[S] is the projection (or restriction) of t to the variables in S. Notice that this definition shows that an MRF is essentially equivalent to a Cost Function Network (or WCSP [START_REF] Cooper | Soft arc consistency revisited[END_REF]). A tuple t ∈ D S will be represented as a set of pairs {(i,

t[i]) | i ∈ S}.
The probability of a tuple t ∈ D X is then defined as:

P (t) = exp(-E(t)) t ∈D X exp(-E(t ))
The normalizing constant below the fraction is usually denoted as Z. The potential φ S are called energies, in relation with statistical physics. An assignment with minimum energy has therefore maximum probability. With pairwise potentials (|S| ≤ 2), an MRF defines a graph with variables as vertices and potential scopes S as edges. In the rest of this paper, for the mere sake of simplicity and w.l.o.g., we assume pairwise MRFs including also unary potential functions and a constant φ ∅ potential function. We denote by d the maximum domain size and e the number of pairwise potential functions. Using table representations, a pairwise MRF requires O(ed 2 ) space to be represented. Note that Bayesian networks can be seen as specific MRFs enforcing a local normalization condition of potentials and a specific DAG-base graph structure, that together guarantee that Z = 1. As soon as evidence (observations) change the domain of the variables however, Bayesian networks become unnormalized and computing Z becomes #P-complete in general.

Computational Protein Design and Binding Affinity

Proteins are linear chains of small molecules called "amino-acids". There are 20 natural different amino-acids. All amino-acids share a common core and the cores of all successive amino-acids in a proteins are linked together to form a linear chain, called the protein backbone. Each amino-acid also has a variable side-chain which chemical nature defined the precise amino-acid used. This lateral chain is highly flexible. The structure of a protein in 3D-space is therefore characterized by the shape of the linear chain itself (the backbone), and the specific spatial orientation of all side-chains, at each position of the chain. Proteins are universally present in the cells of all living organisms and perform a vast array of functions including catalyze, signaling, recognition, transporting, repair. . . Proteins differ from one another primarily in their sequence of amino-acids which usually results in protein folding into a specific 3D structure that determines its function. The characteristic of proteins that also allows their diverse set of functions is their ability to bind other molecules, with high affinity and specificity. See [START_REF] Campeotto | A constraint solver for flexible protein model[END_REF][START_REF] Allouche | Computational protein design as an optimization problem[END_REF] for an intrduction to proteins targeted at the CP audience.

Proteins have a relatively stable general shape. The relative stability of a molecule in a given conformation can be evaluated by computing its energy, lower energy states being more stable. This energy is derived from various molecular forces including bond angles, electrostatic forces, molecular clashes and distances. It can be computed using existing force fields such as Amber [START_REF] Pearlman | Amber, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules[END_REF], the one used in our experiments. Notice that molecular clashes -interpenetrating atoms -may generate infinite energies i.e., determinism.

Despite a plethora of functionalities of proteins, there is still an ever-increasing demand for proteins endowed with specific properties of interest for many applications (in biotechnology, synthetic biology, green chemistry and nanotechnology) which either do not exist in nature or have yet not been found in the biodiversity. To this end, Computational structure-based Protein Design (CPD) has become a key technology. By combining physico-chemical models governing relations between protein amino-acid composition and protein 3D structure with advanced computational algorithms, CPD seeks to identify one or a set of amino-acid sequences that fold into a given 3D structure and possess the targeted properties. This in silico search for the best sequence candidates opens up new possibilities to better guide protein engineering by focusing experimentation on the relevant sequence space for the desired protein function and thereby reducing the size of mutant libraries that need to be built and screened. In recent years, CPD has experienced important success, especially in the design of therapeutic proteins [START_REF] Miklos | Structure-based design of supercharged, highly thermoresistant antibodies[END_REF], novel enzymes [START_REF] Röthlisberger | Kemp elimination catalysts by computational enzyme design[END_REF], protein-protein interfaces [START_REF] Karanicolas | Computational design of affinity and specificity at proteinprotein interfaces[END_REF][START_REF] Sammond | Computational design of secondsite suppressor mutations at protein-protein interfaces[END_REF], and large oligomeric ensembles [START_REF] King | Accurate design of co-assembling multi-component protein nanomaterials[END_REF]. Nevertheless, the computational design of proteins with defined affinity for a given molecule (such as a small organic, a peptide, another protein. . . ) which is essential for large range of applications, continues to present challenges. A traditional approach to model proteins in CPD is to assume that their backbone is totally rigid and that only side-chains move, each side-chain being able to adopt a discrete set of most likely conformations defined in a so-called "rotamer" library (see Figure 1). We use the Penultimate rotamer library [START_REF] Lovell | The penultimate rotamer library[END_REF].

With one variable per side-chain, each with a domain equal to the set of available rotamers for this side-chain and a pairwise decomposable energy function such as Amber force field, a protein naturally defines a pairwise MRF with a rather dense graph. The partition function Z of this MRF captures important properties of the protein. Specifically, the association constant (or binding constant) is used to describe the affinity between a protein and a ligand (a protein or another molecule of interest). This association constant can be estimated by computing the partition function of the two molecules in bound and unbound states. The ratio of these two partition functions being proportional to their affinity.

From a computational point of view, an important property of proteins of interest is that their general shape is stable which means that the proportion of low energy (or high probability) states among the exponential number of possible states is likely to be very small. On the opposite side of the energy scale, the infinite energies created by molecular clashes means that there will be states with 0 probability. This is favorable for exact solvers that can exploit determinism to speedup Z computation. It however means that CPD instances will exhibit unbounded tilt (defined in [START_REF] Chakraborty | Distribution-aware sampling and weighted model counting for sat[END_REF] as the ratio τ = max t∈D X P (t) min t∈D X P (t) ). This situation is not ideal for the WeightMC PAC algorithm which requires a finite upper-bound on τ to run in finite time.

Guaranteed counting

Because it is rarely (if ever) needed to compute a probability or a partition function with an absolute precision (which is also inherently limited by finite representations), we consider the general problem of computing an εapproximation Ẑ of Z, i.e., such that:

Z 1 + ε ≤ Ẑ ≤ Z (1)
Such approximation allows us to compute an estimate P (t) = exp(-E(t)) Ẑ such that

P (t) ≤ P (t) ≤ (1 + ε)P (t).
In the context of #-SAT, it has been shown that providing such relative approximations remains intractable for most of the known SAT polynomial classes [START_REF] Roth | On the hardness of approximate reasoning[END_REF]. As we will see, it can however be exploited to prune during polynomial space depth-first tree-search based counting and sometimes provide important speedups.

Assuming that for any MRF, and any assignment t of some of its variables, we can compute an upper bound U b(t) of the partition function of the MRF where variables are assigned as in t, the Depth First Branch and Bound schema used for exactly solving optimization problems on cost function networks [START_REF] Cooper | Soft arc consistency revisited[END_REF][START_REF] Allouche | Computational protein design as an optimization problem[END_REF] can be adapted to compute Z [START_REF] Viricel | Approximate counting with deterministic guarantees for affinity computation[END_REF].

Function Z * ε (t, V ) if V = ∅ then 1 Ẑ ← Ẑ + exp(-E(t)); 2 else Choose i ∈ V ; 3 for a ∈ D i do t ← t ∪ {(i, a)}; if (U + U b(t ) + Ẑ ≤ (1 + ε) Ẑ) then 4 U ← U + U b(t ); 5 else Z * ε (t , V -{i});
Algorithm 1: Guaranteed approximate counting. Initial call: Z * ε (∅, X). U and Ẑ are global variables initialized to 0.

The algorithm simply explores the tree of all possible assignments of the MRF, starting with the whole set of unassigned variables (in V ), choosing an unassigned variable (line 3), trying all possible values. When all variables are assigned (line 1), the contribution of the complete assignment t is accumulated in a running count which will eventually define the approximation Ẑ (line 2). However, branches which provide a sufficiently small mass of probability (as estimated by U b(t )) are pruned and this overestimation of the neglected mass is accumulated in U (line 5). Because pruning may occur, eventually, Ẑ will be a lower bound of Z.

Theorem 1. Z * ε terminates and returns an ε-approximation of Z.

Proof. The termination follows from the fact that Z * ε explores a finite tree. We now show that the algorithm always provides a ε-approximation. When the algorithm finishes, all the assignments have either been explored (line 2) and counted or pruned (line 5). Since U is the sum of all the upper bounds on the mass of probability in all pruned branches, we have that Ẑ + U ≥ Z. Initially, Ẑ = U = 0 and the invariant Ẑ ≥ Ẑ+U 1+ε holds. The test at line 4 guarantees that this invariant still holds at the end of the algorithm. Therefore Ẑ ≥ Z (1+ε) .

While inspired by Depth First Branch and Bound (DFBB) that provides polynomial space complexity, this algorithm behaves differently from it. In DFBB, for a fixed order of exploration, when the local bound used for pruning (here U b(t)) is tighter, less nodes are explored. This property is lost in Z * ε . Indeed, it is easy to imagine a scenario where a tight bound U b(t) will lead to more nodes being explored than using a weaker U b (t): imagine that search has started and collected a mass Ẑ = 1 and U = 0 for either bounds. Then comes a subtree of small size for which U b(t) = ε while U b (t) ε. This subtree will be pruned by U b(t) leading to U = ε but instead will be enumerated with U b (t) preserving U = 0. In this context, the algorithm using the tight U b(t) is not allowed to prune anymore in the immediate future: if the forthcoming leaves all have very small probability mass, it will be forced to visit all of them while the algorithm using U b (t) preserved some margin and may be able to skip a significant fraction of them.

Indeed, similarly to what happens with the α-β algorithm [START_REF] Knuth | An analysis of alpha-beta pruning[END_REF], the order in which leaves are explored may have a major effect on the algorithm efficiency. Let us assume that we have a perfect U b(t) and that the leaves of the tree have exponentially decreasing mass of probability, the i th visited leaf having a mass of ε i-1 , ε < 1 (such an extreme distribution of probability mass may seem unlikely, but corresponds to linearly increasing energies). In this case, the first leaf bears more mass than all the rest of the tree and the Z * ε algorithm would visit just one leaf. If the inverse ordering of leaves is assumed, the algorithm will have to explore all leaves. It seems therefore important to collect highest masses first. The polynomial space complexity of DFBB comes however with strict constraints on the order of exploration of leaves and best-first algorithms that could overcome this restriction would lead to worst-case exponential space complexity. Interesting future work would be to use the recent highly flexible any-space Branch and Bound algorithm HBFS [START_REF] Allouche | Anytime hybrid best-first search with tree decomposition for weighted csp[END_REF] to improve the leaf ordering within bounded space.

However, contrary to what happens with optimization, even an exact upper bound and a perfect ordering does not guarantee that only one leaf needs to be explored. If we instead assume a totally flat energy landscape, with all leaves having the same energy, Z * ε will have to explore a 1 1+ε fraction of the leaves just to accumulate enough mass in Ẑ to prune.

Overall, it is important to realize that Z * ε needs to achieve two goals: 1. collect probability masses on a potentially very large number of complete assignments to compute a suitable approximation 2. exploit its upper bound to prune the largest possible part of the tree

The first goal could be achieved by existing algorithms producing an exhaustive list of the m-best assignments [START_REF] Dechter | Search algorithms for m best solutions for graphical models[END_REF] or all assignment within a threshold of the optimum (a service that any DFBB-based optimization system provides for free). These algorithms use bounds on the maximum probability instead of the total probability mass which leads to stronger pruning and potentially higher efficiency than Z * ε but do not provide any guarantee since the number m of assignments that would need to be enumerated to provide ε-approximation is unknown.

Because a potentially very large number of probability masses need to be collected, a very fast search is required. To accelerate it, we equip Z * ε with a very simple form of "on the fly" caching: at any node during the search, we eliminate any variable which is either assigned or of bounded degree as proposed initially for optimization [START_REF] Larrosa | Boosting search with variable elimination[END_REF], but using sum-product variable elimination [START_REF] Dechter | Bucket elimination: A unifying framework for reasoning[END_REF]. This caches all the influence of the eliminated variable in a temporary (trailed) potential function. This means that the leaves of the search tree will be sub-problems with bounded tree-width that may represent an exponential number of assignments. This naturally makes Z * ε related to the vec weighted counting algorithm, an anytime MRF counter based on w-cutsets (vertex cutset which if assigned leave a w-tree) and variable elimination over w-trees [START_REF] Dechter | Bucket elimination: A unifying framework for reasoning[END_REF].

The second goal is to prune the largest possible part of the tree search. However, since the first goal requires a very fast search algorithm, using a powerful but computationally expensive bound is probably doomed to fail. For this reason, we have considered simple fast incrementally updated upper bounds by borrowing recent optimization bounds [START_REF] Cooper | Soft arc consistency revisited[END_REF] which are known to work well in conjunction with Depth First Search.

Bounds for guaranteed counting

For any MRF, we define a first upper bound on Z denoted by U b 1 .

Z ≤ U b 1 =   φ S ,|S|<2 t∈D S exp (-φ S (t))   •   φ S ,|S|≥2 exp -min t∈D S φ S (t)   Proof. By definition, we have that Z = t∈D X   φ S ,|S|<2 exp (-φ S (t)) • φ S ,|S|≥2 exp (-φ S (t))  
Trivially, exp (-φ S (t)) ≤ max t∈D S (exp (-φ S (t))) = exp (-min t∈D S φ S (t)) (by monotonicity). Applying this to the right term above, and exploiting the fact that this term now does not depend on t, we get that:

Z ≤   t∈D X φ S ,|S|<2 exp (-φ S (t))   •   φ S ,|S|≥2 exp -min t∈D S φ S (t)  
Since the set {φ S , |S| < 2} contains only unary or constant functions, distributivity allows to swap sum and product and the result follows. Notice that this bound can be computed in linear time.

This bound can be strengthened by selecting a subset of all pairwise potentials in Φ defining a partial spanning k-tree T ⊂ Φ. By applying a sum-product non serial dynamic programming [START_REF] Dechter | Bucket elimination: A unifying framework for reasoning[END_REF] on T = T ∪ {φ S ∈ Φ : |S| < 2}, we can obtain the exact Z T for this sub-MRF in polynomial time. We can multiply Z T by ( φ S ∈Φ\T exp(-min t∈D S φ S (t))) and get a tighter upper bound on Z which we denote U b T :

Z ≤ U b T =   t∈D X φ S ∈T exp (-φ S (t))   Computed using non serial dynamic programming •   φ S ∈Φ\T exp -min t∈D S φ S (t)  
Proof. The proof is essentially similar to the previous one, and obtained by just replacing the set {φ S , |S| < 2} and its complement set {φ S , |S| ≥ 2}, defining the ranges of the products by the sets {φ S ∈ T } (and its complement respectively). The first item can be simply computed in O(nd 2 ) time using non serial dynamic programming.

These bounds alone are very weak. To further strengthen them, we reformulate the MRF using soft arc-consistencies [START_REF] Cooper | Soft arc consistency revisited[END_REF] on its energy representation [START_REF] Allouche | Computational protein design as an optimization problem[END_REF]. Soft arc consistencies essentially shift energy from pairwise potential functions to unary potential functions and eventually to the constant potential function φ ∅ while preserving equivalence. The result of this is an equivalent MRF (defining the same distribution) with increased unary and constant φ ∅ potential functions and pairwise potential functions that satisfy min t∈D S φ S (t) = 0. Besides strengthening the bounds, it removes the need to compute the right term which is always equal to 1. U b 1 and U b T can then be computed in O(nd) and O(nd 2 ) instead of O(ed 2 ) (extension to non pairwise potentials would require the use of partial k-trees instead of trees and change the d 2 into a d k+1 ).

In the rest of the paper we consider spanning trees and try both Existential Directional Arc Consistency (EDAC) and Virtual Arc Consistency (VAC) [START_REF] Cooper | Soft arc consistency revisited[END_REF] as possible ways of strengthening U b 1 and U b T .

Experimental evaluation and comparison

To evaluate the ability of the Z * ε algorithm to provide guaranteed deterministic approximations to Z, we implemented it on the top of the open source toulbar2 solver 4 .

The variable and value ordering used are the default weighted-degree and last conflict variable ordering and the existential support value-ordering [START_REF] Cooper | Soft arc consistency revisited[END_REF]. We enforce EDAC at the root node and during search as usual for optimization. When VAC is used, it is only enforced at the root node because of its computational cost. Instead of the k-way branching described in Algorithm 1, we use a binary branching that either includes or reject a chosen value a at each branching decision. At each node, all variables of degree ≤ 2 are eliminated. The upper bound U b T uses a fixed maximum spanning tree with maximum sum of mean cost after enforcing arc consistencies at the root node. Our implementation is limited to pairwise potentials.

We compared it to different exact weighted counting approaches in terms of efficiency and quality of our guaranteed approximation. Four different exact counters have been considered. The first one is the already described vec exact counter [START_REF] Dechter | Bucket elimination: A unifying framework for reasoning[END_REF]. The second one is the exact SAT based weighted counting tool cachet [START_REF] Sang | Solving bayesian networks by weighted model counting[END_REF] 5 . cachet relies internally on the Zchaff SAT solver to enumerate models with non zero weight and uses context-sensitive independence to cache intermediate counts. We also used the ace 3.0 compiler [START_REF] Chavira | On probabilistic inference by weighted model counting[END_REF], using the UAI competition executable provided in the ace distribution (always using a pseudo-random generator seed of 0). ace computes a tree decomposition and based on the obtained width may either perform tabular variable elimination or encode to CNF and compile in d-DNNF using c2d. We also tested the recent minic2d Sentential Decision Diagram (SDD) compilation package [START_REF] Oztok | A top-down compiler for sentential decision diagrams[END_REF]. SDD are more constrained than d-DNNF and may therefore lead to larger compiled forms than d-DNNF, but since we do not need a compiled form and just the value of Z, we used the -W option of minic2d that performs weighted counting without compilation hoping to trade space for time. minic2d relies on its own internal SAT solver which is provided as a compiled binary in the distributed minic2d package. Because some of the compared solvers (vec, cachet) provide only a double floating representation of Z (or its logarithm), all software has been used in double floating point mode.

All executions have been performed on one core of an Intel Xeon CPU E5-2680 v3 @ 2.50GHz (a Q4 2014 cpu) with a limit of 60 GB on RAM usage.

MRF to #SAT encoding

If ace uses its own internal optimized MRF to SAT encoding, both cachet and minic2d require specific SAT encoding. Exact #SAT weighted counters use weighted literals and define the weight of a model as the product of the weights of all literals which are true in the model. They therefore rely on multiplicative potentials exp(-φ S (t)).

To transform an MRF into a literal-weighted CNF formula with a weighted count equal to the partition function, we use the ENC1 encoding of [START_REF] Chavira | On probabilistic inference by weighted model counting[END_REF], originally described in [START_REF] Darwiche | A logical approach to factoring belief networks[END_REF]. This encoding is the CNF version of the so-called local polytope-based ILP encoding introduced in [START_REF] Schlesinger | Syntactic analysis of two-dimensional visual signals in noisy conditions[END_REF] for MRFs and [START_REF] Koster | Frequency assignment: Models and Algorithms[END_REF] for weighted CSPs [START_REF] Cooper | Soft arc consistency revisited[END_REF]. For each variable i ∈ X, we use one proposition d i,r for each value r ∈ D i . This proposition is true iff variable i is assigned the value r. We encode At Most One (AMO) with hard clauses (¬d i,r ∨ ¬d i,s ) for all i ∈ X and all r < s, r, s ∈ D i , as well as At Least One (f) with one hard clause ( r d i,r ) for each i. These clauses ensure that the propositional encoding allows exactly one value for each variable in each model. For each potential φ S , and each tuple t ∈ D S , we have a propositional variable p S,t . For non-zero energies φ S (t), we have the literal p S,t with weight exp(-φ S (t)). This represents the multiplicative potential to use if the tuple t is used. ¬p S,t is instead weighted by 1, the identity for multiplication. For every variable i ∈ S, we have a hard clause (d i,t[i] ∨ ¬p S,t ). These clauses enforce that if tuple t is used, its values t[i] must be used. Then, for each variable i ∈ S and each value r ∈ D i , we have hard clauses (¬d i,r ∨ t∈D S ,t[i]=r p S,t ) that enforces that if a value r ∈ D i is used, one of the allowed tuples t ∈ D S such that t[i] = r, w S (t) < k must be used.

It is interesting to notice that for pure Constraint Satisfaction Problems (MRFs having only 0/∞ potentials), it is known that Unit Propagation (UP) on this encoding enforces arc consistency in the original CSP [START_REF] Bacchus | GAC via unit propagation[END_REF].

We apply obvious optimization steps, explicitly forbidding local assignments with zero mass (sources of determinism). This encoding can be directly fed into minic2d. Large problems however could not be encoded because minic2d only allows to express weights in a one-line list of maximum 100,000 chars in length 6 .

In cachet, weighted literals l are either such that l and l receive a mass of 1 that has no effect on final mass, or such that the weights of a variable and its negation sum to 1. This is sufficient and convenient to express Bayesian nets because of their local normalization constraint. For arbitrary MRFs, for every p S,t corresponding to a mass m = exp(-φ S (t)) we introduce another propositional variable n S,t with weights m (positive) and 1 -m (negative) and a simple implication clause p S,t → n S,t . This extra variable is connected to the rest of the problem only through this clause and can therefore easily be eliminated, leading to a multiplicative factor m in models where p S,t is true and 1 = m + (1 -m) in models where p S,t is false, as required.

Benchmarks

Two types of benchmarks have been used. The first type of benchmark is made of instances of partition function computation appearing as sub-problems of binding affinity computations on molecular systems defined by a protein interacting with a peptide or an amino-acid. The 3D model of these molecular systems were derived from crystallographic structures of the proteins in complex with their ligands, deposited in the Protein Data Bank. Missing heavy atoms in crystal structures as well as hydrogen atoms were added using the tleap module of the Amber 14 software package [6]. The molecular all-atom ff14SB force field was used for the proteins and the ligands (peptides and amino-acids). The molecular systems were then subjected to 1000 steps of energy minimization with the Sander module of Amber 14. Next, a portion of the proteins including amino-acids at the interface between the protein and the ligand as well as surrounding amino-acids with at least one atom within 8 to 12 Å (according to the molecular system) of the interface was selected. 7To evaluate the effect of the strength of the upper bound on the algorithm efficiency, we applied Z * ε with ε = 10 -3 on a series of 349 systems using the four different bounds. For each system, the most complex partition function, defined on the compound system, is computed. The largest problem has 22 variables and the largest domain size is 34. The gap between our two bounds and the guaranteed approximation of Z determined by Z * ε is shown in Figure 2. The bound U b T is clearly stronger than U b 1 , as expected.

Instances

Log(Ub)-Log(Z)

••••••• • • •• •• • • •• • • • • • • • ••• • •• • • • • ••••• ••• • • • • • •• • • • ••• • • • ••• • • • • • ••• • • • • • • ••••• • • • • • • • • • • • • • • • • • • • • • • • • • •••••••••••• • • • ••• • • • ••• • • • • • • • •• •• • • • • • • • • • • • • • • • • •••••••• • • • • ••• • •• • • • • •• • • • • • • • • • • • • • •• •• • • • • • • •• • ••••• • ••••••••• • • • • • • ••••• • • •• •••••••• • ••••••••••• • •• • •• •• • • •• • • •• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • •• •••• • • • • •• • • • •• • • • • • • • • • • • • • • • • •• • • • • • • • • • • • Ub T Ub 1
Fig. 2. Gap to Z: we represent log(U b) -log( Ẑ) at the root node for U b1 and U bT using EDAC tightening (only negligible difference with VAC on these problems). The instances on the x axis are sorted in increasing gap size for the strongest U bT bound.

We then compare the running times of Z * ε using these 4 bounds in a cactus plot in Figure 3. The best bound in terms of run-time is the lightest U b 1 +EDAC bound confirming that stronger, thus more expensive, bounds may quickly become counter productive.

We represent the same information with the fastest U b 1 +EDAC and two of the three exact counting tools in Figure 4. We omit ace and minic2d. Indeed, ace was able to solve only 17 problems within the time limit and failed on all remaining problems with a memory exception (despite the explicit allocation of 60GB to the JAVA machine). minic2d was instead unable to model 294 systems out of the 349 because of its previously mentioned limitation on the length of the weight line. On the remaining 55 problems, minic2d solved 7 problems in less than one hour. In the rest of the experiments we therefore use the U b 1 +EDAC upper bound which seems the most efficient bound. To see how the Z * ε algorithm performs on other types of problems, we used instances extracted from UAI and PIC'2011 challenge instances (PR task) 8 that use only pairwise potentials as a second set of benchmark. Using the same value of ε = 10 -3 , we again compared Z * ε with vec, cachet, ace and minic2d. The results clearly show there is no single winner: except for cachet which is always dominated by one of the other solvers, each algorithm may outperform others. Specifically, the Z * ε algorithm, despite its lack of sophisticated caching technology, is able to outperform its competitors in various cases. Nevertheless, minic2d outperformed Z * ε on the Grid category (probably because of the combination of Boolean variables and relatively small treewidth), itself outperformed by vec and further outperformed by ace.

Conclusion

Existing solvers providing deterministic guarantees for partition function computation exploit two sources of efficiency. This first one is caching of local counts based on context-sensitive independence [START_REF] Sang | Solving bayesian networks by weighted model counting[END_REF], related to tree-decomposition. The other one is determinism i.e., , the existence of zero probability assignments allowing to prune zero probability mass sub-trees during search. This second source of efficiency will provide significant speedups only when a significant fraction of the search space has 0 probability. Such distributions have very low entropy.

In this paper, motivated by the computation of statistical estimate of affinity between bio-molecules, we have proposed to build upon existing optimization technology to provide a new source of pruning for partition function computation with deterministic guarantees: a branch and bound-based schema equipped with upper bounds derived from soft local consistencies. As existing SAT-based exact approaches, our algorithm exploits determinism and a much simpler and less powerful form of caching than those based on tree-decomposition. It is however able to prune regions of proven negligible mass of probability and is therefore able to exploit relatively low entropy distributions having a much wider support, including those with no determinism. The resulting algorithm offers an adjustable deterministic guarantee on the quality of the computed partition function and, despite its limited caching strategy, may already offer interesting speedups compared to exact solvers.

Z * ε includes two crucial ingredients to quickly gather large number of probability masses: pruning based on very fast incremental upper bounds derived from optimization bounds and on-the-fly sum-prod elimination. An important point is that these ingredients can be easily injected into existing SAT-based counters, including knowledgecompilation based counters using SAT-solver traces. This could be achieved by defining counting upper bounds from existing Max-SAT bounds. These bounds have already been related to soft arc-consistency bounds [START_REF] Larrosa | Resolution in Max-SAT and its relation to local consistency in weighted CSPs[END_REF][START_REF] Bonet | Resolution for max-sat[END_REF][START_REF] Larrosa | A logical approach to efficient max-sat solving[END_REF]. This should extend their range of application to guaranteed approximate probabilistic inference on problems with limited or no determinism.

From an affinity computation point of view, the next step is now to evaluate the actual empirical quality of the association constant estimation provided by the computed ratio of partition functions. Beyond algorithmic approximations, the modeling may also have important effects on the estimated value based on different rotamer discretizations, relative positions of molecules in the complex or weights of different contributions in the energy function. To pursue this target, we intend to use available databases that provide experimental values of the association constant of various protein-ligand complexes following various mutations on one of the partners. To keep the modeling to a reasonable level of complexity, this will be preferably achieved on protein-protein complexes.

Fig. 1 .

 1 Fig. 1. A local view of a protein with a backbone and two acid side-chain reorientations (rotamers) for a given amino-acid (L = Leucine). A typical rotamer library for another amino-acid is shown on the right (ARG = Arginine).

Fig. 3 .

 3 Fig.3. Cactus plot of the running times of Z * ε with our four bounds. A point at (x, y) indicates that the number of solved problems is x if a deadline of y seconds is imposed on each resolution.

Fig. 4 .

 4 Fig. 4. Cactus plot of the running times of Z *ε with the U b1+EDAC bound together with cachet and vec. A point at (x, y) indicates that the number of solved problems is x if a deadline of y seconds is imposed on each resolution.

Table 1 .

 1 Time results for UAI/PIC'2011 instances. Three different categories are represented: Boltzmann machines (rbm) with attractive (ferro) and non attractive coupling, Grids, and graph problems. Running-times are given in seconds. M : Memory Out (60GB), T : Time out (1h). Bold is best.

	Instance	Z * ε	minic2d ace	vec cachet
	smokers 10 < 0.01	0.663 < 0.01 < 0.01 0.264
	smokers 20 < 0.01 312.825 < 0.01 < 0.01 332.168
	rbm 20	7.46	T	3.376 16.17 941.14
	rbm ferro 20 < 0.01	T	3.24	16.18 893.84
	rbm 21	13.75	T	6.854 34.35 2041.64
	rbm ferro 21 < 0.01	T	6.868 34.17 1975.78
	rbm 22	33.75	T	14.411 72.78	T
	rbm ferro 22 < 0.01	T	14.418 72.43	T
	grid10x10.f10 66.32	< 0.01 < 0.01 < 0.01	T
	grid20x20.f10	T	T	13.316 2104.57	T
	grid20x20.f15	T	T	13.665 2099.24	T
	grid20x20.f2	T	T	13.603 2107.92	T
	grid20x20.f5	T	T	13.609 2102.32	T
	GEOM30a 3 < 0.01 < 0.01 < 0.01 < 0.01 2.668
	GEOM30a 4	7.66	78.604 < 0.01 < 0.01 43.77
	GEOM30a 5	67.48 368.361 < 0.01 < 0.01 405.38
	GEOM40 2 < 0.01 < 0.01 < 0.01 < 0.01 < 0.01
	GEOM40 3 < 0.01 < 0.01 < 0.01 < 0.01 3.62
	GEOM40 4	1.42	0.58	< 0.01 < 0.01 616.59
	GEOM40 5	12.67	12.801 < 0.01 < 0.01 828.50
	myciel5g 3	37.2	T	M	T	T
	myciel5g 4	M	T	M	T	T
	myciel5g 5	M	T	M	T	T
	queen5 5 3	367.2	T	83.004 945.20	T
	queen5 5 4 2423.67	T	M	T	T

They may also exploit the fact that counting the number of models of a valid formula is easy. This requires to check for validity, something that modern CDCL solvers do not do anymore.

http://www.inra.fr/mia/T/toulbar2

We thank Jean-Marie Lagniez, CRIL, France for providing us with a patched version of cachet that can be compiled and run without any issue on recent systems.

This parameter could not be changed, being in the non open-source part on minic2d.

Each of these systems requires extensive molecular modeling expertise to be properly defined. We intend to make this benchmark together with the Z * ε implementation available.

http://www.cs.huji.ac.il/project/PASCAL

Acknowledgments : We would like to thank Simon de Givry for his help with toulbar2. We thank the Computing Center of Region Midi-Pyrénées (CALMIP, Toulouse, France) and the Genotoul Bioinformatics Platform of INRA-Toulouse for providing computing resources and support. C. Viricel was supported by a grant from INRA and Region Midi-Pyrénées.