Frank Dondelinger
email: frankd@bioss.ac.uk

Heterogeneous Continuous Dynamic Bayesian Networks with Flexible Structure and Inter-Time Segment Information Sharing

Classical dynamic Bayesian networks (DBNs) are based on the homogeneous Markov assumption and cannot deal with heterogeneity and non-stationarity in temporal processes. Various approaches to relax the homogeneity assumption have recently been proposed. The present paper aims to improve the shortcomings of three recent versions of heterogeneous DBNs along the following lines: (i) avoiding the need for data discretization, (ii) increasing the flexibility over a time-invariant network structure, (iii) avoiding over-flexibility and overfitting by introducing a regularization scheme based in inter-time segment information sharing. The improved method is evaluated on synthetic data and compared with alternative published methods on gene expression time series from Drosophila melanogaster.

Introduction

There has recently been considerable interest in structure learning of dynamic Bayesian networks (DBNs), with a variety of applications in signal processing and computational biology; see e.g. (Robinson & Hartemink, 2009) and [START_REF] Grzegorczyk | Nonstationary continuous dynamic Bayesian networks[END_REF].

The standard assumption underlying DBNs is that time-series have been generated from a homogeneous Markov process. This assumption is too restrictive in many applications and can potentially lead to erroneous conclusions. In the recent past, various research efforts have therefore addressed this issue and proposed more flexible models. (Robinson & Hartemink, 2009) proposed a discrete heterogeneous DBN, which allows for different structures in different segments of the time series, with a regularization term penalizing differences among the structures. [START_REF] Grzegorczyk | Nonstationary continuous dynamic Bayesian networks[END_REF] proposed a continuous heterogeneous DBN, in which only the parameters are allowed to vary, with a common network structure providing information sharing among the time series segments. [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF] proposed an alternative continuous heterogeneous DBN, which is more flexible in that it allows the network structure to vary among the segments. The model proposed in [START_REF] Ahmed | Recovering timevarying networks of dependencies in social and biological studies[END_REF] and [START_REF] Kolar | Sparsistent learning of varying-coefficient models with structural changes[END_REF] can be regarded as a heterogeneous DBN where inference is based on sparse L1-regularized regression (LASSO) of the interaction parameters, and a second L1 regularization term penalizes differences between networks associated with different segments.

Parameter estimation in [START_REF] Ahmed | Recovering timevarying networks of dependencies in social and biological studies[END_REF] and [START_REF] Kolar | Sparsistent learning of varying-coefficient models with structural changes[END_REF] is based on penalized maximum likelihood for fixed regularization parameters, which are optimized using BIC or cross-validation.

In the present paper, we follow (Robinson & Hartemink, 2009), [START_REF] Grzegorczyk | Nonstationary continuous dynamic Bayesian networks[END_REF] and [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF] to infer the network structure, the interaction parameters, as well as the number and location of changepoints in a Bayesian framework by sampling them from the posterior distribution with RJMCMC [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF]. The objective of our paper is to propose a model that addresses the principled shortcomings of the three Bayesian methods mentioned above. Unlike (Robinson & Hartemink, 2009), our model is continuous and therefore avoids the information loss inherent in a discretization of the data. Unlike [START_REF] Grzegorczyk | Nonstationary continuous dynamic Bayesian networks[END_REF], our model allows the network structure to change among segments, leading to greater model flexibility. As an improvement on [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF], our model introduces information sharing among time series segments, which provides an essential regularization effect.

Background

This paragraph summarizes briefly the heterogeneous DBN proposed in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF]. The model is based on the first-order Markov assumption. This assumption is not critical, though, and a generalization to higher orders is straightforward. The value that a node in the graph takes on at time is determined by the values that the node's parents take on at the previous time point, -1, as well as the time series segment. More specifically, the conditional probability of the observation associated with a node at a given time point is a conditional Gaussian distribution, where the conditional mean is a linear weighted sum of the parent values at the previous time point, and the weights themselves depend on the time series segment. The latter dependence adds extra flexibility to the model and thereby relaxes the homogeneity assumption. The interaction weights, the variance parameters, the number of potential parents, the location of changepoints demarcating the time series segments, and the number of changepoints are given (conjugate) prior distributions in a hierarchical Bayesian model. For inference, all these quantities are sampled from the posterior distribution with RJMCMC. Note that a complete specification of all node-parent configurations determines the structure of a regulatory network: each node receives incoming directed edges from each node in its parent set. In what follows, we will refer to nodes as genes and to the network as a gene regulatory network. The method is not restricted to molecular systems biology, though.

Model

Multiple changepoints. Let be the number of observed genes, whose expression values = { ( )} 1≤ ≤ ,1≤ ≤ are measured at time points. ℳ represents a directed graph, i.e. the network defined by a set of directed edges among the genes. ℳ is the subnetwork associated with target gene , determined by the set of its parents (nodes with a directed edge feeding into gene ). The regulatory relationships among the genes, defined by ℳ, may vary across time, which we model with a multiple changepoint process. For each target gene , an unknown number of changepoints define + 1 non-overlapping segments. Segment ℎ = 1, .., +1 starts at changepoint ℎ-1 and stops before ℎ , where = ( 0 , ..., ℎ-1 , ℎ , ..., +1 ) with ℎ-1 < ℎ . To delimit the bounds, 0 = 2 and +1 = + 1. Thus vector has length

| | = + 2.
The set of changepoints is denoted by = { } 1≤ ≤ . This changepoint process induces a partition of the time series, ℎ = ( ( )) ℎ-1 ≤ < ℎ , with different structures ℳ ℎ associated with the different segments ℎ ∈ {1, . . . , + 1}. Identifiability is satisfied by ordering the changepoints based on their position in the time series.

Regression model. For all genes , the random variable ( ) refers to the expression of gene at time . Within any segment ℎ, the expression of gene depends on the gene expression values measured at the previous time point through a regression model defined by (a) a set of ℎ parents (parents /edges) denoted by

ℳ ℎ = { 1 , ..., ℎ } ⊆ {1, . . . , }, |ℳ ℎ | = ℎ , and (b)
a set of parameters (( ℎ ) ∈0.. , ℎ ); ℎ ∈ ℝ, ℎ > 0. For all ∕ = 0, ℎ = 0 if / ∈ ℳ ℎ . For all genes , for all time points in segment ℎ ( ℎ-1 ≤ < ℎ ), random variable ( ) depends on the variables { ( -1)} 1≤ ≤ according to

( ) = ℎ 0 + ∑ ∈ℳ ℎ ℎ ( -1) + ( ) (1)
where the noise ( ) is assumed to be Gaussian with mean 0 and variance ( ℎ ) 2 , ( ) ∼ (0, ( ℎ ) 2 ).

Prior

The + 1 segments are delimited by changepoints, where is distributed a priori as a truncated Poisson random variable with mean and maximum = -2:

( | ) ∝ ! 1l { ≤ } . 1
Conditional on changepoints, the changepoint positions vector = ( 0 , 1 , ..., +1 ) takes nonoverlapping integer values, which we take to be uniformly distributed a priori. There are ( -2) possible positions for the changepoints, thus vector has prior density ( | ) = 1/ ( -2 ) . For all genes and all segments ℎ, the number ℎ of parents for node follows a truncated Poisson distribution with mean Λ and maximum = 5:

( ℎ |Λ) ∝ Λ ℎ ℎ ! 1l { ℎ ≤ } .
(2)

Conditional on ℎ , the prior for the parent set ℳ ℎ is a uniform distribution over all parent sets with cardinality ℎ : (ℳ ℎ |ℳ ℎ | = ℎ ) = 1/( ℎ ). The overall prior on the network structures is given by marginalization:

(ℳ ℎ |Λ) = ∑ ℎ =1 (ℳ ℎ | ℎ ) ( ℎ |Λ) (3)
Conditional on the parent set ℳ ℎ of size ℎ , the ℎ + 1 regression coefficients, denoted by ℳ ℎ = ( ℎ 0 , ( ℎ ) ∈ℳ ℎ ), are assumed zero-mean multivariate Gaussian with covariance matrix ( ℎ ) 2 Σ ℳ ℎ ,

( ℎ |ℳ ℎ , ℎ )=|2 ( ℎ ) 2 Σ ℳ ℎ | -1 2 exp ⎛ ⎝ - † ℳ ℎ Σ -1 ℳ ℎ ℳ ℎ 2( ℎ ) 2 ⎞ ⎠ ,
where the symbol † denotes matrix transposition,

Σ ℳ ℎ = -2 † ℳ ℎ ( ) ℳ ℎ ( ) and ℳ ℎ ( ) is the ( ℎ - ℎ-1 )×( ℎ +1
) matrix whose first column is a vector of 1 (for the constant in model ( 1)) and each ( +1) ℎ column contains the observed values ( ( )) ℎ-1 ≤ < ℎ for all factor gene in ℳ ℎ . Finally, the conjugate prior for the variance ( ℎ ) 2 is the inverse gamma distribution, (( ℎ ) 2 ) = ℐ ( 0 , 0 ). Following [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF], we set the hyper-hyperparameters for shape, 0 = 0.5, and scale, 0 = 0.05, to fixed values that give a vague distribution. The terms and Λ can be interpreted as the expected number of changepoints and parents, respectively, and 2 is the expected signal-to-noise ratio. These hyperparameters are drawn from vague conjugate hyperpriors, which are in the (inverse) gamma distribution family: (Λ) = ( ) = (0.5, 1) and ( 2 ) = ℐ (2, 0.2).

Posterior

Equation (1) implies that

( ℎ | ℎ-1 , ℎ , ℳ ℎ , ℎ , ℎ ) = ( √ 2 ℎ ) -( ℎ -ℎ-1 ) exp - ( ℎ -ℳ ℎ ( ) ℳ ℎ ) † ( ℎ -ℳ ℎ ( ) ℳ ℎ ) 2( ℎ ) 2 ) . ( 4
)
From Bayes theorem, the posterior is given by

( , , ℳ, , 2 , , Λ, 2 | ) ∝ (5) ( 2 ) ( ) (Λ) ∏ =1 ( | ) ( | ) ∏ ℎ=1 (ℳ ℎ |Λ) ([ ℎ ] 2 ) ( ℎ |ℳ ℎ , [ ℎ ] 2 , 2 ) ( ℎ | ℎ-1 , ℎ , ℳ ℎ , ℎ , [ ℎ ] 2 )
where all prior distributions have been defined above.

Inference

An attractive feature of the chosen model is that the marginalization over the parameters and in the posterior distribution of ( 5) is analytically tractable: [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF] for details and an explicit expression. The number of changepoints and their location, , , the network structure ℳ and the hyperparameters , Λ, 2 can be sampled from the posterior distribution ( , , ℳ, , Λ, 2 | ) with RJMCMC. A detailed description can be found in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF].

( , ,ℳ, ,Λ, 2 | ) = ∫ ( , ,ℳ, , , ,Λ, 2 | ) (6) See

Model Improvement

Allowing the network structure to change between segments leads to a highly flexible model. However, this approach faces a conceptual and a practical problem. The practical problem is potential model overflexibility. If subsequent changepoints are close together, network structures have to be inferred from short time series segments. This will almost inevitably lead to overfitting (in a maximum likelihood context) or inflated inference uncertainty (in a Bayesian context). The conceptual problem is the underlying assumption that structures associated with different segments are a priori independent. This is not realistic. For instance, for the evolution of a gene regulatory network during embryogenesis, we would assume that the network evolves gradually and that networks associated with adjacent time intervals are a priori similar.

To address these problems, we propose two methods of information sharing among time series segments. The first method is based on the hierarchical Bayesian model of [START_REF] Werhli | Gene regulatory network reconstruction by Bayesian integration of prior knowledge and/or different experimental conditions[END_REF]. However, rather than sharing information hierarchically -comparing all network structures to a central latent structurewe share information sequentially: a network structure is a priori assumed to be similar to the adjacent ones. The second method uses information from all the other segments to define a prior distribution on the edges for a given segment. We will investigate the relative merits of these two information sharing schemes below.

Sequential information sharing

Denote by := + 1 the number of partitions associated with node , and recall that each time series segment ℎ is associated with a separate subnetwork ℳ ℎ , 1 ≤ ℎ ≤ . We impose a prior distribution (ℳ ℎ |ℳ ℎ-1 , ) on the structures, and the joint probability distribution factorizes according to a Markovian dependence:

( 1 , . . . , , ℳ 1 , . . . , ℳ , ) = ∏ ℎ=1 ( ℎ |ℳ ℎ ) (ℳ ℎ |ℳ ℎ-1 , ) ( ) (7)
Similar to [START_REF] Werhli | Gene regulatory network reconstruction by Bayesian integration of prior knowledge and/or different experimental conditions[END_REF] we define

(ℳ ℎ |ℳ ℎ-1 , ) = exp(-|ℳ ℎ -ℳ ℎ-1 |) ( , ℳ ℎ-1 ) (8)
for ℎ ≥ 2, where is a hyperparameter that defines the strength of the coupling between ℳ ℎ and ℳ ℎ-1 . For ℎ = 1, (ℳ ℎ ) is given by (3). The denominator ( , ℳ ℎ-1 ) in ( 8) is a normalizing constant, also known as the partition function:

( ) = ∑ ℳ ℎ ∈ -|ℳ ℎ -ℳ ℎ-1 | (9)
where is the set of all valid subnetwork structures. If we ignore any fan-in restriction that might have been imposed a priori (via ), then the expression for the partition function can be simplified:

( ) ≈ ∏ ∏ ( ) where ( ) = ∑ 1 ℎ =0 -| ℎ -ℎ-1 | = 1 + - (10) and hence = ( 1 + -) (11)
Inserting ( 11) into (8) gives:

(ℳ ℎ |ℳ ℎ-1 , ) = exp(-|ℳ ℎ -ℳ ℎ-1 |) (1 + -) (12)
It is straightforward to integrate the proposed model into the RJMCMC scheme of [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF]. When

proposing a new network structure ℳ ℎ → M ℎ for segment ℎ, the prior probability ratio has to be replaced by the following one:

(ℳ ℎ+1 | M ℎ , ) ( M ℎ |ℳ ℎ-1 , ) (ℳ ℎ+1 |ℳ ℎ , ) (ℳ ℎ |ℳ ℎ-1 , ) = exp[-(|ℳ ℎ+1 -Mℎ | + | Mℎ -ℳ ℎ-1 |)] exp[-(|ℳ ℎ+1 -ℳ ℎ | + |ℳ ℎ -ℳ ℎ-1 |)] ( 13
)
An additional MCMC step is introduced for sampling the hyperparameters from the posterior distribution. For a proposal move → ˜ with symmetric proposal probability ( ˜ | ) = ( | ˜ ) we get the following acceptance probability:

( ˜ | ) = min { ∏ ℎ=2 exp(-˜ |ℳ ℎ -ℳ ℎ-1 |) exp(-|ℳ ℎ -ℳ ℎ-1 |) ( 1 + -) ( ˜ ) ( 1 + -˜ ) ( ) , 1 } (14)
where in our study the hyperprior ( ) was chosen as the uniform distribution on the interval [0, 10]. Note that the scheme proposed in (Robinson & Hartemink, 2009) can be regarded as a special case of the one we propose, with two simplifications: (1) Changepoints are not allowed to vary between nodes.

(2) The common hyperparameter = ∀ has to be chosen by the user in advance and is not inferred from the data.

Global information sharing

We investigate an alternative scheme, based on ideas presented in [START_REF] Ferrazzi | Population models to learn Bayesian networks from multiple gene expression experiments[END_REF]. Let ℎ ∈ {0, 1} denote the indicator variable for a directed edge from node to node in the ℎth network (i.e. the network corresponding to the ℎth section of the time series), and let ∈ [0, 1] denote the probability that the node pair ( , ) is connected by a directed edge. We assume that for a given node pair ( , ), the edge indicator variables { ℎ } are iid distributed,

( ℎ | ) = ( ) ℎ (1 -) 1-ℎ (15)
with a conjugate beta prior on the parameters :

( ) = Γ( + ) Γ( )Γ( ) -1 (1 -) -1 (16)
where and are hyperparameters. Given the subnetworks ℳ h in all segments h different from the current segment ℎ, the prior probability of the subnetwork in the current segment, ℳ ℎ , is

(ℳ ℎ |{ℳ h } h∕ =ℎ ) = ∏ ( ℎ |{ h } h∕ =ℎ ) (17) ( ℎ |{ h } h∕ =ℎ ) = ∫ ( ℎ | ) ( |{ h } h∕ =ℎ )
where

( |{ h } h∕ =ℎ ) ∝ ({ h } h∕ =ℎ | ) ( ) (18)
We introduce the following sufficient statistics:

ℎ is the number of networks in segments different from the current segment ℎ in which the node pair ( , ) is connected by a directed edge. Conversely, ℎ is the size of the complement set, i.e. the number of networks in segments different from the current segment ℎ without an edge from node to node . Obviously, ℎ + ℎ = -1, and

({ h } h∕ =ℎ | ) = ℎ (1 -) ℎ (19)
Inserting ( 19) and ( 16) into (18) leads to:

( |{ h } h∕ =ℎ ) = Γ( + ℎ + + ℎ ) Γ( ℎ + )Γ( ℎ + ) ( 20
) ℎ + -1 (1 -) ℎ + -1
Inserting ( 15) and ( 20) into (17) yields:

( ℎ |{ h }h ∕ =ℎ ) = Γ( + ℎ + + ℎ ) Γ( ℎ + )Γ( ℎ + ) ∫ ( ) ℎ + ℎ + -1 (1 -) ℎ + ℎ + -1 = Γ( + ℎ + + ℎ ) Γ( ℎ + )Γ( ℎ + ) Γ( ℎ + + ℎ )Γ( ℎ + + ℎ ) Γ( + ℎ + + ℎ + 1) (21)
where we have defined ℎ = 1 -ℎ . Using Γ( + 1) = Γ( ), this expression can be simplified:

( ℎ = 1|{ h }h ∕ =ℎ ) = + ℎ + ℎ + + ℎ (22)
The MCMC scheme is identical to the one described in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF], except that (ℳ ℎ |{ℳ h } h∕ =ℎ ) has to be used as the prior on ℳ ℎ , which is obtained by inserting ( 22) into ( 17). In our study, we have set = = 1, in which case ( ) in ( 16) reduces to the uniform distribution over the unit interval. One can extend this scheme by imposing a hyperprior on and , and sampling these hyperparameters from the posterior distribution with MCMC -this is the subject of our future work.

Setup and Diagnostics

The methods described in this paper have been implemented in R, based on code from [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF], and the software is available from the authors upon request. Our program sets up an RJMCMC simulation to sample the network structure, the changepoints and the hyperparameters from the posterior distribution. As a convergence diagnostic we monitor the potential scale reduction factor (PSRF) [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF], computed from the within-chain and betweenchain variances of marginal edge posterior probabilities. Values of PSRF≤1.1 are usually taken as indication of sufficient convergence. In our simulations, we extended the burn-in phase until a value of PSRF<= 1.05 was reached, and then sampled 1000 network and changepoint configurations in intervals of 200 RJMCMC steps. From these samples we compute the marginal posterior probabilities of all potential interactions, which defines a ranking of the edges in the recovered network. For the synthetic simulation study (see below), the gold standard (i.e. the true interaction network) is known. Therefore, by varying the threshold on the rank, we can construct the Receiver Operating Characteristic, or ROC curve (plotting the sensitivity or recall against the complementary specificity), and the precision-recall or PR curve (plotting the precision against the recall). To assess and succinctly score the network reconstruction accuracy, we follow a three-prong approach and compute three figures of merit that have been widely applied in the literature: the area under the ROC curve (AUROC), the area under the PR-curve (AUPRC), and the true positive rate at a fixed false positive rate of 5% (TPFP5).

Data

Synthetic data

We generated synthetic time series, each consisting of = 10 segments of length 50, as follows. Random networks ℳ ℎ , 1 ≤ ℎ ≤ , are generated stochastically, with the number of edges drawn from a Poisson distribution. Each directed edge from node (the parent) to node (the child) has a weight ℎ that determines the interaction strength, drawn from a Normal distribution. The signal associated with node at time , ( -1), evolves according to the heterogeneous first-order Markov process of equation ( 1). Denote by A ℎ the matrix of all interaction strengths ℎ . To ensure stationarity of the time series, we tested if all eigenvalues of A ℎ had a modulus ≤ 1, and removed edges randomly until this condition was met. The networks ℳ ℎ that generated the time series consisted of 10 nodes, with on average 3 parents per node. To simulate a sequence of networks separated by changepoints, we sampled Δ ℎ from a Poisson distribution and then randomly changed Δ ℎ edges between ℳ ℎ and ℳ ℎ+1 , leaving the total number of edges unchanged. The parameter of the Poisson distribution, which determines the average number of changes between adjacent structures, ℳ ℎ and ℳ ℎ+1 , was varied, as described in more detail in Section 6.1.

Gene expression times course during morphogenesis in Drosophila

We also applied our method to the developmental gene expression time series for Drosophila melanogaster (fruit fly), obtained by [START_REF] Arbeitman | Gene expression during the life cycle of Drosophila melanogaster[END_REF]. Expression values of 4028 genes were measured with microarrays at 67 time points during the Drosophila life cycle, which contains the four distinct phases of embryo, larva, pupa and adult. In our study we concentrated on a subset of 11 genes that regulate muscle development. This dataset has also been used in [START_REF] Guo | Recovering temporally rewiring networks: A modelbased approach[END_REF], [START_REF] Zhao | Inferring gene regulatory networks from time series data using the minimum description length principle[END_REF] and (Robinson & Hartemink, 2009). A paired t-test showed that all differences are significant at the 5% level except for the following. AUROC: HetDBN-GI versus HetDBN-SI, 1 change; AUPRC: HetDBN-0 versus HetDBN-GI, 2 changes; TPFP5: HetDBN-0 versus HetDBN-GI, 1 and 2 changes.

Results and Discussion

Experiments on simulated data

We compared the network reconstruction accuracy of three models: the heterogeneous DBN proposed in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF] (HetDBN-0), the heterogeneous DBN with the sequential information sharing scheme proposed in Section 3.1 (HetDBN-SI), and the heterogeneous DBN with the global information sharing scheme proposed in Section 3.2 (HetDBN-GI). The methods were applied to the synthetic data described in Section 5.1. We repeated the simulations for each experimental setting on 10 independent data instantiations, and scored the network reconstruction accuracy with three separate measures, as discussed in Section 4. We investigated how the average number of changes in network structure between adjacent segments affects the performance. Figure 1 shows boxplots of the score distributions. To test for significance of the discerned trends, we carried out a paired t-test; see the caption of Figure 1, and the supplementary material2 for a table. When there are no changes in the network structure, information sharing results in a considerable performance improvement, and HetDBN-GI outperforms HetDBN-SI. The latter finding is plausible, as HetDBN-GI utilizes information from all the segments, whereas HetDBN-SI only utilizes information from the adjacent segments. When the number of edge changes between segments increases, information sharing achieves a less substantial, yet still significant improvement over HetDBN-0. Also, the performance between the two approaches is inverted, with HetDBN-SI slightly yet significantly outperforming HetDBN-GI. Again, this result is plausible. Larger differences among network structures imply that, per se, less is gained from information sharing. Also, given a segment, a network associated with a remote segment will on average have accumulated a larger number of structure differences than a network associated with a close segment; this explains the superiority of the sequential (HetDBN-SI) over the global (HetDBN-GI) information sharing scheme. To investigate the trend more thoroughly, we reduced the computational costs of the MCMC simulations by reducing the network complexity to 1 target node and 20 potential parents, and keeping the hyperparameters fixed. We then carried out simulations over an extended range of average structure differences. The results are shown in Figure 2 -for space restrictions, we only show the AUROC scores. For small numbers of differences among the network structures associated with different segments, information sharing results in a considerable performance improvement over HetDBN-0. The amount of improvement degrades as the differences among structures increase. For small differences, HetDBN-GI tends to outperform HetDBN-SI. This trend is inverted when the difference among network structures increases. These results thus confirm the patterns found in Figure 1, which have been discussed above.

Gene networks related to morphogenesis in the Drosophila life cycle

The top panel in Figure 3 shows the marginal posterior probability of changepoints during the life cycle of Drosophila melanogaster, inferred with the proposed method HetDBN-SI from the gene expression time series described in Section 5.2. For a comparison, we applied the method proposed in (Ahmed & Xing, Heterogeneous Continuous Dynamic Bayesian Networks 2009), using the authors' software package TESLA. Note that this model depends on various regularization parameters, which were optimized by maximizing the BIC score, as in [START_REF] Ahmed | Recovering timevarying networks of dependencies in social and biological studies[END_REF]. The results are shown in the bottom panel of Figure 3, where the graph shows the L1-norm of the difference of the regression parameter vectors associated with adjacent time points. (Robinson & Hartemink, 2009) applied their discrete heterogeneous DBN to the same data set, and a plot corresponding to the top panel of Figure 3 can be found in their paper. A comparison of these plots suggests that our method is the only one that clearly detects all three morphogenic transitions: embryo → larva, larva → pupa, and pupa → adult. The bottom panel of Figure 3 indicates that the last transition, pupa → adult, is less clearly detected with TESLA, and it is completely missing in (Robinson & Hartemink, 2009). Both our method, HetDBN-SI, as well as TESLA detect additional transitions during the embryo stage, which are missing in (Robinson & Hartemink, 2009). We would argue that a complex gene regulatory network is unlikely to transit into a new morphogenic phase all at once, and some pathways might have to undergo activational changes earlier in preparation for the morphogenic transition. As such, it is not implausible that additional transitions at the gene regulatory network level occur. However, a failure to detect known morphogenic transitions can clearly be seen as a shortcoming of a method, and on these grounds our model appears to outperform the two alternative ones. Top panel: HetDBN-SI, posterior probability of a changepoint occurring for any node at a given time (vertical axis) plotted against time (horizontal axis). Bottom panel: TESLA, L1-norm of the difference of the regression parameter vectors associated with two adjacent time points (vertical axis) plotted against time (horizontal axis). The vertical dotted lines indicate the three morphogenic transitions.

In addition to the changepoints, we have inferred network structures for the morphogenic stages of embryo, larva, pupa and adult. An objective assessment of the reconstruction accuracy is not feasible due to the limited existing biological knowledge and the absence of a gold standard. However, our reconstructed networks show many similarities with the networks discovered by (Robinson & Hartemink, 2009), [START_REF] Guo | Recovering temporally rewiring networks: A modelbased approach[END_REF] and [START_REF] Zhao | Inferring gene regulatory networks from time series data using the minimum description length principle[END_REF]. For instance, we recover the interaction between two genes, eve and twi. This interaction is also reported in [START_REF] Guo | Recovering temporally rewiring networks: A modelbased approach[END_REF] and [START_REF] Zhao | Inferring gene regulatory networks from time series data using the minimum description length principle[END_REF], while (Robinson & Hartemink, 2009) seem to have missed it. We also recover a cluster of interactions among the genes myo61f, msp300, mhc, prm, mlc1 and up during all morphogenic phases. This result is not implausible, as all genes (except up) belong to the myosin family. However, unlike (Robinson & Hartemink, 2009), we find that actn also participates as a regulator in this cluster. There is some indication of this in [START_REF] Zhao | Inferring gene regulatory networks from time series data using the minimum description length principle[END_REF], where actn is found to regulate prm. As far as changes between the different stages are concerned, we found an important change in the role of twi. This gene does not have an important role as a regulator during the early phases, but functions as a regulator of five other genes during the adult phase: mlc1, gfl, actn, msp300 and sls. The absence of a regulatory role for twi during the earlier phases is consistent with [START_REF] Elgar | mef2 activity levels differentially affect gene expression during Drosophila muscle development[END_REF], who found that another regulator, mef2 (not included in the dataset) controls the expression of mlc1, actn and msp300 during early development.

Conclusions

We have proposed a novel heterogeneous DBN, which has various advantages over existing schemes: it does not require the data to be discretized (as opposed to (Robinson & Hartemink, 2009)); it allows the network structure to change with time (as opposed to [START_REF] Grzegorczyk | Nonstationary continuous dynamic Bayesian networks[END_REF]); it includes a regularization scheme based on inter-time segment information sharing (as opposed to [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF])); and it allows all hyperparameters to be inferred from the data via a consistent Bayesian inference scheme (as opposed to [START_REF] Ahmed | Recovering timevarying networks of dependencies in social and biological studies[END_REF]). An evaluation on synthetic data has demonstrated an improved performance over [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF]. The application of our method to gene expression time series taken during the life cycle of Drosophila melanogaster has revealed better agreement with known morphogenic transitions than the methods of (Robinson & Hartemink, 2009) and [START_REF] Ahmed | Recovering timevarying networks of dependencies in social and biological studies[END_REF], and we have detected changes in gene regulatory interactions that are consistent with independent biological findings. We have carried out a comparison between two alternative paradigms of information sharing -global versus sequential -and we have discussed the relative merits and shortcomings.

Figure 1 .

 1 Figure 1. Network reconstruction accuracy, measured with three scoring schemes, as discussed in Section 4. Left panel: AUROC; centre panel: AUPRC; right panel: TPFP5. The boxplots show the distributions of these scores, where the horizontal bar shows the median, the box margins show the 25th and 75th percentiles, the whiskers indicate data within 2 times the interquartile range, and circles are outliers. The grey shading indicates the method. Unshaded boxes: HetDBN-0. Light shading: HetDBN-SI. Dark shading: HetDBN-GI. The numbers on the horizontal axes indicate the average number of network structure changes per node between adjacent time series segments. A paired t-test showed that all differences are significant at the 5% level except for the following. AUROC: HetDBN-GI versus HetDBN-SI, 1 change; AUPRC: HetDBN-0 versus HetDBN-GI, 2 changes; TPFP5: HetDBN-0 versus HetDBN-GI, 1 and 2 changes.

Figure 2 .

 2 Figure2. Network reconstruction accuracy for different methods. The plots show mean AUROC scores (vertical axis) plotted against the average number of network structure changes per node between adjacent time series segments (horizontal axis). Mean values and standard errors were obtained from 10 independent time series.

  Figure 3. Changepoints during morphogenesis in Drosophila melanogaster.Top panel: HetDBN-SI, posterior probability of a changepoint occurring for any node at a given time (vertical axis) plotted against time (horizontal axis). Bottom panel: TESLA, L1-norm of the difference of the regression parameter vectors associated with two adjacent time points (vertical axis) plotted against time (horizontal axis). The vertical dotted lines indicate the three morphogenic transitions.

Appearing in Proceedings of the 27 ℎ International Conference on Machine Learning, Haifa, Israel, 2010. Copyright 2010 by the author(s)/owner(s).

A restrictive Poisson prior encourages sparsity of the network, and is therefore comparable to a sparse exponential prior, or an approach based on the LASSO.

The supplementary material can be found at: http://www.bioss.ac.uk/staff/dirk/Supplements/


Zhaohui Fu

Yinlei Yu

Sharad Malik
email: sharad@princeton.edu

Considering Circuit Observability Don't Cares in CNF Satisfiability *

 that enables the use of the CNF representation even for circuit verification tasks. However, this process loses some information regarding the structure of the circuit. One example of such structural information is the Circuit Observability Don't Cares. Several recent papers [6] [7] [8] [9] [11] [13] have addressed the issue of handling circuit unobservability in CNF-based SAT. However, as we will demonstrate, none of these accurately captures the conditions for use of this information in all stages of a CNF-based SAT solver. In this paper, we propose a broader approach to take such Don't Care information into consideration in a CNF-based SAT solver. It does so by adding certain don't care literals to clauses in the CNF representation. These don't care literals are treated differently at different times during the solution process, much like don't cares in logic synthesis. The major merit of this scheme, unlike other recently proposed techniques, is that the solver can continue to use this don't care information during the learning process, which is an important part of contemporary SAT solvers. We have implemented this approach in the zChaff SAT solver and experiments show that significant performance gain can be obtained through their use.

INTRODUCTION

Boolean Satisfiability (SAT) is probably one of the most well studied combinatorial optimization problems. Researchers have devoted significant effort to developing efficient SAT solvers. For decades, Electronic Design Automation (EDA), and in particular synthesis and verification, has been one of the major drivers for SAT research.

Most SAT solvers work on a CNF representation of a formula. The primary motivation is that the iterative application of the unit clause rule (referred to as Boolean Constraint Propagation or BCP) can be done very efficiently using this data structure, and in fact is the workhorse of almost all SAT solvers. This is not a limitation for using SAT on formulas that arise from circuits as there is a simple conversion for capturing a circuit as a CNF. This involves taking the conjunction of the consistency conditions for each gate in the circuit expressed in CNF form. This is linear in the size of the circuit. However, the disadvantage of using a CNF SAT solver over a circuit based verification technique is the lack of circuit structural information. Such structural information that includes the direction of gates and Circuit Observability Don't Cares (Cir-ODCs), is lost after the circuit is translated into a CNF formula. This information is potentially useful in the solution process. This paper attempts to bridge this gap by providing an approach that takes Cir-ODCs into account during the SAT solving process.

The proposed approach identifies the Cir-ODCs in a given circuit statically, i.e. prior to the start of the SAT solving process. This information is passed to the SAT solver, which in turn adjusts its decision heuristic, BCP and conflict driven procedures based on this additional information. Most importantly, our approach distinguishes itself from other methods of handling observability don't cares by propagating this don't care information during the learning process that is an integral part of modern SAT solvers. This ability is extremely useful because the number of clauses that the SAT solver learns during the search for a solution is usually much larger than the number of original clauses.

While in general, the unobservability condition for a signal is represented as a set of cubes, we (and all other competing approaches) limit ourselves to literals for efficient use in SAT solvers. We also consider an extension of the static scheme presented in the paper that can potentially exploit a larger set of don't cares by allocating them dynamically during the SAT solving process rather than statically. However, in practice, the overhead of managing them generally overweighs their performance advantage.

This paper is organized as follows. Section 2 describes how Cir-ODCs arise in circuits. Our static approach is discussed in detail in Section 3. Section 4 describes the modification of a CNF SAT solver for handling Cir-ODCs. Next, we very briefly explain the dynamic approach for handling Cir-ODCs in Section 5. Experimental results are presented and discussed in Section 6. We give a review of related work in Section 7 and then provide conclusions and future directions in Section 8.

CIRCUIT OBSERVABILITY DON'T CARES

The notion of circuit observability don't cares is as follows. Some signal s, under certain conditions C, no longer affects the outputs of the circuit. These conditions, C, are referred to as observability don't care conditions for s.

We consider circuits with simple gates (AND, OR, NOT, NOR, NAND) as these are typically used as the starting point in either CNF or circuit based SAT solvers due to the ease of deduction. Consider Figure 1 with a simple two input AND gate G, which takes the output value of a logic block, or a logic "cone", C of gates as one of its inputs. The output value of the AND gate G is fixed to be 0 if its input from B is zero, regardless of the output value of C, and hence the entire logic cone of C need not be considered any more if it has only one fan-out that goes into G. Correspondingly in SAT, all the clauses corresponding to C can be ignored in the SAT solver in subsequent search. These clauses are labeled as inactive in [START_REF] Gupta | Dynamic detection and removal of inactive clauses in SAT with application in image computation[END_REF]. We treat these clauses as Cir-ODC clauses and those gate outputs as Cir-ODC variables, since they are unobservable at the outputs of the circuit. Note that the Cir-ODC status of a clause or a variable is a dynamic property in SAT search since it depends on the current assignment during the search.

Ignoring the Cir-ODC clauses in CNF will not change the satisfiability of the original problem. However, it potentially leaves the variables in these clauses under-constrained. This is not an issue if the instance is unsatisfiable, since even with being underconstrained the instance has already been shown to be unsatisfiable. Consider Figure 1 again. When B out = 0, the output of G is 0 regardless of the output value of C. If the original problem is unsatisfiable, the unsatisfiability must not involve any clauses in C. This is because the output of C is no longer being constrained when B out = 0. However, if the instance is satisfiable, then to get a complete assignment, we need to ensure that the final assignment of variables in the Cir-ODC clauses is consistent with the circuit. This is easily accomplished by no longer ignoring the Cir-ODCs once the instance is determined to be SAT. Consider the following Algorithm 1, where v is a variable and V is the set of variables in the CNF.

Algorithm 1 Re-validate the truth assignment.

1: for all variables v ∈ V (CNF) do 2:

if v only appears in ignored Cir-ODC clauses then 3:

Unassign the truth value of v 4:

end if 5: end for 6: Re-solve the CNF instance with the current partial assignment

The main idea behind Algorithm 1 is that the only inconsistency in the satisfying assignment of the CNF formula comes from the clauses that are completely ignored due to Cir-ODC.

A STATIC MECHANISM FOR HANDLING CIR-ODCS IN CNF-SAT

We start with a basic review of the use of Cir-ODCs in SAT solvers. Consider a portion of a typical logic circuit in Figure 2. The lower case letters are the outputs associated with each gate and they also correspond to the variables in the CNF formula of the circuit. The value of a is unobservable if either b or c is 0 since the output f is determined independent of a. We denote that { b, c} (i.e. b+ c) is the Cir-ODC condition to make a unobservable. Thus, if we propagate this Cir-ODC condition to the gates G 2 and G 3 , all clauses associated with these gates become Cir-ODC when this condition is true, i.e. either b or c is set to 0. The backward propagation of Cir-ODC condition only stops at some gate that has fan-out outside the transitive fan-in cone of f .

Total Ordering of Gate Inputs

Consider the following case in Figure 2. First a is set to 0, making both b and c unobservable. Then c is set to 0. Setting c to 0 by itself makes a and b unobservable. There is a problem in considering both a and c as unobservable as at least one of a and c needs to take responsibility for setting f to 0. This can be taken care of by a total ordering on the inputs of a, b and c, i.e. only a lower order signal can appear in the Cir-ODC condition of a higher order one. This issue is discussed in the context of using ODCs in logic synthesis [START_REF] Bartlett | Multilevel logic optimization using implicit don't cares[END_REF]. In fact, any arbitrary ordering π of gate inputs is sufficient to ensure that one gate input is responsible for generating the controlling value at the gate output. For a two input logic AND gate with inputs u and v, we force that v can appear in the Cir-ODC condition of u if and only if π(v) < π(u). We can use the circuit in Figure 2 as an example and order the inputs alphabetically, which results in π(a) < π(b) < π(c). The Cir-ODC conditions for C 1 and C 2 are { ā} and { ā, b} respectively. Gate G 2 has an empty set as its Cir-ODC condition. Clearly, such a static scheme loses some of the Cir-ODC condition at each gate. This is because in a static total ordering of a logic AND gate with inputs u and v, either ū can make v unobservable, or v can make u unobservable, but not both. Thus if we choose either one, say ū is the Cir-ODC condition of v, and in the SAT solving process there is a chance that v = 0 is assigned earlier, which can make u unobservable, but the SAT solver cannot use that since such condition was thrown away in the translation phase just to ensure the total ordering.

Propagation of the Cir-ODC Conditions

As we have mentioned in the previous section, each gate input has a set of don't care literals that captures the condition for it to be unobservable. Each gate also inherits a set of Cir-ODC conditions from its successors. Recall that a Cir-ODC condition is a set of don't care literals. These Cir-ODC conditions are then propagated towards the primary inputs. A Depth First Search (DFS) from the primary inputs is sufficient for such propagation in the circuit. We start by initializing the set of don't care literals, i.e. the Cir-ODC condition, to be empty for each gate. Then we perform a DFS from each of the inputs. The exact recursive DFS procedures are given in Algorithm 2 and 3.

Algorithm 2 Cir-ODC-Gate (Gate g) 1: if g is a primary output then 2: Return / 0; 3: end if 4: C = Universal Set; 5: for all fan out edge [g, h] of g do 6: C := C Cir-ODC-Gate(h) Cir-ODC-Edge([g, h]) ; 7: end for 8: Return C;

Algorithm 3 Cir-ODC-Edge (Edge [g, h]) 1: C = / 0; 2: for all fan in gate i of h, i = g do 3: if π(i) < π(h) then 4: C := C {i = controlling value for h} 5:
end if 6: end for 7: Return C;

We illustrate the Cir-ODC condition propagation using the circuit in Figure 3. Assuming an alphabetical ordering of the inputs, i.e. π(r) < π(s) < π(t) and π(p) < π(q), all logic gates in cone C 1 inherit a Cir-ODC condition {r, p}. Similarly, all gates in cone C 2 inherit {r, s, p}. Gate G 2 has two successors, one in C 1 and the other in C 2 . The Cir-ODC condition that G 2 inherits is the intersection of the two Cir-ODC conditions from them, which is {r, p}, i.e. a normal set intersection.

A Greedy Ordering Heuristic for Gate Inputs

As we have discussed earlier in Section 3.1, a total order of the gate inputs is needed to ensure that at least one input will provide the controlling value of the gate when needed. For each gate, any given total order of the inputs completely determines the total number of don't care literals contributed by this gate. In order to utilize these don't care literals in CNF-SAT, we are interested in finding a good total order such that it gives the most don't care literals in the circuit. This is because a more "popular" don't care literal is likely to cause a larger portion of the circuit to be unobservable when it is set to be true.

We propose a greedy ordering heuristic that works in a similar manner as the propagation of don't care literals discussed in previous section. For each gate that could generate new don't care literals, the maximal set of don't care literals for each input is propagated towards the primary input, disregarding the total order con-straint, i.e. ignoring the π(i) < π(h) condition in Algorithm 3. Consider Figure 3 as an example. Don't care set { s, t, p} is propagated through input pin r; don't care set {r, t, p} is propagated through input pin s; and so on. Obviously, this violates the total ordering constraint, but such propagation gives an estimate of the number of gates that each don't care literal can eventually reach. After the propagation is finished, the gate input whose don't care literal that reaches more gates is assigned a lower order rank. A tie in the order is broken by random. When the circuit is to be translated in to CNF format, the don't care literals for each clause are also ordered in the same way as above, i.e. a don't care literal with larger use appears earlier in the clause. This helps the SAT solver pick the more important literals when it cannot accept all of them because of efficiency reasons.

SOLVING SAT WITH CIR-ODCS

There are several Cir-ODC related modifications that have to be made in a CNF-based SAT solver. These modifications include devising an augmented CNF format for Cir-ODC encoding, changing the BCP process to ignore the inactive Cir-ODC clauses, handling the don't care literals in conflict analysis and conflict driven clause learning and using Cir-ODC condition to guide the decision heuristics.

Format of Cir-ODC-CNF

In a CNF file using DIMACS format [START_REF]Satisfiability suggested format[END_REF], a clause is represented with a line of literals that ends in a 0. To minimize the change in the CNF format, we define our new file format for Cir-ODC-CNF such that each clause also has one line. This can be represented using a regular expression as follows.

Clause = l + 0 d * 0
where l stands for a normal literal with phase (+/-), this is the same as in conventional CNF. l + means that there are 1 or more normal literals. d is a don't care literal and d * means that there may be 0 or more don't care literals. 0s are used as delimiters and a line starting with a leading 0 ends the input. A clause can then be ignored if any of the don't care literals is set to be true. Consider the example in Figure 3 and G 2 has a Cir-ODC condition (don't care set) of {r, p}. Suppose the variable indices are p = 1, r = 2, s = 3,t = 4, x = 5, y = 6, z = 7 and u = 8. Also suppose p occurs more frequently as a don't care literal in the circuit. The following four lines will be generated in the Cir-ODC-CNF file.

5 -8 0 -1 -2 0 6 -8 0 -1 -2 0 7 -8 0 -1 -2 0 -5 -6 -7 8 0 -1 -2 0
Note that in the above example, p and r appear as don't care literals in these clauses. They may appear as regular, i.e. not don't care, literals in other clauses.

Enabling zChaff with Cir-ODCs

BCP in zChaff and other state-of-the-art CNF SAT solvers benefits significantly from the two literal watching scheme [START_REF] Moskewicz | Chaff: Engineering an efficient SAT solver[END_REF]. We will assume familiarity with this and will not be reviewing this for sake of brevity.

With the augmented clause description, the don't care literals are not watched but examined during two literal watching. If a don't care literal is satisfied, the clause is considered to be satisfied. This terminates the search for a new literal to watch and thus can potentially terminate BCP quicker. Essentially the don't care literals have a conditional don't care status. They can only satisfy a clause early without needing to be zero before an implication is derived from the clause. This is completely different from Velev's use of unobservability variables [START_REF] Velev | Encoding global unobservability for efficient translation to SAT[END_REF], in which the don't care literals can no longer be differentiated from the regular ones after the CNF translation.

The checking of don't care literals in a clause is performed in a lazy way. The clauses are not marked inactive or ignored until they are accessed during BCP. This is an advantage over the other approaches of explicitly marking clauses active/inactive, which need an explicit bookkeeping overhead. Consider the example given in Figure 4. The clause shown has literals ā, b, c,..., ē, f , g with don't care literals r, s,t delimited by a 0. Suppose one of the two watch pointers is initially pointing to literal f . When f = 1, this clause is examined in BCP since f = 0 and we need to find a new unassigned literal to watch. The watch pointer then moves to the next literal after f , which is g. If g is assigned to be 0, the watch pointer continues moving right and it reaches the don't care delimiter 0. The watch pointer continues to move right and examines the don't care literals. If any of the don't care literals is found to be true, e.g. s = 1, the clause can be ignored and the BCP for this clause is immediately terminated without scanning the rest of the clause, i.e. ā, b, c,..., ē. Thus, if a clause is ignored (inactive), there will be no implications based on this clause. Note that even though the don't care literals are directly integrated into zChaff's literal pool, they must be differentiated from other normal literals since a don't care literal is never allowed to be watched. In other words, there should be no implication on a don't care literal. When a clause is ignored, we leave the watch pointer unchanged from its position before the clause was ignored. Note that this is slightly different from the original zChaff implementation of the two literal watching scheme. We use Figure 4 as an example. In the ordinary two literal watching scheme, the watch pointer needs to be updated such that it points to the recently discovered true literal, in this example, s. However in our case, s cannot be watched since it is a don't care literal. We leave the watch pointer on f unchanged. This modification does ensure that the original invariant (i.e. each "not yet satisfied" clause has two watch pointers pointing to two different unassigned literals) holds after a backtrack. This is because dlevel(s) ≤ dlevel( f )1 , i.e. s = 1 is assigned at a decision level no later than f = 0. This clause can only resume its "not yet satisfied" status on a backtrack to a decision level d < dlevel(s), by which time f has already been unassigned since d < dlevel(s) ≤ dlevel( f ).

zChaff uses the Variable State Independent Decaying Sum (VSIDS) [START_REF] Moskewicz | Chaff: Engineering an efficient SAT solver[END_REF] decision heuristic. VSIDS is based on the number of occurrences of a literal in the CNF, known as literal count, and this is updated periodically. For each don't care literal in the clause, we also increase its literal count. This has a direct impact on the decision heuristic, where a frequently occurring don't care literal may be chosen and assigned to its controlling value. This, in turn, sets clauses to be ignored. The original SAT problem is thus potentially simplified due to the decrease in the total number of unsatisfied clauses. We believe that the performance gain for the static approach comes partially from the increasing scores of the don't care literals, as they are more likely to be branched on earlier in the search.

Handling this Cir-ODC condition in conflict clause generation is the main contributor in the performance gain, as we will demonstrate in Section 6. Recall that learned clauses are particularly useful in generating implications (i.e. forced assignments) during the search process. However, we would like these implications to be useful, i.e. arise from observable active parts of the circuits. The following algorithm provides for this. The normal literals in the learned conflict clause are obtained in the same way as before using resolution on a set of clauses. The don't care literals of the learned clause are obtained as follows. A set union operation of all the don't care literals from the clauses involved in the resolution captures all the don't care literals for the learned clause. The rationale for this is as follows. The learned clause is satisfied if any of its source clauses is satisfied. These source clauses are satisfied when any of the don't care literals is true. Thus the don't care literal set of the learned clause is the set union of the don't care literals of the source clauses. Consider the following example of resolving two clauses (a + b + c) and (b + d + ē) with sets of don't cares {r, s} and { x, r, s} respectively. The resolved (learned) clause is (a + c + d + ē) with the set of don't cares {r, s, x, s}. It is worth mentioning that in classical SAT if a clause contains a variable in both phases, this clause is a tautology and can be deleted from the formula. However, this is not true for a clause with don't care literals. A variable can appear in both phases in the don't care literals, as either of them may trigger the unobservability of this clause.

For implementation efficiency, when a learned clause is added to the clause database, the don't care literals are ordered according to the non-increasing order of their occurrence count and only a certain amount (proportional to the number of its normal literals) of the don't care literals are added to the learned clause.

Conflict driven clause learning with observability don't care information is a major contribution of our approach. This is not handled by either the circuit based approaches of marking inactive sections of the circuits or clauses (e.g. Gupta et al. [START_REF] Gupta | Dynamic detection and removal of inactive clauses in SAT with application in image computation[END_REF]) or Velev's unobservability variables. In our approach, the observability don't care information can be learned during the searching process. This is significantly different from any other method, which can only make limited use of Cir-ODCs information provided with the original circuit. Learning of the Cir-ODC condition is very important because the majority of the clauses that the SAT solver deals with during the search process are learned conflict clauses; this is especially true for the very hard instances. Learning with Cir-ODC also contributes to the VSIDS scoring as certain key don't care literals will appear in more clauses, which results in an increase in their VSIDS scores.

A DYNAMIC APPROACH FOR CIR-ODC IN SAT

The static approach enforces the input ordering during the translation phase. This unavoidably reduces the total amount of Cir-ODC conditions that a CNF SAT solver can utilize. To overcome this we also considered a more dynamic technique for a SAT solver to encode the Cir-ODC conditions. In the dynamic approach, the assignment of variables during the search process is used to order the inputs of a gate. The inputs are ordered in increasing order of their being set to controlling values. Now the Cir-ODC condition of a logic gate is encoded using a set of tuples instead of just lit-erals. These tuples capture the dynamic ordering information. For the example in Figure 2 the Cir-ODC condition for all clauses in the transitive fan-in of signal a is if b = 0 or c = 0 before a is set to be 0. This is represented by the tuple ({ ā}, { b, c}) in each of the these clauses.

Though the rationale behind Cir-ODC condition propagation and conflict driven clause learning in the dynamic approach is same as the one used in the static approach, the implementation is much more complicated than the static case because of the use of tuple representations. These details are omitted here for sake of brevity. As a result of more complicated operations, it is empirically observed that the overhead of the dynamic approach outweighs any benefit received from them.

EXPERIMENTS

In order to perform the experiments with Cir-ODCs, we need to have the combinational circuit description for the benchmarks. The only two benchmark families we found with such descriptions are iscas85 [START_REF] Brglez | A neutral netlist of 10 combinational benchmark circuits and a target translator in Fortran[END_REF] and itc99 [START_REF] Corno | RT-Level ITC 99 benchmarks and first ATPG results[END_REF]. All test case used whose names start with a c are from the iscas85 family, where c1908 is generated from a ECAT circuit; c2670, c3540, c5315 and c7552 are all generated from ALU and control circuits; and c6288 is taken from a 16-bit multiplier circuit. We use five test cases from the itc99 family, namely b14 from a subset of Viper processor circuit; b15 from a subset of 80386 processor circuit; b20, b21 and b22 are composed of multiple copies of circuits that are similar to b14. b14-opt is the optimized circuit for b14, and similarly for b15-opt, b20-opt, b21-opt and b22-opt. For each test case in the benchmark, we generate a miter circuit [START_REF] Brand | Verification of large synthesized designs[END_REF] that consists of two identical copies of the test circuit. A miter circuit outputs 0 if the two test circuits always output the same values and it outputs 1 otherwise. Table 1 shows the number of don't care literals using both the random ordering and the greedy ordering heuristic. Row 1 of this table corresponds to the miter version of c1908 with 1819 logic gates. There are 4,882 CNF clauses generated using this circuit. If we use a random total ordering of the gate inputs, a total of about eight thousand don't care literals are captured, which gives an average of 1.53 don't care literals per CNF clause generated. However, if we use the greedy ordering heuristic to order the gate inputs, we could identify a total of seventeen thousand don't care literals and this gives an average of 3.46 don't care literals per CNF clause generated. Clearly, the greedy order-ing heuristic is able to capture many more don't care literals than the random ordering. Besides, the greedy ordering heuristic could also order the don't care literals in a single CNF clause according to their importance, i.e. the most frequently appearing literals are placed in front of the others.

Given the numbers in Table 1, we can see that the total number of don't care literals in some circuits is extremely large. Scanning all these literals could be slow. Instead of using all the Cir-ODC condition, we restrict the number of don't care literals for each clause to not exceed a certain threshold, which is usually set to be half of the total number of normal literals in this clause.

All the experiments are conducted on a Dell PowerEdge 700 running Linux Fedora Core 1.0 (g++ GCC 3.3.2) with single Pentium 4 2.8GHz, 1MB L2 cache CPU on 800MHz main bus. Table 2 tabulates the running time of both random and heuristic ordering in the static approach (columns labeled Static Random and Static Heuristic). In addition, to study the benefit of don't care literals through the learning process, we consider the case when their use is turned off during learning (column labeled Static NL). The running time of the dynamic (column labeled Dynamic) approach and the latest zChaff is also presented for comparison. Clearly, the static approach outperforms both the dynamic approach and zChaff. The static approach without don't care literal learning performs only slightly better (sometimes worse, as for c2670) than original zChaff. The main performance gain for the static method comes from the conflict clause learning with Cir-ODC. This is particularly true for large and hard benchmarks as they tend to have a large clause pool, most of which are learned conflict clauses. Don't care literals in these large number of clauses increase the VSIDS score dramatically as well as provide benefit during BCP. The dynamic approach seems to suffer from the huge overhead introduced by the tuple representation and consequent bookkeeping operations.

RELATED WORK

The concept of circuit observability don't cares was first introduced by Bartlett et al. in their work on multilevel logic minimization, where a Boolean network is optimized into a prime, irredundant R-minimal form [START_REF] Bartlett | Multilevel logic optimization using implicit don't cares[END_REF].

Various circuit-based SAT solvers have been implemented. Safarpour et al. have dealt with both Observability Don't Cares and Controllability Don't Cares using both a CNF and circuit based SAT solver. A variable, which is also an input to a logic gate, is marked lazy [START_REF] Safarpour | Managing don't cares in Boolean satisfiability[END_REF] if it is being dominated by others. They then guide the decision heuristic of the SAT solver not to branch on the lazy variables. The laziness then propagates through the circuit. They sometimes also mark an assigned variable lazy to erase some old bad decisions. In addition to the CNF formula, their SAT solver also needs circuit layout information in order to propagate the laziness of variables. The non-branching on lazy variables is a decision heuristic that uses some Cir-ODC information, but in a limited way, since those Cir-ODC related clauses are still active in the sense that they are generating implications, most of which are on lazy variables. Gupta et al. [START_REF] Gupta | Dynamic detection and removal of inactive clauses in SAT with application in image computation[END_REF] implemented a CNF-based SAT solver that identifies and masks unobservable gates during SAT search. A circuit-based modification [START_REF] Ganai | Combining strengths of circuit-based and CNF-based algorithms for a high-performance SAT solver[END_REF] of the zChaff SAT solver to avoid this directly works on the logic level representation of the circuit, it only uses CNF for the learned clauses. Thus, this technique is a hybrid CNF and circuit based SAT solver. Other more circuitbiased SAT solvers also monitor the circuit signals and utilize the correlation of different signals to prune the search space [START_REF] Kuehlmann | Circuit-based Boolean reasoning[END_REF] [START_REF] Lu | A signal correlation guided ATPG solver and its applications for solving difficult industrial cases[END_REF].

Velev uses a different approach, which focuses on the translation from circuits to CNF [START_REF] Tseitin | On the complexity of derivation in propositional calculus[END_REF] and is thus the closest to our approach. This is illustrated by the circuit given in Figure 1. Velev's approach will add an unobservability literal u to every CNF clause generated by the logic block C. Since the value of C is unobservable at the primary outputs if the output x of gate B is 0, Velev's approach adds a new clause (x + u) during the CNF translation. This clause implies u to be 1 when x = 0, which makes all clauses generated by C satisfied. Another clause ( x + ū) is also added to ensure the observability of the logic block C, where x = 1 implies u = 0 and the output of C is observable at gate G. However, it is not clear how the mutual exclusion in exploiting the don't cares is avoided for the multiple inputs of a gate. Also, once the circuit has been translated into CNF, it is impossible to differentiate the unobservability literals from normal literals. Thus all literals are treated equally and the unobservability literals will also be required to generate implications. Velev's approach is also limited since it only considers the Cir-ODC condition from the logic blocks with only one fan-out. In the experiments, the unobservability variable is only introduced for ITE-trees [START_REF] Velev | Encoding global unobservability for efficient translation to SAT[END_REF] with fan-out count of 1. Velev's approach increases the size of the translated CNF file, but does not change its format. So any stand-alone SAT solver could still be used. However, such an implementation may significantly increase both the number of variables and the number of clauses, which add to the burden of the SAT solver.

CONCLUSIONS AND FUTURE DIRECTIONS

In this paper we describe a technique for handling logic circuit unobservability conditions in a CNF SAT solver. The key contribution is that the don't care conditions are captured as don't care literals in the clauses describing the logic circuit. We show how the don't care information is captured from the local unobservability conditions at a gate and then propagated through the circuit. The don't care literals differ from previous attempts to handle ODCs in SAT solvers in the following ways:

-They enable these literals to be treated differently in different contexts. They can enable a clause to be ignored (effectively satisfied) early without needing to be set while deriving an implication during BCP. This is a significant difference from Velev's unobservability literals [START_REF] Velev | Encoding global unobservability for efficient translation to SAT[END_REF], which are indistinguishable from the normal literals in CNF clauses.

-We show how the don't care information is propagated and inherited during conflict driven clause learning without further referencing any circuit description. This is absent in previous approaches.

Experiments show that most circuits generate a large number of don't care literals. Efficient utilization of this information may tremendously speed up the SAT solver by both speeding up implications and guiding the decision heuristic. There is room for further improving the ideas introduced in this paper and some of the possible areas include engineering a more efficient implementation of the dynamic approach, exploiting other total orderings for the static approach, and fine-tuning the Cir-ODC parameters, some of which may be correlated to other parameters used in zChaff.

Figure 1 :

 1 Figure 1: An AND gate G takes as input the only fan-out from a logic cone C.

Figure 2 :

 2 Figure 2: Cir-ODC condition associated with a three input AND gate.

Figure 3 :

 3 Figure 3: Two inputs of gate G 1 converge at G 2 , thus G 2 has two Cir-ODC condition sets coming from the same successor, G 1 .

Figure 4 :

 4 Figure 4: A typical clause in zChaff's literal pool using the static approach.

Table 1 : Total and average number of don't care literals cap- tured by random ordering and greedy ordering heuristic

 1

	bench	# of	# of	Random	Heuristic
	mark	gates	clauses	# of dc	dc / cls	# of dc	dc / cls
	c1908	1819	4882	8K	1.53	17K	3.46
	c2670	2608	6859	14K	2.10	134K	19.49
	c3540	3411	9327	9K	1.00	44K	4.77
	c5315	4916	14002	91K	6.50	389K	27.78
	c6288	4897	14593	4K	0.29	77K	5.25
	c7552	7338	19848	121K	6.11	600K	30.22
	b14	20057	58594	1490K	25.44	3422K	58.39
	b14-opt 11217	35510	732K	20.64	1787K	50.33
	b15	17669	53468	2277K	42.58	4670K	85.46
	b15-opt 14979	47862	1948K	40.69	3880K	81.06
	b20	40399 118351	7903K	66.78 13801K 116.61
	b20-opt 24949	79389	4690K	59.08	7929K	99.88
	b21	41089 120601	7955K	65.96 14136K 117.21
	b21-opt 25303	80157	4678K	58.36	8131K 101.44
	b22	59849 175632 19059K 108.52 29452K 167.69
	b22-opt 36119 114122 10274K	90.03 15777K 138.25

Table 2 : Running time (seconds) of the static with random and heuristic ordering, static without using Cir-ODCs in learning (Static NL) and dynamic approaches compared to zChaff.

 2

	bench	zChaff		Static Ordering		Dynamic
	mark		Heuristic	Random	NL	Ordering
	c1908	0.72	0.78	0.78	0.49	1.95
	c2670	0.56	0.49	0.51	1.07	0.63
	c3540	18.75	11.21	10.45	17.42	15.63
	c5315	6.36	4.96	5.51	6.78	7.79
	c7552	23.52	12.63	13.41	18.95	25.76
	b14	6620.57	4079.24	3989.55	6037.52	10622.50
	b14-opt	4602.08	3745.72	4005.39	4683.84	4735.81
	b15	152.39	81.20	129.08	150.24	171.61
	b15-opt	161.14	114.02	159.81	160.28	171.69
	b20	20225.20 11014.80 17810.40 25187.20 > 24 hours
	b20-opt 33067.40 26662.70 28580.20 32039.30 > 24 hours
	b21	19625.40 10220.60 17168.30 19034.60 > 24 hours
	b21-opt 34711.90 12109.10 28137.70 33093.90 > 24 hours
	b22	22625.00 15074.50 15085.00 22032.40 > 24 hours
	b22-opt 34224.90 17811.10 28579.70 34948.20 > 24 hours

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

dlevel(v) refers to the decision level in the search tree that variable v is assigned a value.

* This work is supported by grant 2002-TJ-1025 from SRC Corporation and grant 2001-111DT from Intel Corporation. † Z.