Jie S Hu

Feihui Li

Vijay Degalahal

Mahmut Kandemir

N Vijaykrishnan

Mary J Irwin

Compiler-Directed Instruction Duplication for Soft Error Detection *

In this work, we experiment with complier-directed instruction duplication to detect soft errors in VLIW datapaths . In the proposed approach, the compiler determines the instruction schedule by balancing the permissible performance degradation with the required degree of duplication. Our experimental results show that our algorithms allow the designer to perform tradeoff analysis between performance and reliability.

Introduction

Transient faults due to soft errors are becoming an important concern in current computing systems. Consequently, techniques to improve the reliability of application execution in VLIW processors are of interest. In this work, we focus on detecting the transient errors in the VLIW datapath through complier directed instruction duplication. Since the compiler has full control of instruction scheduling in VLIW architectures, the duplication is controlled at compilation time. This has a distinct advantage over runtime techniques as the compiler can balance power, performance, and code size requirements when performing the duplication. Prior work [START_REF] Holm | Low cost concurrent error detection in a vliw architecture using replicated instructions[END_REF][START_REF] Bolchini | A software methodology for detecting hardware faults in vliw data path[END_REF] on soft error detection in VLIW architectures fully duplicates instructions and includes additional comparison instructions to check the results of the duplicate executions. Our experiments reveal that full duplication can be costly from both performance and energy perspectives. In fact, while full duplication is essential in safety-critical applications, for large segments of embedded markets, factors such as power and performance will continue to be as important as reliability.

Compiler-Directed Instruction Duplication

The main idea behind our algorithm is to fill empty execution slots (cells) with duplicate instructions under a performance bound. The instruction duplication mechanism in our approach is supported by an hardware enhancement for efficient result comparison that significantly improves our full-duplication scheme with a performance loss of 10% compared to 70.9% reported in [START_REF] Bolchini | A software methodology for detecting hardware faults in vliw data path[END_REF]. We implemented our instruction duplication approach within the Trimaran framework, which is an integrated compiler/simulation environment. Our results clearly show that our algorithms can be used for improving reliability against soft errors under performance bounds. The following subsections present the architecture and compiler support needed to implement the scheme.

Architectural Support

Architectural support is required for associating a particular instruction with its duplicate. We maintain a integer/floating-point register value queue (IRVQ/FRVQ) and a load/store address queue (LSAQ) for verifying the outputs of instructions (see Figure 1). When a register-output instruction (generates register value) completes, it writes the value into the output register as well as in the RVQ. When the duplicate of this instruction completes its execution, its results are compared with the contents of the RVQ entry associated with the original. A single bit (B1) is used for identifying the instruction as either original (1) or duplicate (0), and the ordering requirement is achieved by appropriate labeling of this bit based on the relative ordering between the original and the duplicate in the instruction schedule. There is an additional bit (B2) associated with each instruction indicating the required action to the RVQ. If the original instruction precedes its duplicate, B2 bit of the original instruction is set. In the duplicate instruction, B2 is also set and B1 is reset to 0. After finishing execution, the original instruction writes to the tail slot of the RVQ, and the duplicate compares its output with this corresponding entry in the RVQ. If both the original and duplicate instructions are scheduled at the same cycle, both B2 bits are reset. These two instructions compare their outputs immediately after execution without invoking any Write/Read to the RVQ.

Compiler Support

We see the compiler-scheduled code as a twodimensional table, where the rows correspond to schedule cycles and the columns represent functional units. Each entry in this table is referred to as a cell. An empty cell represents an idle cycle. Duplication range of an instruction i is the range of cycles within which a duplicate of the instruction can be scheduled. This range is determined by the instructions that i depends on as well as the instructions that overwrite the register read by i. More specifically, we cannot schedule the duplicate before the source operands for the instruction are ready. Similarly, we need to schedule the duplicate before any of the source operands are overwritten (see Figure 2).

Duplication under Performance Bound

For duplicating under the presence of a performance bound, given that the original (i.e., without any duplicates) schedule length is C cycles, our objective is to maximize the number of instructions with duplicates so that the new schedule length, C , is smaller than or equal to (1 + f ) × C, where f is an input parameter, set by the user, that enforces the tolerable limit in performance degradation. One can conduct a reliability-performance tradeoff analysis by varying the value of f . Basically, this algorithm considers each instruction in turn, identifies its duplication range, and cre- ates a duplicate for it if doing so does not cause the schedule length to exceed (1 + f ) × C.

Results

We evaluate the proposed technique using both compiletime and run-time metrics. The compile-time reliability metric is the fraction of instructions with duplicates, or FWD for short. Intuitively, more duplication should result in more reliable execution assuming a uniform distribution of potential soft errors. The run-time reliability metric, on the other hand, is the fraction of instructions executed with their duplicates. As for the performance metric, the compile-time one is the percentage increase in the static schedule length, and the run-time one is the percentage increase in the program execution time with respect to the case without any duplication. Throughout our discussion we use the term original schedule to refer to the version without any duplication.

Compile-Time Evaluation

Our experimental results show that the average FWD values across our benchmarks are 86.8% and 41.4% for the full-duplication scheme and the scheme with zero performance penalty. The encouraging news is that our approach is able to duplicate more than 40% of the instructions without an increase in the original schedule cycles. Note that, the full-duplication scheme also incurs an average increase of 42.2% in the original schedule length. These FWD and schedule length results clearly illustrate the tradeoffs between reliability and performance. Figure 3 shows the FWD values under different bounds on the percentage increase in the original schedule length. For example, a point x% on the x-axis indicates that we can tolerate at most x% increase in the original schedule length. We see that, as expected, the FWD value increases as we relax the performance bound. The curves in Figure 4 shows the actual increase in the original schedule length. The curve labeled "reference" represents the worst case scenario where the schedule for each basic block is increased by x%. We see that the actual schedule length increases are always lower than the maximum allowable limit.

Run-Time Evaluation

Our experimental results show that the runtime behavior matches very well the static behavior extracted by the compiler. Figure 5 gives the percentage of instructions that are executed with their duplicates under different maximum allowable increases in the original schedule length. We see that these values are similar to the corresponding values given in Figure 3 and4, indicating that the compile-time reliability estimate is reasonable. The graph in Figure 6 shows the percentage increase in the original execution cycles. Note that these values are lower than the corresponding bounds for all benchmarks.

Concluding Remarks

In this paper, we present a study of instruction duplication for detecting soft errors in VLIW datapaths and present compiler algorithms for this purpose. The characteristic of the algorithm is that it improves reliability under performance constraints. Overall, this paper shows that our algorithm allows the designer to conduct tradeoff analyses between performance and reliability.

Fig 2 :

 2 Fig 1: Architectural support

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

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

* This work was supported by NSF grant no. 0093082, 0093085 and 0202007


Nikolaos Kavvadias

Spiridon Nikolaidis

Hardware support for arbitrarily complex loop structures in embedded applications

In this paper, the program control unit of an embedded RISC processor is enhanced with a novel zerooverhead loop controller (ZOLC) supporting arbitrary loop structures with multiple-entry/exit nodes. The ZOLC has been incorporated to an open RISC processor core to evaluate the performance of the proposed unit for alternative configurations of the selected processor. It is proven that speed improvements of 8.4% to 48.2% are feasible for the used benchmarks.

Introduction

Last years, the embedded processor market is dominated by new 32-bit RISC architectures (ARM, MIPS32), and embedded DSPs (Motorola 56300, ST120, TMS320C54x) featuring architectural and power consumption characteristics suitable to portable platforms. Some of these new features provide better means for the execution of loops, by surpassing the significant overhead of the loop overhead instruction pattern which consists of the required instructions to initiate a new iteration of the loop.

Generally, looping cycle overheads are confronted by using branch-decrement instructions, zero-overhead loops or customized units for more complex loop nests [START_REF] Campi | IP-reusable 32-bit VLIW RISC core[END_REF], [START_REF] Talla | Bottlenecks in multimedia processing with SIMD style extensions and architectural enhancements[END_REF]. For the XiRisc processor [START_REF] Campi | IP-reusable 32-bit VLIW RISC core[END_REF], branch-decrement instructions can be configured prior synthesis. For the DSP56300 supporting up to 7 nesting levels there is a 5cycle overhead applied even to the innermost loops. Also, a single-cycle multiple-index update unit for perfect loop nests has been described in [START_REF] Talla | Bottlenecks in multimedia processing with SIMD style extensions and architectural enhancements[END_REF]. Its main advantage is that successive last iterations of nested loops are performed in a single cycle. In contrast to our approach, only perfect loop nests are supported and the area requirements grow proportionally to the considered number of loops.

In our approach, a ZOLC method is introduced that eliminates the loop overheads and can be applied to an arbitrary combination of loops. The initialization of ZOLC presents only a very small cycle overhead since it occurs outside of loop nests. The unit has been incorporated to the XiRisc 32-bit processor [START_REF] Campi | IP-reusable 32-bit VLIW RISC core[END_REF], which is distributed as a VHDL soft-core.

Incorporating ZOLC to a programmable processor

A block diagram indicating how the proposed ZOLC architecture is incorporated in the control path of a typical RISC processor is shown in Fig. 1. The purpose of ZOLC is to provide a proper candidate program counter (PC) target address to the PC decoding unit for each substituted looping operation. The instruction decoder, the PC decoding unit and the general-purpose register file communicate with the ZOLC hardware. ZOLC is composed from the task selection unit, which determines the appropriate next PC value when execution resides in a loop structure, the loop parameter tables where the loop bound values are kept and the index calculation unit.

Two modes of operation are distinguished for the ZOLC. In "initialization" mode, the ZOLC storage resources are initialized with the known loop bound values and the loop structure encoding by a special instruction sequence. In "active" mode, the ZOLC: a) determines the following task, b) issues a new target PC value and a set of candidate exit values for the case of multiple-exit loops to PC decode, c) loop indices are updated and written back to the integer register file. The task sequencing information (tasks are defined as control-flow graph regions among loop boundaries) is stored in a LUT within the task selection unit. On completion of a task, a task end signal is issued from PC decode, and an entry is selected from the LUT to address the succeeding task and the loop parameter blocks, based on which task has completed and the current loop status.

We have implemented three different configurations of a ZOLC engine. ZOLCfull refers to a ZOLC supporting 32 task switching entries, and 8-loop structure with up to 4 entries/exits/entries per loop. ZOLClite lacks support for multiple-entry/exit and uZOLC, is usable for single loops. Along with the three variations of ZOLC, two instances of the XiRisc processor family are invoked, the unmodified core noted as XRdefault, and XRhrdwil employing branch-decrement instructions.

It was found for uZOLC, ZOLClite and ZOLCfull, that the requirements in storage resources are 30, 258 and 642 storage bytes and in combinational area 298, 4056, and 4428 equivalent gates, respectively. The processor cycle time is not affected due to ZOLC and corresponds to about 170MHz on a 0.13 m ASIC process. The relative cycle measurements of Fig. 2 for comparing ZOLClite against two XiRisc configurations have been acquired for a set of 12 benchmark applications, collected from the XiRisc validation suite [START_REF] Campi | IP-reusable 32-bit VLIW RISC core[END_REF], and software implementations of motion estimation kernels. The use of branch-decrement instructions provides a cycle reduction of up to 27.5% and about 11.1% in average, while incorporating the ZOLC unit is responsible for improvements of up to 48.2% and about 26.2% in average.

Conclusion

In this paper, a zero-overhead loop controller is introduced and incorporated to the XiRisc processor. The presented architecture is able to execute structured algorithms for an arbitrary combination of loops, with no cycle overheads incurred for task switching. For a representative benchmark suite, execution time improvements up to 48% are reported against alternative XiRisc configurations.

Figure 1 .Figure 2 .

 12 Figure 1. Incorporating the ZOLC architecture to programmable RISC processors

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


Jairo Cugliari
email: jairo.cugliari@univ-lyon2.fr

José G Gómez

Assessing Imputation of Extreme Data on Climatological Time Series

Techniques for imputation of missing data are generally well known. However, they are often applied and validated on non-extreme data. Even for specific approaches specialized on extreme data, the question raises on how to assess the imputation. Because of this, we propose to use known estimators of extreme value theory, such as the extremogram, which we generalize to time series with the presence of missing data.

Introduction

Recently extreme temperature events took important relevance in global scale studies. However, there is a lack of information and/or studies on Africa and South America where the geographical inhomogeneity of the results is caused mainly due to the absence of data. A special case is Uruguay, where long time series are available, but incomplete, with missing records on several consecutive years. On the other hand, climate change, annual seasonality and exogenous factors are some of the well known non stationary patterns present in climatological records. We use a flexible non parametric model to extract these factors from a non stationary time series. The remaining noise will then be assumed as a strictly stationary R d -valued time series, denoted by X := (X i ) i∈Z .

Preliminaries : The extremogram [2]

Def. 1.The extremogram of X, for two sets A and B bounded away from zero, is defined (provided the limit exists) as

ρ X A,B (h) := lim x→∞ P   x -1 X h ∈ B|x -1 X 0 ∈ A   , h = 0, 1, . . . . ( 1
)
For the observations X 1 , . . . , X n , the sample extremogram is given by

ρ X A,B,n (h) := n-h i=1 I {u -1 m X i+h ∈B,u -1 m X i ∈A} n i=1 I {u -1 m X i ∈A} , ( 2
)
where u m is the

Our framework

The aim of this work is to study the extremogram of X when in the time series X some observations might be missed. For this, we suppose that the observations of the data (X i ) i=1,...,n occur at times

1 = i 1 < i 2 < • • • < i m = n.
From [START_REF] Parzen | On spectral analysis with missing observations and amplitude modulation[END_REF], we use the amplitude modulated observations of X, which is defined by

Y i := b i X i , i = 1, . . . , n (3)
where b i = 1 if X i is observed and b i = 0 if X i is missing. Moreover, we will assume that (X i ) i and (b i ) i are independent sequences.

Results

Proposition 1 Suppose that (b i ) i satisfy the condition (A.1) and that X is regularly varying. Then, for two sets A and B bounded away zero, ρ X A,B (h) and ρ Y A,B (h) exist for h = 0, 1, 2, . . ., and the following equality holds

ρ Y A,B (h) = ρ b (h) ρ X A,B (h), (6)
for each h = 0, 1, 2, . . ., where ρ b (h

) := P(b h = 1|b 0 = 1) = ν b (h)/µ b .
From Proposition 1, we can naturally provide an estimator of the extremogram of the sequence X = (X i ) i through its amplitude modulated version (Y i ) i and the sequence (b i ) i as follows:

ρ X A,B,n (h) = ρ Y A,B,n (h) ρ b,n (h) ,
provided that ρ b,n (h) = 0, where ρ b,n (h) := ν b,n (h)/ µ b .

Theorem 1 Suppose that (b i ) i satisfies the conditions (A.1)-(A.4) and assume that X = (X i ) i is regularly varying with index α > 0.

If condition (D) holds, α r = o(m/n) and m = o(n 1/3 ), then       n m       1/2    ρ X A,B,n (h) -ρ X A,B,n (h)    h=0,1,...,H D -→ N (0, Σ A,B ), (7)
where the asymptotic covariance is defined in [2, Corollary 3.3] and ρ X A,B,n (h After following the imputation strategy presented in [START_REF] Cugliari | Imputation of Temperature Extremes using Generalized Additive Models[END_REF], we obtain the extremogram for the original and imputed data. The graphic of this object is in Figure 2. Both lines (black and red) follow a reasonable similar pattern. Then, both series present the same extreme behavior and thus we conclude that the imputation scheme is appropriate.

) := P(u -1 m X h ∈ B|u -1 m X 0 ∈ A).

Application on real data

Figure 2: Extremogram for the original and imputed series, both by the method proposed in [START_REF] Cugliari | Imputation of Temperature Extremes using Generalized Additive Models[END_REF] and using a naive approach.

Conclusion & future work

• Use measures of extremes constructed from cluster functionals [START_REF] Gómez | Dependent Lindeberg central limit theorem for the fidis of empirical processes of cluster functionals[END_REF].

• Adaptation for usage on climatological random fields.

• Modeling the absence/presence of missing data through a 2-state Markov process.

• Assume only local stationarity on the X processes.

( 1 -

 1 1/m)-quantile of the distribution of |X 0 |. In order to have a consistent result, we require m = m n -→ ∞ with m = o(n) as n → ∞, and u m is the sequence used in the definition of regularly varying time series in [2, § 1.2].

P

  Assumptions for the process (b i ) (A.1) (b i ) has finite second moments and, for h = 0, 1, . . ., γ b (h) := Cov (b i , b i+h ) is independent of i ∈ Z. (A.2) µ b,n := 1 b := Eb 0 = P(b 0 = 1) (A.3) ν b,n (h) := 1 b (h) := Eb 0 b h , h = 0, 1, . . .. (A.4) µ b = 0 and ν b (h) = 0 for each h = 0, 1, . . .. Assumptions for the process (X i ) (D) (X i ) i∈Z is α-mixing with rate function (α l ) l∈N . Moreover, there exist m = m n and r = r n -→ ∞ with m/n -→ 0 and r/m -→ 0, such that lim n→∞ m (|X h | > a m , |X 0 | > a m ) = 0, ∀ > 0. (5)

Data.Figure 1 :

 1 Figure 1: Geographical location of measuring stations (l.) and missing values patterns (r.).


Abhishek Somani
email: somani@vlsi.iitkgp.ernet.in

P P Chakrabarti

Amit Patra

Mixing Global and Local competition in Genetic Optimization based Design Space Exploration of Analog circuits

The knowledge of optimal design space boundaries of component circuits can be extremely useful in making good subsystem-level design decisions which are aware of the parasitics and other second-order circuit-level details. However, direct application of popular Multi-objective genetic optimization algorithms were found to produce Pareto fronts with poor diversity for analog circuits problems. This work proposes a novel approach to control the diversity of solutions by paritioning the solution space, using Local Competition to promote diversity and Global competition for convergence, and by controlling the proportion of these two mechanisms by a Simulated Annealing based formulation. The algorithm was applied to extract numerical results on analog switched capacitor integrator circuits with a wide range of tight specifications. The results were found to be significantly better than traditional GA based uncontrolled optimization methods.

Introduction

Including an optimization framework in the design of analog circuits has become fairly important in recent times. Typical analog circuit sizing optimization problems require handling of multiple, usually non-commensurate and often conflicting cost criteria. Earlier efforts in this area focussed on methodologies which converted some of the design objectives into constraints and combined the remaining into a single objective function [START_REF] Alpaydin | An evolutionary approach to automatic synthesis of high-performance analog integrated circuits[END_REF]. However, comparatively newer analog systems design methodologies [START_REF] Bernardinis | A methodology for system-level analog design space exploration[END_REF] require design of component circuits with optimal design surfaces instead of an optimal design point. At the circuit level, this translates into a problem of optimization to obtain parameterized behavior of circuits. An example of such a parameter is the load capacitance for analog components in MOS integrated analog systems. This necessitates a methodology for generating multiple optimal solutions for the circuits in the entire range of the parameterized variables/specifications. Traditionally, multi-objective optimization methods try to address this problem of generating near-optimal, non-dominated design points (Pareto Front) that are well distributed over the design space [START_REF] Deb | Multi-Objective Optimization using Evolutionary Algorithms[END_REF][4] [START_REF] Smedt | Watson: Design space boundary exploration and model generation for analog and RFIC design[END_REF] [START_REF] Smedt | Holmes: Capturing the yield optimized design space boundaries of analog and RF integrated circuits[END_REF].

One method of solving a multi-objective circuit optimization problem is to transform it into a set of scalarized single objective optimization problems by the weighted sum approach or the Normal-Boundary Intersection method [START_REF] Stehr | Performance trade-off analysis of analog circuits by normal-boundary intersection[END_REF]. However, methods applying deterministic algorithms like steepest descent or quadratic programming for solving the single objective optimization problems have been known to get trapped in local optima. Recently, application of sequential stochastic refinement methods like Genetic Algorithms(GA) for analog circuits design space exploration have been successfully reported in [START_REF] Smedt | Watson: Design space boundary exploration and model generation for analog and RFIC design[END_REF] [START_REF] Smedt | Holmes: Capturing the yield optimized design space boundaries of analog and RF integrated circuits[END_REF]. The population based approach of GA makes it very suitable for multiobjective optimization problems, where multiple solutions are sought in a single run of the algorithm. This work explores the possibility of obtaining near-optimal design surfaces using GA methods for analog circuits. An example analog circuit is described in Section 2. Results obtained by application of known GA approaches to this problem are presented in Section 3. We then propose a variation based on mixture of local and global competition in GA and present the associated results in Section 4. For the sake of clarity of expression and explanation, these algorithms have been described for two-objective problems. The extension to an arbitrary number of objective functions is straightforward. Section 5 concludes the paper with a brief discussion on the results in the previous section. Our chosen example sub-circuit is a correlated double sampling (CDS) offset-compensated switched-capacitor integrator, which is the basic building block for sigma-delta modulators used extensively to build filters and data converters. We wish to use the optimal design surface of this circuit for the construction of a fourth-order sigma-delta modulator. Fig. 1 shows a typical CDS offset-compensated integrator with sampling, feedback and offset cancellation capacitors in the differential configuration.

Integrator Circuit

Standard two-stage opAmp topology was used in the integrator for this work. Apart from the gate capacitances of MOSFETs, we consider bottom-plate parasitic capacitances of standard integrated capacitors and drain diffusion and overlap capacitances of MOSFETs for more accurate description and prediction of circuit behaviour. We combine velocity saturation effect and advanced mobility models to arrive at the MOSFET model (shown in eqn.1) for deep submicron MOS circuit optimization.

I D = 1 2 µC ox W L (V GS -V T ) 2 1 -VGS -VT EsatL (1 + λV DS ) 1 + θ 1 (V GS + V T -V K ) 1 3 + θ 2 (V GS + V T -V K ) n
(1) where θ 1 , θ 2 and V K are fitting parameters, and the rest of the symbols have their usual meanings. Here, n = 1 for NMOS and 2 for PMOS. The target technology is an industry-standard 0.18µm, 1.8 Volts, n-well digital CMOS process.

After initial topology based reduction of design parameters, the optimization problem is framed with 15 design parameters. The two objectives are Power dissipation and the Load capacitance. The other circuit performance parameters, viz., Dynamic Range (DR), Settling Time (ST), Settling error (SE), Area, Output voltage range (OR) and Yield Calculation [START_REF] Smedt | Holmes: Capturing the yield optimized design space boundaries of analog and RF integrated circuits[END_REF] (Robustness) were fixed as circuit constraints. Further the problem was constrained by the requirement of all the transistors to be in the proper DC operating region and matching constraints across all manufacturing process corners. It may be noted here that the effect of non-dominant poles and zeros was included in the equations for ST, SR and DR, which makes them more non-linear than those obtained by standard dominant pole analysis based derivations. This helps in making the entire search-space available to the optimizer by means of a small number of analytical equations. 20 different specifications of the circuit graded by their level of difficulty were chosen and the Pareto Fronts were obtained for all of them using the algorithms described in the following sections. The results obtained for all of the cases showed some common trends, which we will present with explicit illustrations for only one of them due to space limitations in this paper. The specifications for the chosen case are : DR ≥ 96 dB, OR ≥ 1.4 volts, ST ≤ 0.24 µsec, SE ≤ 7 × 10 -4 , Robustness ≥ 0.85.

Applying known GA approaches

A popular Multi-objective Genetic optimization algorithm called Elitist non-dominated sorting GA (NSGA-II) [START_REF] Deb | Multi-Objective Optimization using Evolutionary Algorithms[END_REF] with a real-parameter continuous space was chosen for approaching the above problem. However, the results obtained were not very encouraging. Fig. 2 shows the Pareto-optimal front obtained after 800 iterations of NSGA-II. Here, it was desired that the algorithm produce Paretooptimal solutions which are well-distributed over the entire range of Load capacitance from 0 to 5 pF, but the solutions were found to cluster mostly between 4 and 5 pF. Investigations into such clustered behavior of the Paretooptimal solutions revealed the following :

• Loss of diversity : If in the early phase of GA, a comparatively higher number of solutions are produced in a small part of the objective function space in comparison to the rest of the space, the "cross-over" operator tends to produce more solutions "in between" the solutions in the cluster and thus tends to make the cluster more dense. • Premature Termination : Comparatively weaker solutions from the rest of the space are gradually lost due to failure to win in the global competition. This leads to further loss in diversity and ultimately premature termination of evolution. • Population Size: Clustering effect can be reduced by increasing population size considerably, but this increases the computational cost also.

Partition based Multiple Objectives Optimization Methodology

The core problem identified in the previous section was to introduce mechanisms which can preserve diversity in the population while not sacrificing the speed of convergence to the near-optimal Pareto Front in a resource-efficient manner. A first approach for increasing diversity is to increase the mutation probability. Increased mutation helps upto a certain extent beyond which the entire optimization process becomes random and loses the focus required for convergence in reasonable time. Also, the optimal value of mutation probability is heavily problem-dependent and its estimation is non-trivial.

Contribution of this work

A known method of diversity preservation is parallel population GA with inter-population migration controlled in a tribe or island based framework [START_REF] Homayounfar | An island-based GA for dynamic optimization problems[END_REF], which can be extended for Multi-objective GA. However, in this work, we try to establish that this objective can be accomplished by a simple modification in the traditional single-population GA. We argue that if the relatively weaker solutions from sparsely populated regions of the objective function space can be protected from global competition and evolved until they become strong enough to compete globally, one would get reasonably diverse solutions in the final Pareto optimal front. To achieve this, a Simulated Annealing driven Competition in Genetic Algorithm (SACGA) is proposed in this paper, where the objective function space is partitioned and a purely Local Competition in each partition transforms into a purely Global competition in the entire space in the course of the evolution of the GA. Further, a Multi-phase Expanding-partitions SACGA (MESACGA) is proposed to speed up the rate of advancement of the Pareto Front.

Performance Metric

The performance metric for the results obtained from a Multi-objective Optimization Algorithm should evaluate both the diversity of the solutions and their convergence to the near-optimal Pareto front. One such metric for minimization problems is the Hypervolume Metric [START_REF] Deb | Multi-Objective Optimization using Evolutionary Algorithms[END_REF]. Mathematically, for each solution i, a hypercube v i is constructed with the origin of the objective function space and the solution i as the diagonal corners of the hypercube. Thereafter, a union of all hypercubes is found and its hypervolume is calculated. Obviously, a lower value of the hypervolume indicates a better pareto front. This method produces reasonably diverse solutions and keeps the eligible candidates from partitions falling within the feasible objective function space alive for large number of iterations, thus allowing them improved chance to participate in the Global Mating Pool. However, after the initial stage of a couple of hundred iterations, the Global Pareto Front was found to progress towards the near-optimal Pareto front extremely slowly and would require large number iterations to reach there. The reason for this observation can be better explained by defining a few terms: The population members on the Pareto front of the global population are called Globally Superior solutions and the remaining are termed Globally Inferior solutions. Similarly, the best individuals obtained from Local competition of individuals in a subset of the population form the Local Pareto Front and are termed Locally Superior solutions. Thus, in a GA, the cross-over operation has a better chance of producing more Globally Superior solutions when operating on two Globally Superior solutions, than when operating on two Globally Inferior solutions or a Globally Superior solution with a Globally Inferior solution. In absence of any Global Competition, there are typically many Locally Superior solutions which are Globally Inferior. This causes slow production of Globally Superior solutions.

Local Competition only

Mixing Local and Global Competition

It is evident that while Local competition improves the diversity of the solutions along the Global Pareto Front, Global Competition improves the speed of movement of the front towards the global near-optimal front. We propose a Simulated Annealing driven Competition in Genetic Algorithm (SACGA) to allow the evolution process in the algorithm to gradually transit from pure Local competition in the initial phase to pure Global Competition in the final phase of the duration of solving the problem. Fig. 3 shows the flow of this algorithm. The phase of mixed local and global competition is preceded by a phase of pure local competition which ends after each partition has evolved at least one solution which satisfies all the constraints. However, some partitions may be in the infeasible region. So, we define an upper limit on the number of iterations allowed in this phase, after which partitions with no constraintsatisfying solutions are discarded. Let n be the desired number of globally non-dominated (i.e., superior) solutions per partition at the end of the second phase and gen t be the number of iterations required for the first phase. Let span be the total number of iterations allowed in the second phase.

In any iteration in the second phase, let m p be the number of locally superior solutions in partition p and c be the cost of allowing a locally superior solution to participate in the global competition. Then, we formulate :

c = k 1 e k2 i n-1 (2)
where, k 1 ,k 2 are constants to be chosen for shaping the cost functions and i assumes the values 1, 2, . . . , m p for partition p. Thus, the locally superior solution i will participate in

prob = 1 -e -α c×T A (3)
where α is a chosen constant and Annealing temperature T A is given by the temperature schedule :

T A = T init e -k3× ln(T init ) span ×(gen-gent) (4)
where, gen is the current iteration number and T A is cooled down from T init to a value of 1. This formulation of probability prob provides the following features :

1. It encourages pure local competition in the initial iterations and pure global competion in the final iterations of this phase. 2. The more the number of locally superior solutions produced by a partition, the better the partition is. However, we do not allow all such "good" solutions to participate in the Global Competition till the very final phase of the evolution. This is necessary because if some of these solutions are retained, then the partition will maintain its representation in the Global pool even if all its participants in the Global Competition are dominated by solutions from other partitions. This is accomplished by assigning different probabilities to different locally superior solutions. Thus, if the m p "good" solutions in a partition p are considered in a randomly determined sequence for evaluation of their candidature for global competition, the solutions considered earlier have a higher probability than those considered later. 3. The shapes of the probability curves can be easily controlled by selecting the parameters k 1 ,k 2 and k 3 for desired values of probability at iteration gen = gen t + span 2 for i = 1, n and that at iteration gen = gen t + span for i = 1, . . . , n. Fig. 4 shows the probability curves for n = 5 and span = 100. In the example circuit problem, the objective function space partitioning is induced by the division of the range space of the Load Capacitance. Fig. 5 presents the Pareto Front obtained after 800 iterations of an 8-partition SACGA and its comparison with that obtained earlier by NSGA-II using purely global competition. A prominent issue which affects the efficiency of SACGA is the problem of selecting the optimal number of partitions with respect to each objective function and determining their (generally, and no method is known of finding them. A simplified approach may be to choose partitions of equal sizes. However, it is intuitive and our studies have shown that the optimality of the solutions obtained by a finite number of iterations depends upon the chosen number of partitions. Fig. 6 shows the quality of results obtained after 1200 iterations for the example problem. Dividing the objective function space into 16 equal partitions was found to be optimal for this problem instance. At present, no alternative to complete experimentation is known for finding this number. It is to circumvent this problem coupled with a few other advantages that we suggest a modified version of SACGA in the next subsection.

Multi-phase Expanding-partitions SACGA (MESACGA)

A better option to chosing a set of static partitions would be to start with a high number of partitions and gradually expand the partition-size, thereby reducing their number. Such a method would essentially "grow" the individual Local Pareto Fronts till they combine together to form the Global Pareto Front in the final phase having just one partition covering the entire objective function space. Essentially, MESACGA is a SACGA running in multiple phases where the number of partitions is reduced and their size increased at the end of each phase. A brief overview is shown in Fig. 7. MESACGA provides a systematic method for striking a balance between preserving Pareto front diversity and convergence to global near-optimal. In MESACGA, at the end of each phase of SACGA, the partitions start to cover a wider region in the objective function space in the current phase as compared to that in the previous phase. This causes some Locally Superior but Globally Inferior solutions to be thrown away and there is a tendency of decrease in the proportion of Globally Inferior solutions in the pool of Locally Superior solutions and hence, faster movement of the pareto front. However, the presence of such solutions over periods of iterations is necessary for preserv-ing diversity. The parameter span, i.e., the number of iterations for each phase of MESACGA provides a direct control to maintain the above mentioned balance. Thus, a smaller value of span would give quick results which may not cover a significant portion of the objective function space, whereas, a higher value of span would make the conversion of Local Pareto Fronts to the Global Pareto Front more gradual and would, in general, produce solutions with higher diversity. In addition, MESACGA was found to produce results comparable in quality to those produced by selecting the best possible number of equal partitions in SACGA.

Discussions and Conclusion

The results obtained for all the 20 different circuit specifications showed similar trends with both SACGA and MESACGA, taking, on an average, 18% more computa-tional time compared to NSGA-II, due to additional overheads of these algorithms. In brief, the trends were as follows :

1. In all cases, where the evolution was continued for more than 650 iterations, the quality of the solutions obtained for equal number of iterations for TPG (Traditional Purely Global competition based GA), SACGA and MESACGA were found to be in the order MESACGA ≥ SACGA ≥ TPG. Fig. 8 shows the Pareto Fronts obtained for 800 iterations from the three approaches. However, no clear trends were visible for lesser number of iterations. 2. It was observed in SACGA, that the differential gain obtained by increasing the number of iterations beyond 700 turned out to be very low.

Figure 1 .

 1 Figure 1. CDS Offset-compensated Switched-Capacitor Integrator

Figure 2 .

 2 Figure 2. Pareto-optimal Front after 800 iterations of NSGA-II

  Given two objective functions, we propose to partition the associated function space into m equal partitions induced by the division of the range space of any one of the objective functions into m equal, disjoint partitions. In each iteration, only Local competition is allowed within each partition for non-dominated ranking of individuals. A Global Mating Pool is created by rank-based selection of individuals from the entire population of solutions at the end of each iteration and the candidates for next iteration are created by performing cross-over and mutation operations in the Global Mating Pool. After a fixed number of iterations, Global Competition is performed once on the entire population to find the Global Pareto Front containing the set of globally non-dominated solution points.

Figure 3 .

 3 Figure 3. SACGA Algorithm flow global competition with probability prob given by :

Figure 4 .Figure 5 .

 45 Figure 4. Probability curves for n=5

Figure 6 .

 6 Figure6. Determination of Optimal number of partitions unequal) sizes. They are dependent upon the solution space and no method is known of finding them. A simplified approach may be to choose partitions of equal sizes. However, it is intuitive and our studies have shown that the optimality of the solutions obtained by a finite number of iterations depends upon the chosen number of partitions. Fig.6shows the quality of results obtained after 1200 iterations for the example problem. Dividing the objective function space into 16 equal partitions was found to be optimal for this problem instance. At present, no alternative to complete experimentation is known for finding this number. It is to circumvent this problem coupled with a few other advantages that we suggest a modified version of SACGA in the next subsection.

Figure 7 .Figure 8 .

 78 Figure 7. MESACGA Algorithm FlowMESACGA was applied to the example problem with 7 successive phases of SACGA containing 20, 13, 8, 5, 3, 2, 1 partitions respectively. Fig.8shows the MESACGA results in comparison to the both traditional and SACGA approaches.

Fig. 9

 9 plots the Hypervolume metric calculated for the Pareto fronts obtained by setting the span parameter to progressively increasing values in an 8-partition SACGA. It is clear that not much improvement of the Pareto front is obtained for span > 1000. 3. Fig.10 shows the Hypervolume metric calculated for the Pareto fronts at the end of each phase of SACGA in a 7-phase MESACGA. The value of span was kept the same for each phase and the results were found to improve by increasing the value of span. It was shown in Fig.6 that a 1200 iteration SACGA produces best results for 16 partitions. It can be clearly observed in Fig.11 that a 1250 iteration MESACGA having a purely local competition phase of 200 iterations followed by 7 phases of 150 iterations each produces comparable results. The Hypervolume metric calculated for the Pareto fronts are 22.19 and 21.83 for SACGA and MESACGA respectively.

Figure 9 .

 9 Figure 9. Performance of SACGA for various preset values of total number of iterations Thus, we conclude that the poor diversity of the Pareto fronts obtained from application of traditional, purely global competition based Multi-objective genetic optimization algorithms can be remedied to a large extent by the proposed MESACGA approach.

Figure 11 .

 11 Figure 10. Progress of Pareto Front across various SACGA phases of MESACGA

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)


Dirk Husmeier

Frank Dondelinger

Sophie Lèbre
email: sophie.lebre@lsiit-cnrs.unistra.fr

Inter-time segment information sharing for non-homogeneous dynamic Bayesian networks

Conventional dynamic Bayesian networks (DBNs) are based on the homogeneous Markov assumption, which is too restrictive in many practical applications. Various approaches to relax the homogeneity assumption have recently been proposed, allowing the network structure to change with time. However, unless time series are very long, this flexibility leads to the risk of overfitting and inflated inference uncertainty. In the present paper we investigate three regularization schemes based on inter-segment information sharing, choosing different prior distributions and different coupling schemes between nodes. We apply our method to gene expression time series obtained during the Drosophila life cycle, and compare the predicted segmentation with other state-of-the-art techniques. We conclude our evaluation with an application to synthetic biology, where the objective is to predict a known in vivo regulatory network of five genes in yeast.

Introduction

There is currently considerable interest in structure learning of dynamic Bayesian networks (DBNs), with a variety of applications in signal processing and computational biology; see e.g. [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF][START_REF] Grzegorczyk | Non-stationary continuous dynamic Bayesian networks[END_REF][START_REF] Ahmed | Recovering time-varying networks of dependencies in social and biological studies[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. While there have been various efforts to relax the homogeneity assumption for undirected graphical models [START_REF] Talih | Structural learning with time-varying components: Tracking the cross-section of financial time series[END_REF][START_REF] Xuan | Modeling changing dependency structure in multivariate time series[END_REF], relaxing this restriction in DBNs is a more recent research topic [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF][START_REF] Grzegorczyk | Non-stationary continuous dynamic Bayesian networks[END_REF][START_REF] Ahmed | Recovering time-varying networks of dependencies in social and biological studies[END_REF][START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF][START_REF] Kolar | Sparsistent learning of varying-coefficient models with structural changes[END_REF]. At present, none of the proposed methods is without its limitations, leaving room for further methodological innovation. The method proposed in [START_REF] Ahmed | Recovering time-varying networks of dependencies in social and biological studies[END_REF][START_REF] Kolar | Sparsistent learning of varying-coefficient models with structural changes[END_REF] is non-Bayesian. This requires certain regularization parameters to be optimized "externally", by applying information criteria (like AIC or BIC), cross-validation or bootstrapping. The first approach is suboptimal, the latter approaches are computationally expensive 1 . In the present paper we therefore follow the Bayesian paradigm, like [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF][START_REF] Grzegorczyk | Non-stationary continuous dynamic Bayesian networks[END_REF][START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF]. These approaches also have their limitations. The method proposed in [START_REF] Grzegorczyk | Non-stationary continuous dynamic Bayesian networks[END_REF] assumes a fixed network structure and only allows the interaction parameters to vary with time. This assumption is too rigid when looking at processes where changes in the overall regulatory network structure are expected, e.g. in morphogenesis or embryogenesis. The method proposed in [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF] requires a discretization of the data, which incurs an inevitable information loss. These limitations are addressed in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF], where the authors propose a method for continuous data that allows network structures associated with different nodes to change with time in different ways. However, this high flexibility causes potential problems when applied to time series with a low number of measurements, as typically available from systems biology, leading to overfitting or inflated inference uncertainty. The objective of the work described in our paper is to propose a model that addresses the principled shortcomings of the three Bayesian methods mentioned above. Unlike [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF], our model is continuous and therefore avoids the information loss inherent in a discretization of the data. Unlike [START_REF] Grzegorczyk | Non-stationary 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][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF], our model introduces information sharing among time series segments, which provides an essential regularization effect.

Background: non-homogeneous DBNs without information coupling

This section summarizes briefly the non-homogeneous DBN proposed in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF], which combines the Bayesian regression model of [START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] with multiple changepoint processes and pursues Bayesian inference with reversible jump Markov chain Monte Carlo (RJMCMC) [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF]. 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 p be the number of observed genes, whose expression values y = {yi(t)} 1≤i≤p,1≤t≤N are measured at N time points. M represents a directed graph, i.e. the network defined by a set of directed edges among the p genes. M i is the subnetwork associated with target gene i, determined by the set of its parents (nodes with a directed edge feeding into gene i). The regulatory relationships among the genes, defined by M, may vary across time, which we model with a multiple changepoint process. For each target gene i, an unknown number k i of changepoints define k i + 1 non-overlapping segments. Segment h = 1, .., k i + 1 starts at changepoint ξ h-1 i and stops before ξ h i , where ξi = (ξ

0 i , ..., ξ h-1 i , ξ h i , ..., ξ k i +1 i ) with ξ h-1 i < ξ h i . To delimit the bounds, ξ 0 i = 2 and ξ k i +1 i = N + 1. Thus vector ξ i has length |ξi| = ki + 2.
The set of changepoints is denoted by ξ = {ξi} 1≤i≤p . This changepoint process induces a partition of the time series,

y h i = (yi(t)) ξ h-1 i ≤t<ξ h i , with different structures M h i associated with the different segments h ∈ {1, . . . , ki + 1}.
Identifiability is satisfied by ordering the changepoints based on their position in the time series.

Regression model: For all genes i, the random variable Y i (t) refers to the expression of gene i at time t. Within any segment h, the expression of gene i depends on the p gene expression values measured at the previous time point through a regression model defined by (a) a set of s h i parents denoted by

M h i = {j1, ..., j s h i } ⊆ {1, . . . , p}, |M h i | = s h i , and (b) a set of parameters ((a h ij )j∈0..p, σ h i ); a h ij ∈ R, σ h i > 0. For all j = 0, a h ij = 0 if j / ∈ M h i .
For all genes i, for all time points t in segment h (ξ h-1 i ≤ t < ξ h i ), the random variable Yi(t) depends on the p variables {Yj(t -1)} 1≤j≤p according to

Yi(t) = a h i0 + X j∈M h i a h ij Yj(t -1) + εi(t) (1)
where the noise εi(t) is assumed to be Gaussian with mean 0 and variance

(σ h i ) 2 , εi(t) ∼ N (0, (σ h i ) 2 ).
We define a h i = (a h ij )j∈0..p.

Prior

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

P (ki|λ) ∝ λ k i k i ! 1l {k i ≤k} . Conditional on k i changepoints, the changepoint positions vector ξi = (ξ 0 i , ξ 1 i , ..., ξ k i +1 i
) takes non-overlapping integer values, which we take to be uniformly distributed a priori. There are (N -2) possible positions for the k i changepoints, thus vector ξi has prior density

P (ξi|ki) = 1/ " N -2 k i "
. For all genes i and all segments h, the number s h i of parents for node i follows a truncated Poisson distribution2 with mean Λ and maximum s = 5:

P (s h i |Λ) ∝ Λ s h i s h i ! 1l {s h i ≤s} . Conditional on s h
i , the prior for the parent set M h i is a uniform distribution over all parent sets with cardinality s h i :

P (M h i ˛|M h i | = s h i ) = 1/( p s h i ).
The overall prior on the network structures is given by marginalization:

P (M h i |Λ) = X s s h i =1 P (M h i |s h i )P (s h i |Λ) (2)
Conditional on the parent set M h i of size s h i , the s h i + 1 regression coefficients, denoted by

a M h i = (a h i0 , (a h ij ) j∈M h i
), are assumed zero-mean multivariate Gaussian with covariance matrix

(σ h i ) 2 Σ M h i , P (a h i |M h i , σ h i )=|2π(σ h i ) 2 Σ M h i | -1 2 exp 0 @ - a † M h i Σ -1 M h i a M h i 2(σ h i ) 2 1 A (3)
where the symbol † denotes matrix transposition,

Σ M h i = δ -2 D † M h i (y)D M h i (y) and D M h i (y) is the (ξ h i -ξ h-1 i ) × (s h i + 1)
matrix whose first column is a vector of 1 (for the constant in model ( 1)) and each (j + 1) th column contains the observed values (yj(t)) ξ h-1 i -1≤t<ξ h i -1 for all factor gene j in M h i . This prior was also used in [START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] and is motivated in [START_REF] Zellner | On assessing prior distributions and Bayesian regression analysis with g-prior distributions[END_REF]. Finally, the conjugate prior for the variance (σ h i ) 2 is the inverse gamma distribution, P ((σ h i ) 2 ) = IG(υ0, γ0). Following [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[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: P (Λ) = P (λ) = Ga(0.5, 1) and P (δ 2 ) = IG(2, 0.2).

Posterior

Equation (1) implies that

P (y h i |ξ h-1 i , ξ h i , M h i , a h i , σ h i ) = " √ 2πσ h i " -(ξ h i -ξ h-1 i ) exp 0 @ - (y h i -D M h i (y)a M h i ) † (y h i -D M h i (y)a M h i ) 2(σ h i ) 2 1 A (4)
From Bayes theorem, the posterior is given by the following equation, where all prior distributions have been defined above:

P (k, ξ, M, a, σ, λ, Λ, δ 2 |y) ∝ P (δ 2 )P (λ)P (Λ) p Y i=1 P (ki|λ)P (ξi|ki) k i Y h=1 P (M h i |Λ) (5)
P ([σ h i ] 2 )P (a h i |M h i , [σ h i ] 2 , δ 2 )P (y h i |ξ h-1 i , ξ h i , M h i , a h i , [σ h i ] 2 )

Inference

An attractive feature of the chosen model is that the marginalization over the parameters a and σ in the posterior distribution of ( 5) is analytically tractable:

P (k,ξ,M,λ,Λ,δ 2 |y) = Z P (k,ξ,M,a,σ,λ,Λ,δ 2 |y)dadσ (6)
See [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] for details and an explicit expression. The number of changepoints and their location, k, ξ, the network structure M and the hyperparameters λ, Λ, δ 2 can be sampled from the posterior P (k, ξ, M, λ, Λ, δ 2 |y) with RJMCMC [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF]. A detailed description can be found in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF].

Model improvement: information coupling between segments

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 over-flexibility. 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 three methods of information sharing among time series segments, as illustrated in Figure 1. The first method is based on hard information coupling between the nodes, using the exponential distribution proposed in [START_REF] Werhli | Gene regulatory network reconstruction by Bayesian integration of prior knowledge and/or different experimental conditions[END_REF]. The second scheme is also based on hard information coupling, but uses a binomial distribution with conjugate Beta prior. The third scheme is based on the same distributional assumptions as the second scheme, but replaces the hard by a soft information coupling scheme.

Hard information coupling based on an exponential prior

Denote by Ki := ki + 1 the total number of partitions in the time series associated with node i, and recall that each time series segment y h i is associated with a separate subnetwork

M h i , 1 ≤ h ≤ Ki.
We impose a prior distribution

P (M h i |M h-1 i
, β) on the structures, and the joint probability distribution factorizes according to a Markovian dependence:

P (y 1 i , . . . , y K i i , M 1 i , . . . , M K i i , β) = K i Y h=1 P (y h i |M h i )P (M h i |M h-1 i , β)P (β) (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

P (M h i |M h-1 i , β) = exp(-β|M h i -M h-1 i |) Zi(β, M h-1 i ) (8)
for h ≥ 2, where β is a hyperparameter that defines the strength of the coupling between M h i and M h-1 i , and |.| denotes the Hamming distance. For h = 1, P (M h i ) is given by (2). The denominator Z(β, M h-1 i ) in ( 8) is a normalizing constant, also known as the partition function:

Z(β) = P M h i ∈M e -β|M h i -M h-1 i
| where M is the set of all valid subnetwork structures. If we ignore any fan-in restriction that might have been imposed a priori (via s), then the expression for the partition function can be simplified:

Z(β) ≈ Q p j=1 Zj(β), where Zj(β) = P 1 e h j =0 e -β|e h j -e h-1 j | = 1 + e -β
and hence Z(β) = `1 + e -β ´p. Inserting this expression into (8) gives:

P (M h i |M h-1 i , β) = exp(-β|M h i -M h-1 i |) (1 + e -β ) p (9)
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][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF] as described in Section 2.4. When proposing a new network structure M h i → Mh i for segment h, the prior probability ratio has to be replaced by:

P (M h+1 i | Mh i ,β)P ( Mh i |M h-1 i ,β) P (M h+1 i |M h i ,β)P (M h i |M h-1 i ,β
) . An additional MCMC step is introduced for sampling the hyperparameter β from the posterior distribution. For a proposal move β → β with symmetric proposal probability Q( β|β) = Q(β| β) we get the following acceptance probability:

A( β|β) = min  P ( β) P (β) Q p i=1 Q K i h=2 exp(-β|M h i -M h-1 i |) exp(-β|M h i -M h-1 i |) (1+e -β ) p (1+e -β ) p , 1 ff
where in our study the hyperprior P (β) was chosen as the uniform distribution on the interval [0, 10].

Hard information coupling based on a binomial prior

An alternative way of information sharing among segments and nodes is by using a binomial prior:

P (M h i |M h-1 i , a, b) = a N 1 1 [h,i] (1 -a) N 0 1 [h,i] b N 0 0 [h,i] (1 -b) N 1 0 [h,i] (10)
where we have defined the following sufficient statistics:

N 1 1 [h, i] is the number of edges in M h-1 i that are matched by an edge in M h i , N 0 1 [h, i] is the number of edges in M h-1 i for which there is no edge in M h i , N 1 0 [h, i]
is the number of edges in M h i for which there is no edge in M h-1 i , and

N 0 0 [h, i]
is the number of coinciding non-edges in M h-1 i and M h i . Since the hyperparameters are shared, the joint distribution can be expressed as:

P ({M h i }|a, b) = p Y i=1 P (M 1 i ) K i Y h=1 P (M h i |M h-1 i , a, b) = a N 1 1 (1 -a) N 0 1 b N 0 0 (1 -b) N 1 0 p Y i=1 P (M 1 i ) (11)
where we have defined

N l k = P p i=1 P K i h=2 N l k [h, i],
and the right-hand side follows from Eq. ( 10). The conjugate prior for the hyperparameters a, b is a beta distribution, P (a, b|α, α, γ, γ) 1) , which allows the hyperparameters to be integrated out in closed form:

∝ a (α-1) (1 - a) (α-1) b (γ-1) (1 -b) (γ-
P ({M h i }|α, α, γ, γ) = P ({M h i }|a, b)P (a, b|α, α, γ, γ)dadb (12)
∝ Γ(α + α) Γ(α)Γ(α) Γ(N 1 1 + α)Γ(N 0 1 + α) Γ(N 1 1 + α + N 0 1 + α) Γ(γ + γ) Γ(γ)Γ(γ) Γ(N 0 0 + γ)Γ(N 1 0 + γ) Γ(N 0 0 + γ + N 1 0 + γ)
The level-2 hyperparameters α, α, γ, γ are given a uniform hyperprior over [0, 20]. The MCMC scheme of Section 2.4 has to be modified as follows. When proposing a new network structure for node i and segment h, M h i → Mh i , the structures M h i and Mh i enter the prior probability ratio via the expression P ({M h i }|α, α, γ, γ), as

P ({M 1 i ,..., Mh i ,...,M K i i } p i=1 |α,α,γ,γ) P ({M 1 i ,...,M h i ,...,M K i i } p i=1 |α,α,γ,γ)
.

Note that as a consequence of integrating out the hyperparameters, all network structures become interdependent, and information about the structures is contained in the sufficient statistics N 1 1 , N 0 1 , N 1 0 , N 0 0 . A new proposal move for the level-2 hyperparameters is added to the existing RJMCMC scheme of Section 2.4. New values for the level-2 hyperparameters x ∈ {α, α, γ, γ} are proposed from a uniform distribution over a fixed interval. For a move x → x, the acceptance probability is:

A(x|x) = min  P ({M 1 i ,...,M K i i } p i=1 |x,{α,α,γ,γ}\x) P ({M 1 i ,...,M K i i } p i=1 |x,{α,α,γ,γ}\x) , 1 ff
where {α, α, γ, γ} \ x corresponds to {α, γ, γ} if

x designates hyperparameter α, and similarly for α, γ, γ.

Soft information coupling based on a binomial prior

We can relax the information sharing scheme from a hard to a soft coupling by introducing node-specific hyperparameters a i , b i that are softly coupled via a common level-2 hyperprior, P (ai, bi|α, α, γ, γ) ∝ a 1) , as illustrated in Figure 1(b):

(α-1) i (1 -ai) (α-1) b (γ-1) i (1 -bi) (γ-
P (M h i |M h-1 i , a i , b i ) = (a i ) N 1 1 [h,i] (1 -a i ) N 0 1 [h,i] (b i ) N 0 0 [h,i] (1 -b i ) N 1 0 [h,i] (13)
This leads to a straightforward modification of eq. ( 11) -replacing a, b by a i , b i -from which we get as an equivalent to (13), using the definition

N l k [i] = Ki h=2 N l k [h, i]: P (M 1 i , . . . , M K i i |α, α, γ, γ) ∝ Γ(α + α) Γ(α)Γ(α) Γ(N 1 1 [i] + α)Γ(N 0 1 [i] + α) Γ(N 1 1 [i] + α + N 0 1 [i] + α) Γ(γ + γ) Γ(γ)Γ(γ) Γ(N 0 0 [i] + γ)Γ(N 1 0 [i] + γ) Γ(N 0 0 [i] + γ + N 1 0 [i] + γ) (14)
As in Section 3.2, we extend the RJMCMC scheme from Section 2.4 so that when proposing a new network structure, M h i → Mh i , the acceptance probability has to be updated with the prior ratio:

P (M 1 i ,..., Mh i ,...,M K i i |α,α,γ,γ) P (M 1 i ,...,M h i ,...,M K i i |α,α,γ,γ)
. In addition, we have to add a new level-2 hyperparameter update move x → x, where the prior and proposal probabilities are the same as in Section 3.2, and the acceptance probability becomes:

A(x|x) = min  Q p i=1 P (M 1 i ,...,M K i i |x,{α,α,γ,γ}\x) P (M 1 i ,...,M K i i |x,{α,α,γ,γ}\x) , 1
ff .

Results

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][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF]. 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 between-chain 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. When the true network is known, this allows us to construct the Receiver Operating Characteristic (ROC) curve (plotting the sensitivity or recall against the complementary specificity) and the precisionrecall (PR) curve (plotting the precision against the recall), and to assess the network reconstruction accuracy in terms of the areas under these graphs (AUROC and AUPRC, respectively); see [START_REF] Davis | The relationship between precision-recall and ROC curves[END_REF].

Comparative evaluation on simulated data

We randomly generated 10 networks with 10 nodes each, with the number of parents per node drawn from a Poisson distribution with mean λ = 3. To simulate changes in the network structure, we created 4 different network segments by drawing the number of changes from a Poisson distribution and applying the changes uniformly at random to edges and non-edges in the previous segment. For each segment, we generated a time series of length 15 using a linear regression model. The regression weights were drawn from a Gaussian N (0, 1), and Gaussian observation noise N (0, 1) was added. We compared the network reconstruction accuracy of the non-homogeneous DBN without information sharing proposed in [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF] (HetDBN-0) with the three information sharing approaches, based on the exponential prior from Section 3.1 (HetDBN-Exp), the binomial prior with hard node coupling from Section 3.2 (HetDBN-Bino1), and the binomial prior with soft node coupling from Section 3.3 (HetDBN-Bino2). Figures 2(a) and 2(b) shows the network reconstruction performance of the different information sharing methods in terms of AUROC and AUPRC scores. All information sharing methods show a clear improvement in network reconstruction over HetDBN-0, as confirmed by paired t-tests (p < 0.01). We investigated two different situations, the case where all segment structures are the same (although edge weights are allowed to vary) and the case where changes are applied sequentially to the segments3 . Information sharing is most beneficial for the first case, but even when we introduce changes we still see an increase in the network reconstruction scores compared to HetDBN-0. When all segments are the same, HetDBN-Bino1 and HetDBN-Bino2 outperform HetDBN-Exp (p < 0.05), but there is no significant difference between the two binomial methods. Paired t-tests showed that all other differences in mean are significant. When the segments are different, all information sharing methods outperform HetDBN-0 (p < 0.05), but the difference between the information sharing methods is not significant.

Morphogenesis in Drosophila melanogaster

We applied our methods to a gene expression time series for eleven genes involved in the muscle development of Drosophila melanogaster [START_REF] Arbeitman | Gene expression during the life cycle of Drosophila melanogaster[END_REF]. The microarray data measured gene expression levels during all four major stages of morphogenesis: embryo, larva, pupa and adult. We investigated whether our methods were able to infer the correct changepoints corresponding to the known transitions between stages. Figure 3(a) shows the posterior probabilities of inferred changepoints for any gene using HetDBN-0, while Figure 3(c) shows the posterior probabilities for the information shar- ing methods. For comparison, we applied the method proposed in [START_REF] Ahmed | Recovering time-varying networks of dependencies in social and biological studies[END_REF], using the authors' software package TESLA (Figure 3(b)). Robinson and Hartemink applied the discrete non-homogeneous DBN in [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF] to the same data set, and a plot corresponding to Figure 3(b) can be found in their paper.

Our non-homogeneous DBN methods are generally more successful than TESLA, in that they recover changepoints for all three transitions (embryo → larva, larva → pupa, and pupa → adult). Figure 3(b) indicates that the last transition, pupa → adult, is less clearly detected with TESLA, and it is completely missing in [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF]. Both our method as well as TESLA detect additional transitions during the embryo stage, which are missing in [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF]. We would argue that a complex gene regulatory network is unlikely to transition 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. We note that the main effect of information sharing is to reduce the size of the smaller peaks, while keeping the three most salient peaks (corresponding to larva → pupa, and pupa → adult, and an extra transition in the embryo phase). This reflects the fact that these changepoints are associated with significant changes in network structure, and adds to the interpretability of the results. The drawback is that the third morphological transition (embryo → larva) is less pronounced.

Reconstruction of a synthetic gene regulatory network in Saccharomyces cerevisiae

The highly topical field of synthetic biology enables biologists to design known gene regulatory networks in living cells. In the work described in [START_REF] Cantone | A synthetic network for in vivo assessment of reverse-engineering and modeling approaches[END_REF], a synthetic regulatory network of 5 genes was constructed in Saccharomyces cerevisiae (yeast), and gene expression time series were measured with RT-PCR for 16 and 21 time points under two experimental conditions, related to the carbon source: galactose ("switch on") and glucose ("switch off"). The authors tried to reconstruct the known gold-standard network from these time series with two established state-of-the-art methods from computational systems biology, one based on ordinary differential equations (ODEs), called TSNI, the other based on conventional DBNs, called Banjo; see [START_REF] Cantone | A synthetic network for in vivo assessment of reverse-engineering and modeling approaches[END_REF] for details. Both methods are optimization-based and output a single network. By comparison with the known gold standard, the authors obtained the precision (proportion of predicted interactions that are correct) and recall (proportion of predicted true interactions) scores. In our study, we merged the time series from the two experimental conditions under exclusion of the boundary point 4 , and applied the four nonhomogeneous DBNs described before. Figure 3(d) shows the inferred marginal posterior probability of potential changepoints. The most significant changepoint is at the boundary between "switch on" and "switch off" data, confirming that the known true changepoint is consistently identified.

The biological mechanism behind the other peaks is not known, and they are potentially spurious. Interestingly, the application of the proposed information-coupling schemes reduces the height of these peaks, with the binomial models having a stronger effect than the exponential one.

As we pursue a Bayesian inference scheme, we also obtain a ranking of the potential gene interactions in terms of their marginal posterior probabilities. From this we computed the precision-recall curves [START_REF] Davis | The relationship between precision-recall and ROC curves[END_REF] shown in Figure 4. Our non-homogeneous DBNs with information sharing outperform Banjo and TSNI both in the "switch on" and the "switch off" phase. They also perform better than HetDBN-0 on the "switch off" data, but are slightly worse on the "switch on" data. Note that the reconstruction accuracy on the "switch off" data is generally poorer than on the "switch on" data [START_REF] Cantone | A synthetic network for in vivo assessment of reverse-engineering and modeling approaches[END_REF]. Our results are thus plausible, suggesting that information sharing boosts the reconstruction accuracy on the poorer time series segment at the cost of a degraded performance on the stronger one. This effect is more pronounced for the exponential prior than for the binomial one, indicating a tighter coupling. The average areas under the PR curves, averaged over both phases ("switch on and off"), are as follows. HetDBN-0= 0.70, HetDBN-Exp= 0.77, HetDBN-Bino1= 0.75, HetDBN-Bino2= 0.75. Hence, the overall effect of information sharing is a performance improvement.

Conclusions

We have described a non-homogeneous DBN, which has various advantages over existing schemes: it does not require the data to be discretized (as opposed to [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF]); it allows the network structure to change with time (as opposed to [START_REF] Grzegorczyk | Non-stationary continuous dynamic Bayesian networks[END_REF]); it includes three different regularization schemes 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][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[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 time-varying networks of dependencies in social and biological studies[END_REF]). An evaluation on simulated data has demonstrated an improved performance over [START_REF] Lèbre | Stochastic process analysis for Genomics and Dynamic Bayesian Networks inference[END_REF][START_REF] Lèbre | Statistical inference of the time-varying structure of gene-regulation networks[END_REF] when information sharing is introduced. 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 [START_REF] Robinson | Non-stationary dynamic Bayesian networks[END_REF] and [START_REF] Ahmed | Recovering time-varying networks of dependencies in social and biological studies[END_REF]. We have carried out a comparative evaluation of different information coupling schemes: a binomial versus an exponential prior, and hard versus soft coupling. In an application to data from a topical study in synthetic biology, our methods have outperformed two established network reconstruction methods from computational systems biology.

Figure 1 :

 1 Figure 1: Hierarchical Bayesian models for inter-segment and inter-node information coupling. 1(a): Hard coupling between nodes with common hyperparameter Θ regulating the strength of the coupling between structures associated with adjacent segments, M h i and M h+1 i . This corresponds to the models in Section 3.1, with Θ = β, Ψ = [0, 10], and no Ω, and Section 3.2, with Θ = {a, b}, Ψ = {α, α, γ, γ}, and Ω = [0, 20]. 1(b): Soft coupling between nodes, with node-specific hyperparameters Θi coupled via level2-hyperparameters Ψ. This corresponds to the model in Section 3.3, with Θi = {ai, bi}, Ψ = {α, α, γ, γ}, and Ω = [0, 20].

Figure 2 :

 2 Figure2: Network reconstruction performance comparison of AUROC and AUPRC reconstruction scores for the four methods, HetDBN-0 (white), HetDBN-Exp (light grey), HetDBN-Bino1 (dark grey, left), HetDBN-Bino2 (dark grey, right). The boxplots show the distributions of the scores for 10 datasets with 4 network segments each, 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. "Same Segs" means that all segments in a dataset have the same structure, while "Different Segs" indicates that structure changes are applied to the segments sequentially.

Figure 3 :

 3 Figure 3: Changepoints inferred on gene expression data related to morphogenesis in Drosophila melanogaster, and synthetic biology in Saccharomyces cerevisiae (yeast). All figures using HetDBN plot the posterior probability of a changepoint occurring for any node at a given time plotted against time. 3(a): HetDBN-0 changepoints for Drosophila (no information sharing) 3(b): TESLA, L1norm of the difference of the regression parameter vectors associated with two adjacent time points plotted against time. 3(c): HetDBN changepoints for Drosophila with information sharing; the method is indicated by the legend. 3(d) HetDBN changepoints for the synthetic gene regulatory network in yeast. In 3(a)-3(c), the vertical dotted lines indicate the three morphogenic transitions, while in 3(d) the line indicates the boundary between "switch on" and "switch off" data.

  Reconstruction of a known gene regulatory network from synthetic biology in yeast. The network was reconstructed from two gene expression time series obtained with RT-PCR in two experimental conditions, reflecting the switch in the carbon source from galactose ("switch on") to glucose ("switch off"). The reconstruction accuracy of the methods proposed in Section 3, where the legend is explained, is shown in terms of precision (vertical axis) -recall (horizontal axis) curves. Results were averaged over 10 independent MCMC simulations. For comparison, fixed precision/recall scores are shown for two state-of-the-art methods reported in[START_REF] Cantone | A synthetic network for in vivo assessment of reverse-engineering and modeling approaches[END_REF]: Banjo, a conventional DBN, and TSNI, a method based on ODEs.

			Precision-Recall for Switch On				Precision-Recall for Switch Off
		1.0						1.0
		0.8						0.8
				TSNI					Banjo
									and
	Precision	0.4 0.6	Banjo				Precision	0.4 0.6	TSNI
		0.2	HetDBN-0 HetDBN-Exp					0.2	HetDBN-0 HetDBN-Exp
			HetDBN-Bino1						HetDBN-Bino1
		0.0	HetDBN-Bino2					0.0	HetDBN-Bino2
		0.0	0.2	0.4	0.6	0.8	1.0	0.0	0.2	0.4	0.6	0.8	1.0
				Recall						Recall
	Figure 4:

See[START_REF] Larget | Markov chain Monte Carlo algorithms for the Bayesian analysis of phylogenetic trees[END_REF] for a demonstration of the higher computational costs of bootstrapping over Bayesian approaches based on MCMC.

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.

We chose to draw the number of changes from a Poisson with mean 1 for each node.

When merging two time series (x1, . . . , xm) and (y1, . . . , yn), only the pairs xi → xj and yi → yj are presented to the DBN, while the pair xm → y1 is excluded due to the obvious discontinuity.